本發(fā)明涉及地球物理學中物探開發(fā)
技術領域：
，具體涉及基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合方法及系統(tǒng)。
背景技術：
：在油藏數(shù)值模擬中，為了使動態(tài)預測能夠盡量接近實際情況，通常需要對油藏數(shù)據進行歷史擬合，根據所觀測到的實際油藏動態(tài)來調整油藏模型參數(shù)，使得模型預測值與實際觀測值的誤差在允許范圍內，為后續(xù)油藏開采服務。傳統(tǒng)的歷史擬合方法通過手工來調整模型參數(shù)，工作量大、繁瑣，且效率低下。自動歷史擬合方法采用優(yōu)化算法自動調整油藏模型參數(shù)，大大縮短歷史擬合時間，提高擬合精度。因此，研究快速的自動歷史擬合方法是實現(xiàn)油藏歷史擬合的急切需求。歷史擬合問題是通過調整敏感參數(shù)(如孔隙度、滲透率等)，使得數(shù)值模擬計算的量如壓力、油氣比、含水率等都接近實際測量值，實質上是一個最優(yōu)化問題。在求解歷史擬合問題上，常見的有三種方法：梯度類方法、數(shù)據同化方法和隨機類方法。1、梯度類方法：梯度類算法中應用較多的是牛頓型方法。T.B.Tan和N.Kalogera設計了全隱式的三維三相的數(shù)值模擬器，將其應用在小型油藏模型中。但是Gauss-Newton方法不適合應用到大型的油藏模型中，因為Gauss-Newton方法在Hessian矩陣方面不便于計算。1995年孟雅杰在Gauss-Newton方法的基礎上提出了改進的牛頓法，這個方法并不需要計算Hessian矩陣。2002年Razza和Reynolds又對此進行了修正，加入有限儲存BFGS的策略，使得算法不再需要存儲矩陣，僅僅需要計算前一步的梯度值和目標值。該方法解決了Gauss-Newton方法不適合于處理大規(guī)模的油藏歷史擬合問題的弊端。2010年Razza和Reyonlds采用奇異值分解方法對算法參數(shù)進行降維，并將其應用到有限儲存策略中，以此提出新的自動歷史擬合的思路。梯度類算法是一種解決自動歷史擬合問題的有效算法，然而因為該方法對伴隨矩陣的計算的依賴性很高，并且其計算量大，不具備良好的可移植性。2、數(shù)據同化方法：集合卡爾曼濾波方法(ENKF)是一種十分重要的數(shù)據同化方法，該方法最早主要被使用在氣象學和海洋動力學中，因為集合卡爾曼濾波方法沒有利用梯度類算法中伴隨梯度的運算，所以在算法實現(xiàn)上較為方便，并且優(yōu)化后的油藏模型能夠體現(xiàn)真實油藏的不確定性。ENKF方法也存在同化循環(huán)中的濾波發(fā)散問題和在計算過程中不滿秩的問題。3、隨機類算法：隨機類算法是目前發(fā)展較快的一種算法，該類算法在計算過程中以隨機概率和搜索策略來求解問題，它能夠解決目標函數(shù)復雜和梯度求解困難的問題。2004年Tokuda和Takahashi將遺傳算法應用巖心驅替的歷史擬合中，實驗結果表明雖然遺傳算法能夠有效的求解歷史擬合問題，但是存在計算效率較低的問題，并且在歷史擬合中可能陷入局部收斂。雖然遺傳算法在計算過程中能夠搜索到較優(yōu)的解，但是在當油藏模型較大時計算效率較低。2009年YasinHajizadeh將ACO算法引入到歷史擬合問題的求解中，實驗結果表明該算法相對于傳統(tǒng)的遺傳算法求解效率更高，同年YasinHajizad將DE算法引入到歷史擬合問題的求解中，該算法僅需要少量的參數(shù)就能夠實現(xiàn)油藏自動歷史擬合，但是上述兩種算法在大型油藏模型中難以實現(xiàn)，并且存在如遺傳算法易陷入早熟收斂且計算速率緩慢，模擬退火計算量大等問題。另外，傳統(tǒng)方法往往采用高斯分布得到滲透率等模型參數(shù)初值，但由于油藏的強非均質性，特別是在多次注水,注化學藥劑的多次采油之后,儲層中各物性的不確定性強，模型參數(shù)特征一般呈現(xiàn)尖峰厚尾的非動力學特征，并不滿足高斯分布。技術實現(xiàn)要素：本發(fā)明所要解決的技術問題是提供基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合方法及系統(tǒng)，能夠采用t分布得到初始的油藏模型參數(shù)，并采用基于馬爾科夫鏈的蒙特卡洛法不斷優(yōu)化油藏模型參數(shù)擬合生產實際動態(tài)，得到盡可能接近真實模型的油藏數(shù)值模型，使油田開發(fā)動態(tài)預測的結果更加接近實際生產動態(tài)。本發(fā)明解決上述技術問題的技術方案如下：一方面，本發(fā)明提供了基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合方法，其特征在于，所述方法包括：S1、采用t分布隨機初始化得到初始的油藏靜態(tài)參數(shù)；S2、根據貝葉斯公式構造油藏模型的目標函數(shù)；S3、采用馬爾科夫鏈蒙特卡洛歷史擬合法對所述油藏靜態(tài)參數(shù)進行迭代優(yōu)化，得到優(yōu)化油藏靜態(tài)參數(shù)；S4、對迭代優(yōu)化得到的所有優(yōu)化油藏靜態(tài)參數(shù)對應的所有目標函數(shù)值進行最大后驗估計，得到最優(yōu)目標函數(shù)值，并輸出所述最優(yōu)目標函數(shù)值及其對應的最優(yōu)油藏靜態(tài)參數(shù)。本發(fā)明的有益效果：本發(fā)明提供的一種基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合方法，采用t分布得到初始的油藏靜態(tài)參數(shù)，然后采用馬爾科夫鏈蒙特卡洛歷史擬合法對所述油藏靜態(tài)參數(shù)進行迭代優(yōu)化，得到優(yōu)化油藏靜態(tài)參數(shù)，并對所有的優(yōu)化油藏靜態(tài)參數(shù)對應的所有目標函數(shù)值進行最大后驗估計，得到最優(yōu)目標函數(shù)值。本發(fā)明根據油藏的強非均質性，采用t分布得到初始的油藏靜態(tài)參數(shù)，符合油藏模型參數(shù)特征呈現(xiàn)尖峰厚尾的非動力學特征，基于概率統(tǒng)計基本思想，采用基于馬爾科夫鏈的蒙特卡洛法不斷優(yōu)化模型滲透率等油藏靜態(tài)參數(shù)擬合生產實際動態(tài)，對參數(shù)空間的不確定性進行量化，使預測值與真實值盡可能接近，得到盡可能接近真實模型的油藏數(shù)值模型，使油田開發(fā)動態(tài)預測的結果更加接近實際生產。本發(fā)明自動調整油藏模型參數(shù)，以縮短擬合時間，提高歷史擬合的效率和精度，研究對后期油藏開采方案的制定，以及后續(xù)生產過程優(yōu)化具有十分重要的意義。進一步的，所述S2具體包括：S21、根據貝葉斯公式得到油藏靜態(tài)參數(shù)的后驗分布函數(shù)的正比公式，所述油藏靜態(tài)參數(shù)的后驗分布函數(shù)正比于油藏靜態(tài)參數(shù)的先驗t分布的概率函數(shù)與油藏生產動態(tài)數(shù)據的正態(tài)分布的似然函數(shù)的乘積；S22、根據t分布公式和正態(tài)分布公式，得到所述油藏靜態(tài)參數(shù)的后驗分布函數(shù)的等式公式，具體由先驗項函數(shù)和似然項函數(shù)相加得到；S23、將所述油藏靜態(tài)參數(shù)的后驗分布函數(shù)的公式作為油藏模型的目標函數(shù)。采用上述進一步方案的有益效果：可以將求解油藏靜態(tài)參數(shù)的問題轉化為求解使目標函數(shù)的最小值，便于求解合適的油藏靜態(tài)參數(shù)。進一步的，所述S3具體包括：S31、設置馬爾科夫鏈的鏈長和優(yōu)化停止條件，將t分布得到的初始的油藏靜態(tài)參數(shù)作為當前狀態(tài)對應的優(yōu)化油藏靜態(tài)參數(shù)，并將所述當前狀態(tài)放入馬爾科夫鏈中；S32、根據當前狀態(tài)的對應的優(yōu)化油藏靜態(tài)參數(shù)，計算得到當前狀態(tài)對應的后驗分布函數(shù)值；S33、迭代產生下一個狀態(tài)對應的油藏靜態(tài)參數(shù)，并根據所述下一個狀態(tài)對應的油藏靜態(tài)參數(shù)，計算得到所述下一個狀態(tài)對應的后驗分布函數(shù)值；S34、根據所述下一個狀態(tài)對應的后驗分布函數(shù)值和當前狀態(tài)對應的后驗分布函數(shù)值，計算得到接受率R；S35、從0～1的均勻分布中隨機取一個數(shù)y，如果y≤R,則接受所述下一個狀態(tài)，并替代當前狀態(tài)成為新的當前狀態(tài)放入馬爾科夫鏈中，所述下一個狀態(tài)對應的油藏靜態(tài)參數(shù)成為新的當前狀態(tài)對應的優(yōu)化油藏靜態(tài)參數(shù)；否則不接受所述下一個狀態(tài)，依然將當前狀態(tài)放入馬爾科夫鏈中；S36、判斷是否滿足所述優(yōu)化停止條件，若滿足則結束流程，否則返回步驟S33。采用上述進一步方案的有益效果：采用馬爾科夫鏈蒙特卡洛歷史擬合法對所述油藏靜態(tài)參數(shù)進行迭代優(yōu)化，縮短了擬合時間，提升了擬合精度，克服了傳統(tǒng)隨機類方法運行計算量大的問題。進一步的，根據油藏靜態(tài)參數(shù)計算得到分布函數(shù)值具體包括：根據油藏靜態(tài)參數(shù)計算得到目標函數(shù)中先驗項函數(shù)對應的先驗項函數(shù)值；對所述油藏靜態(tài)參數(shù)采用油藏模擬器進行油藏擬合模擬計算，得到油藏生產動態(tài)數(shù)據；根據所述油藏靜態(tài)參數(shù)和所述油藏生產動態(tài)數(shù)據計算得到所述目標函數(shù)中似然項函數(shù)對應的似然項函數(shù)值；根據所述先驗項函數(shù)值和所述似然項函數(shù)值計算得到所述后驗分布函數(shù)值。采用上述進一步方案的有益效果：計算得到后驗分布函數(shù)值，以便用于計算得到接受率，并且所述后驗分布函數(shù)值即為目標函數(shù)值，以便對目標函數(shù)值進行最大后驗估計。另一方面，本發(fā)明提供了基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合系統(tǒng)，所述系統(tǒng)包括：初始化模塊，用于采用t分布隨機初始化得到初始的油藏靜態(tài)參數(shù)；構造模塊，用于根據貝葉斯公式構造油藏模型的目標函數(shù)；優(yōu)化模塊，用于采用馬爾科夫鏈蒙特卡洛歷史擬合法對所述油藏靜態(tài)參數(shù)進行迭代優(yōu)化，得到優(yōu)化油藏靜態(tài)參數(shù)；后驗估計模塊，用于對迭代優(yōu)化得到的所有優(yōu)化油藏靜態(tài)參數(shù)對應的所有目標函數(shù)值進行最大后驗估計，得到最優(yōu)目標函數(shù)值；輸出模塊，用于輸出所述最優(yōu)目標函數(shù)值及其對應的最優(yōu)油藏靜態(tài)參數(shù)。本發(fā)明的有益效果：本發(fā)明提供的一種基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合系統(tǒng)，采用t分布得到初始的油藏靜態(tài)參數(shù)，然后采用馬爾科夫鏈蒙特卡洛歷史擬合法對所述油藏靜態(tài)參數(shù)進行迭代優(yōu)化，得到優(yōu)化油藏靜態(tài)參數(shù)，并對所有的優(yōu)化油藏靜態(tài)參數(shù)對應的所有目標函數(shù)值進行最大后驗估計，得到最優(yōu)目標函數(shù)值。本發(fā)明根據油藏的強非均質性，采用t分布得到初始的油藏靜態(tài)參數(shù)，符合油藏模型參數(shù)特征呈現(xiàn)尖峰厚尾的非動力學特征，基于概率統(tǒng)計基本思想，采用基于馬爾科夫鏈的蒙特卡洛法不斷優(yōu)化模型滲透率等油藏靜態(tài)參數(shù)擬合生產實際動態(tài)，對參數(shù)空間的不確定性進行量化，使預測值與真實值盡可能接近，得到盡可能接近真實模型的油藏數(shù)值模型，使油田開發(fā)動態(tài)預測的結果更加接近實際生產。本發(fā)明自動調整油藏模型參數(shù)，以縮短擬合時間，提高歷史擬合的效率和精度，研究對后期油藏開采方案的制定，以及后續(xù)生產過程優(yōu)化具有十分重要的意義。進一步的，所述構造模塊具體包括：第一構造單元，用于根據貝葉斯公式得到油藏靜態(tài)參數(shù)的后驗分布函數(shù)的正比公式，所述油藏靜態(tài)參數(shù)的后驗分布函數(shù)正比于油藏靜態(tài)參數(shù)的先驗t分布的概率函數(shù)與油藏生產動態(tài)數(shù)據的正態(tài)分布的似然函數(shù)的乘積；第二構造單元，用于根據t分布公式和正態(tài)分布公式，得到所述油藏靜態(tài)參數(shù)的后驗分布函數(shù)的等式公式，具體由先驗項函數(shù)和似然項函數(shù)相加得到；目標函數(shù)構造單元，用于將所述油藏靜態(tài)參數(shù)的后驗分布函數(shù)的公式作為油藏模型的目標函數(shù)。采用上述進一步方案的有益效果：可以將求解油藏靜態(tài)參數(shù)的問題轉化為求解使目標函數(shù)的最小值，便于求解合適的油藏靜態(tài)參數(shù)。進一步的，所述優(yōu)化模塊具體包括：設置單元，用于設置馬爾科夫鏈的鏈長和優(yōu)化停止條件，將t分布得到的初始的油藏靜態(tài)參數(shù)作為當前狀態(tài)對應的優(yōu)化油藏靜態(tài)參數(shù)，并將所述當前狀態(tài)放入馬爾科夫鏈中；分布函數(shù)計算單元，用于根據當前狀態(tài)的對應的優(yōu)化油藏靜態(tài)參數(shù)，計算得到當前狀態(tài)對應的后驗分布函數(shù)值，以及用于根據下一個狀態(tài)對應的油藏靜態(tài)參數(shù)，計算得到下一個狀態(tài)對應的后驗分布函數(shù)值；迭代單元，用于迭代產生下一個狀態(tài)對應的油藏靜態(tài)參數(shù)；接受率計算單元，用于根據所述下一個狀態(tài)對應的后驗分布函數(shù)值和當前狀態(tài)對應的后驗分布函數(shù)值，計算得到接受率R；替換單元，用于從0～1的均勻分布中隨機取一個數(shù)y，如果y≤R,則接受所述下一個狀態(tài)，并替代當前狀態(tài)成為新的當前狀態(tài)放入馬爾科夫鏈中，所述下一個狀態(tài)對應的油藏靜態(tài)參數(shù)成為新的當前狀態(tài)對應的優(yōu)化油藏靜態(tài)參數(shù)；否則不接受所述下一個狀態(tài)，依然將當前狀態(tài)放入馬爾科夫鏈中；判斷單元，用于判斷是否滿足所述優(yōu)化停止條件，若滿足則結束流程，否則轉至所述迭代單元。采用上述進一步方案的有益效果：采用馬爾科夫鏈蒙特卡洛歷史擬合法對所述油藏靜態(tài)參數(shù)進行迭代優(yōu)化，縮短了擬合時間，提升了擬合精度，克服了傳統(tǒng)隨機類方法運行計算量大的問題。進一步的，所述分布函數(shù)計算單元具體包括：先驗項計算單元，用于根據油藏靜態(tài)參數(shù)計算得到目標函數(shù)中先驗項函數(shù)對應的先驗項函數(shù)值；油藏生產動態(tài)數(shù)據計算單元，對所述油藏靜態(tài)參數(shù)采用油藏模擬器進行油藏擬合模擬計算，得到油藏生產動態(tài)數(shù)據；似然項計算單元，用于根據所述油藏靜態(tài)參數(shù)和所述油藏生產動態(tài)數(shù)據計算得到所述目標函數(shù)中似然項函數(shù)對應的似然項函數(shù)值；函數(shù)值計算單元，用于根據所述先驗項函數(shù)值和所述似然項函數(shù)值計算得到所述后驗分布函數(shù)值。采用上述進一步方案的有益效果：計算得到后驗分布函數(shù)值，以便用于計算得到接受率，并且所述后驗分布函數(shù)值即為目標函數(shù)值，以便對目標函數(shù)值進行最大后驗估計。附圖說明圖1為本發(fā)明實施例1的基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合方法流程圖；圖2為本發(fā)明實施例1的概率密度函數(shù)隨自由度v的變化圖；圖3為本發(fā)明實施例1的基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合法詳細流程圖；圖4為本發(fā)明實施例1的PUNQS3模型的水平滲透率分布圖；圖5為本發(fā)明實施例1的井1的井底壓力、氣油比、含水率擬合圖；圖6為本發(fā)明實施例1的井4的井底壓力、氣油比、含水率擬合圖；圖7為本發(fā)明實施例1的井5的井底壓力、氣油比、含水率擬合圖；圖8為本發(fā)明實施例1的井11的井底壓力、氣油比、含水率擬合圖；圖9為本發(fā)明實施例1的井12的井底壓力、氣油比、含水率擬合圖；圖10為本發(fā)明實施例1的井15的井底壓力、氣油比、含水率擬合圖；圖11為本發(fā)明實施例1的井1的井底壓力、氣油比、含水率擬合圖；圖12為本發(fā)明實施例1的井4的井底壓力、氣油比、含水率擬合圖；圖13為本發(fā)明實施例1的井5的井底壓力、氣油比、含水率擬合圖；圖14為本發(fā)明實施例1的井11的井底壓力、氣油比、含水率擬合圖；圖15為本發(fā)明實施例1的井12的井底壓力、氣油比、含水率擬合圖；圖16為本發(fā)明實施例1的井15的井底壓力、氣油比、含水率擬合圖；圖17為本發(fā)明實施例2的基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合系統(tǒng)示意圖；圖18為本發(fā)明實施例2的基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合系統(tǒng)中構造模塊的結構示意圖；圖19為本發(fā)明實施例2的基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合系統(tǒng)中的優(yōu)化模塊的結構示意圖；圖20為本發(fā)明實施例2的基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合系統(tǒng)中的優(yōu)化模塊中的分布函數(shù)計算單元的結構示意圖。具體實施方式以下結合附圖對本發(fā)明的原理和特征進行描述，所舉實例只用于解釋本發(fā)明，并非用于限定本發(fā)明的范圍。實施例1、基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合方法。下面結合圖1至圖16對本實施例提供的方法進行詳細說明。參見圖1至圖3，S1、采用t分布隨機初始化得到初始的油藏靜態(tài)參數(shù)。具體的，油藏數(shù)值模型中需要優(yōu)化的參數(shù)是各劃分網格的油藏靜態(tài)參數(shù)如滲透率、孔隙度等，通過某種概率分布模型隨機賦初始值，傳統(tǒng)方法往往采用高斯分布得到滲透率等油藏靜態(tài)參數(shù)初始值，但由于油藏的強非均質性，特別是在多次注水、注化學藥劑的多次采油之后,儲層中各物性的不確定性強，模型參數(shù)特征一般呈σ現(xiàn)尖峰厚尾的非動力學特征，不滿足高斯分布。當實際單個變量的邊際分布比正常的邊際分布的尾部大的時候，可使用t分布來代替正態(tài)分布。t分布曲線的形狀與自由度v的大小相關，如圖2所示，自由度v越小，t分布曲線越平坦，曲線中間的值越低，曲線雙側尾部越高；自由度v越大，t分布曲線越接近正態(tài)分布曲線，當自由度v→∞時，t分布曲線逐步趨近于標準正態(tài)分布曲線。因而采用t分布隨機初始化得到初始的油藏靜態(tài)參數(shù)。所述油藏靜態(tài)參數(shù)可為個生產井的滲透率、孔隙度等參數(shù)，也可以為區(qū)塊各時間片的滲透率、孔隙度等參數(shù)。S2、根據貝葉斯公式構造油藏模型的目標函數(shù)。具體的，所述S2具體包括以下步驟：S21、根據貝葉斯公式得到油藏靜態(tài)參數(shù)的后驗分布函數(shù)的正比公式，所述油藏靜態(tài)參數(shù)的后驗分布函數(shù)正比于油藏靜態(tài)參數(shù)的先驗t分布的概率函數(shù)與油藏生產動態(tài)數(shù)據的正態(tài)分布的似然函數(shù)的乘積。具體的，傳統(tǒng)的貝葉斯方法應用在儲層數(shù)值模擬的時候，通過評估“最可能模型”來建立模型，其中先驗分布主要用于描述油藏靜態(tài)參數(shù)如孔隙度和滲透率等是否符合某種概率分布，后驗估計在抽樣以后可以得到，求解油藏靜態(tài)參數(shù)m的問題可以轉化為使目標函數(shù)O(m)取得最小值的問題。由貝葉斯公式得到油藏靜態(tài)參數(shù)m的后驗分布函數(shù)的正比公式，如公式(1)所示：p(m|dobs)∝p(dobs|m)·p(m)(1)其中，dobs為油藏生產動態(tài)數(shù)據，即為含水率、井底壓力和氣油比等參數(shù)；m為不確定的參數(shù)，即為滲透率等待優(yōu)化的油藏靜態(tài)參數(shù)；p(m)為m的先驗t分布的概率函數(shù)；p(dobs|m)為油藏生產動態(tài)數(shù)據的正態(tài)分布似然函數(shù)；p(m|dobs)為m的后驗分布函數(shù)。S22、根據t分布公式和正態(tài)分布公式，得到所述油藏靜態(tài)參數(shù)的后驗分布函數(shù)的等式公式，具體由先驗項函數(shù)和似然項函數(shù)相加得到。具體的，根據t分布的概率密度公式和正態(tài)分布的概率密度公式，得到所述油藏靜態(tài)參數(shù)m的后驗分布函數(shù)的具體的等式公式，其中，正態(tài)分布的似然函數(shù)具體為公式(2)所示：f(x)=1|Σ|(2π)dexp{-12{[g(m)-dobs]TCD-1[g(m)-dobs]}}---(2)]]>其中，d為向量m的維度；g(m)為符合不確定性參數(shù)的先驗概率分布；Σ為協(xié)方差矩陣,dobs為油藏生產動態(tài)數(shù)據。t分布的概率密度函數(shù)如公式(3)所示：f(x)=1|Σ|1/21(vπ)dΓ((v+d)/2)Γ(v/2)(1+m′Σ-1mv)-(v+d)/2---(3)]]>其中，x為向量，v為自由度,Σ為協(xié)方差矩陣,d為向量m的維度。因而，所述油藏靜態(tài)參數(shù)m的后驗分布函數(shù)的具體的等式公式如公式(4)所示：P(m|dobs)=1|Σ|1/21(vπ)dΓ((v+d)/2)Γ(v/2)(1+m′Σ-1mv)-(v+d)/2+1|Σ|(2π)dexp{-12{[g(m)-dobs]TCD-1[g(m)-dobs]}}---(4)]]>其中，v為自由度；d為向量m的維度；g(m)為符合不確定性參數(shù)的先驗概率分布；Σ為協(xié)方差矩陣。S23、將所述油藏靜態(tài)參數(shù)的后驗分布函數(shù)的公式作為油藏模型的目標函數(shù)。具體的，將所述油藏靜態(tài)參數(shù)的后驗分布函數(shù)的公式作為油藏模型的目標函數(shù)O(m)，具體如公式(5)所示：0(m)=1|Σ|1/21(vπ)dΓ((v+d)/2)Γ(v/2)(1+m′Σ-1mv)-(v+d)/2+1|Σ|(2π)dexp{-12{[g(m)-dobs]TCD-1[g(m)-dobs]}}]]>其中，μ為先驗值，為先驗項函數(shù)；為似然項函數(shù)。S3、采用馬爾科夫鏈蒙特卡洛歷史擬合法對所述油藏靜態(tài)參數(shù)進行迭代優(yōu)化，得到優(yōu)化油藏靜態(tài)參數(shù)。具體的，在連續(xù)的歷史擬合過程中，采用基于馬爾科夫鏈的蒙特卡洛法更新模型參數(shù)。其原理是采用先驗t分布隨機采樣得到油藏靜態(tài)參數(shù)的初始狀態(tài)，并循環(huán)進行狀態(tài)轉移，當隨機取的數(shù)小于等于接受率時，接收下一狀態(tài)，否則舍棄，將當前狀態(tài)放入鏈中；重復上述操作以得到優(yōu)化的油藏靜態(tài)參數(shù)。對于每一條馬爾科夫鏈需要考慮四個參數(shù)：(1)、馬爾科夫鏈的初始狀態(tài)initial，表征馬爾科夫鏈隨機取樣的起點；(2)、先驗項logprior，表征計算先驗項；(3)、后驗項loglikelihood，表征計算后驗項；(4)、馬爾科夫鏈鏈長mccount，表征馬爾科夫鏈轉移狀態(tài)的長度。基于先驗t分布隨機生成初始狀態(tài)，所述初始狀態(tài)對應t分布隨機初始化得到的油藏靜態(tài)參數(shù)的初始值，根據馬爾科夫鏈生成下一個狀態(tài)，計算得到接受率，從均勻分布中隨機取一個數(shù)，當隨機取的數(shù)小于等于接受率時，接收下一狀態(tài)，否則舍棄下一狀態(tài)，將當前狀態(tài)放入鏈中。從而進行不斷地循環(huán)，修改油藏參數(shù)，求解得到與生產歷史擬合匹配的最優(yōu)解。主要輸入的數(shù)據包括各類靜動態(tài)數(shù)據，如生產動態(tài)數(shù)據包括各油井的產液量、日產油、含水率等，網格數(shù)據，相對滲透率、毛管壓力以及油藏流體的PVT屬性數(shù)據，油、水、氣的地面密度，巖石壓縮系數(shù)等物性參數(shù)。具體的，所述S3具體包括以下步驟：S31、設置馬爾科夫鏈的鏈長和優(yōu)化停止條件，將t分布得到的初始的油藏靜態(tài)參數(shù)作為當前狀態(tài)對應的優(yōu)化油藏靜態(tài)參數(shù)，并將所述當前狀態(tài)放入馬爾科夫鏈中。所述優(yōu)化停止條件具體為到達馬爾科夫鏈的鏈長，也可根據情況設置其他停止條件。S32、根據當前狀態(tài)的對應的優(yōu)化油藏靜態(tài)參數(shù)，計算得到當前狀態(tài)對應的后驗分布函數(shù)值。S33、迭代產生下一個狀態(tài)對應的油藏靜態(tài)參數(shù)，并根據所述下一個狀態(tài)對應的油藏靜態(tài)參數(shù)，計算得到所述下一個狀態(tài)對應的后驗分布函數(shù)值。S34、根據所述下一個狀態(tài)對應的后驗分布函數(shù)值和當前狀態(tài)對應的后驗分布函數(shù)值，計算得到接受率R。具體的，根據所述下一個狀態(tài)對應的后驗分布函數(shù)值和當前狀態(tài)對應的后驗分布函數(shù)值，計算得到接受率R，所述接收率R的計算公式如公式(6)所示：R=P(mti+1|dobsti+1)P(mti|dobsti)---(6)]]>其中，所述為下一個狀態(tài)對應的后驗分布函數(shù)值，為當前狀態(tài)對應的后驗分布函數(shù)值。S35、從0～1的均勻分布中隨機取一個數(shù)y，如果y≤R,則接受所述下一個狀態(tài)，并替代當前狀態(tài)成為新的當前狀態(tài)放入馬爾科夫鏈中，所述下一個狀態(tài)對應的油藏靜態(tài)參數(shù)成為新的當前狀態(tài)對應的優(yōu)化油藏靜態(tài)參數(shù)；否則不接受所述下一個狀態(tài)，依然將當前狀態(tài)放入馬爾科夫鏈中。具體的，只有在一個狀態(tài)被接受后，該狀態(tài)對應的油藏靜態(tài)參數(shù)才會成為迭代得到的優(yōu)化油藏靜態(tài)參數(shù)；如果一個狀態(tài)沒有被接受，則該狀態(tài)對應的油藏靜態(tài)參數(shù)不會成為迭代得到的優(yōu)化油藏靜態(tài)參數(shù)，該油藏靜態(tài)參數(shù)會被舍棄。S36、判斷是否滿足所述優(yōu)化停止條件，若滿足則結束流程，否則返回步驟S33。具體的，根據油藏靜態(tài)參數(shù)計算得到分布函數(shù)值具體包括以下步驟：a、根據油藏靜態(tài)參數(shù)計算得到目標函數(shù)中先驗項函數(shù)對應的先驗項函數(shù)值。b、對所述油藏靜態(tài)參數(shù)進行油藏擬合模擬計算，得到油藏生產動態(tài)數(shù)據。具體為調用油藏數(shù)值模擬器進行油藏擬合模擬計算得到油藏生產動態(tài)數(shù)據。c、根據所述油藏靜態(tài)參數(shù)和所述油藏生產動態(tài)數(shù)據計算得到所述目標函數(shù)中似然項函數(shù)對應的似然項函數(shù)值。d、根據所述先驗項函數(shù)值和所述似然項函數(shù)值計算得到所述后驗分布函數(shù)值。所述后驗分布函數(shù)值與目標函數(shù)值相同。無論是當前狀態(tài)還是下一個狀態(tài)對應的后驗分布函數(shù)值均用上述方法進行計算。算法馬爾科夫鏈的蒙特卡洛法具體包括以下步驟：首先設置輸入輸出參數(shù)：輸入參數(shù)為：鏈長t為N(N正整數(shù))、先驗值μ、先驗項的計算函數(shù)、似然項的計算函數(shù)以及輸出結果設置；輸出參數(shù)為：馬爾科夫鏈。步驟1、初始化馬爾科夫鏈的初始狀態(tài)即鏈長t＝1時的初始狀態(tài)，即第1個狀態(tài)。步驟2、對鏈長t＝2，3，4,…,N；i＝1，2，3，…，N,循環(huán)以下過程進行采樣：2.1、從第i個狀態(tài)到第i+1個狀態(tài)，計算下一個狀態(tài)的值；2.2、計算接收率R；2.3、從0～1的均勻分布中隨機取一個數(shù)y，如果y≤R,則接受下一狀態(tài)，并放入馬爾科夫鏈中；否則又將當前狀態(tài)放入馬爾科夫鏈中?；趖分布的馬爾科夫鏈蒙特卡洛油藏參數(shù)自動歷史擬合法具體包括以下步驟：首先設置輸入輸出參數(shù)：輸入參數(shù)為：油藏模型數(shù)據文件、鏈長t為N(N正整數(shù))、先驗值μ、輸出結果位置和格式；輸出參數(shù)為：滲透率優(yōu)化參數(shù)和目標函數(shù)值。步驟1、基于t分布初始化馬爾科夫鏈的初始狀態(tài)并計算得到初始狀態(tài)對應的后驗分布函數(shù)值；步驟2、對t＝2，3，4，…，N；i＝1，2，3，…，N,直到達到馬爾科夫鏈的鏈長，循環(huán)以下過程進行滲透率的參數(shù)歷史擬合；2.1、產生可能的下一個狀態(tài)，并計算下一個狀態(tài)對應的先驗項函數(shù)值；2.2、運行油藏數(shù)值模擬程序ECLIPSE，計算得到油藏生產動態(tài)數(shù)據；2.3、計算所述下一個狀態(tài)對應的似然函數(shù)值；2.4、計算下一個狀態(tài)對應的后驗分布函數(shù)值；2.5、使用馬爾科夫鏈的蒙特卡洛法判斷是否跳轉到下一個狀態(tài)。步驟3、對目標函數(shù)進行最大后驗估計，得到滲透率分布。S4、對迭代優(yōu)化得到的所有優(yōu)化油藏靜態(tài)參數(shù)對應的所有目標函數(shù)值進行最大后驗估計，得到最優(yōu)目標函數(shù)值，并輸出所述最優(yōu)目標函數(shù)值及其對應的最優(yōu)油藏靜態(tài)參數(shù)。具體的，在采用所述馬爾科夫鏈蒙特卡洛法自動歷史擬合迭代更新結束后，會將在迭代更新過程中所有的優(yōu)化油藏靜態(tài)參數(shù)值以及其對應的所有的目標函數(shù)值進行輸出，所述目標函數(shù)值即所述后驗分布函數(shù)值，然后對所有的目標函數(shù)值進行最大后驗估計，得到最優(yōu)目標函數(shù)值，并輸出所述最優(yōu)目標函數(shù)值及其對應的最優(yōu)油藏靜態(tài)參數(shù)，所述最優(yōu)油藏靜態(tài)參數(shù)即滲透率分布。所有的優(yōu)化油藏靜態(tài)參數(shù)值具體指在迭代優(yōu)化過程中被接受的狀態(tài)的對應的油藏靜態(tài)參數(shù)，即被放入馬爾科夫鏈中的狀態(tài)對應的油藏靜態(tài)參數(shù)。另外，也可以在每次迭代更新得到一個優(yōu)化油藏靜態(tài)參數(shù)，就對所述優(yōu)化油藏靜態(tài)參數(shù)對應的目標函數(shù)值進行最大后驗估計，然后在迭代結束之后直接輸出得到的最優(yōu)目標函數(shù)值及其對應的最優(yōu)油藏靜態(tài)參數(shù)。綜上所述，采用基于t分布的馬爾科夫鏈蒙特卡洛油藏自動歷史擬合方法，通過調用油藏數(shù)值模擬軟件計算，使預測值與真實值盡可能接近，得到和真實油藏模型較為一致的數(shù)值模型。具體實例：1、主要通過對基于t分布的馬爾科夫鏈蒙特卡洛油藏參數(shù)自動歷史擬合法進行實驗以檢驗其效果。實驗采用的是PUNQ-S3油藏數(shù)據模型，PUNQ-S3油藏數(shù)據模型是一個三維三相的油藏工程模型，由19*28*25個網格塊構成，共分為五層，每層為2660個網格塊，每個網格塊的大小一致，其中包含1761個有效網格塊。如圖2所示，空白部分表示的是無效網格，線段特點的網格表示的是不同數(shù)值的水平滲透率，對于模型的每一層可以將水平滲透率分為不同的區(qū)塊，綜上所述，可將PUNQS3油藏模型的1761個網格，分為5*9共45個區(qū)塊，達到歷史擬合分區(qū)分塊的目的。PUNQS3模型每層的水平滲透率分布圖如圖4所示。2、油藏單井歷史擬合情況比較根據基于t分布的馬爾科夫鏈蒙特卡洛油藏自動歷史擬合法的實驗結果，對歷史擬合進行了相關的實驗和分析，其中鏈長設為500。目標函數(shù)值越小說明擬合值與實際測量值之間的差異程度越小，即擬合效果越好，效果越優(yōu)。為進一步比較說明基于正態(tài)分布和t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合方法的效果，將計算出的單口井的含水率(WWCT)、井底壓力(WBHP)和氣油比(WGOR)等參數(shù)同模型真實值分別進行對比，如圖5到圖10所示。其中，圖5-10分別是生產井1、4、5、11、12和15六口井的參數(shù)擬合圖，其中節(jié)點為圓圈的曲線為基于正態(tài)分布的馬爾科夫鏈的蒙特卡洛自動歷史擬合方法的擬合曲線，節(jié)點為三角的曲線為基于t分布的馬爾科夫鏈的蒙特卡洛自動歷史擬合方法的擬合曲線，節(jié)點為正方形的曲線為模型真實值。為進一步說明擬合效果，對擬合結果采用均方根誤差(RE)及整體誤差(EE)進行了計算和統(tǒng)計，計算公式如公式(7)和(8)所示：RE=1NΣi=1N[Di]2---(7)]]>EE=Σi=1N[|Ni-Ni′|]---(8)]]>其中，N為維度，Di＝Ni-N′i，Ni為真實值，N′i為擬合值?；谡龖B(tài)分布馬爾科夫鏈的蒙特卡洛自動歷史擬合方法的擬合誤差統(tǒng)計結果如表1所示，基于t分布馬爾科夫鏈的蒙特卡洛自動歷史擬合方法的擬合誤差統(tǒng)計結果如表2所示：表1基于正態(tài)分布馬爾科夫鏈的蒙特卡洛自動歷史擬合方法的誤差表表2基于t分布馬爾科夫鏈的蒙特卡洛自動歷史擬合方法的誤差表3、預測為進一步說明馬爾科夫鏈的蒙特卡洛自動歷史擬合方法生產預測的效果，取每口井前1460天的生產數(shù)據作為訓練集，訓練得到模型，然后采用訓練模型預測后1540天的生產數(shù)據并與真實數(shù)據進行對比，如圖11-16所示，其中，圖11-16分別是生產井1、4、5、11、12和15六口井的含水率(WWCT)、井底壓力(WBHP)和氣油比(WGOR)的預測結果。其中節(jié)點為圓圈的曲線為真實值，節(jié)點為三角的曲線為基于正態(tài)分布的馬爾科夫鏈的蒙特卡洛法的擬合曲線，節(jié)點為正方形的曲線為基于t分布的馬爾科夫鏈的蒙特卡洛法的擬合曲線。采用均方根誤差(RE)及整體誤差(EE)進一步計算和統(tǒng)計，基于正態(tài)分布馬爾科夫鏈的蒙特卡洛自動歷史擬合方法的擬合誤差統(tǒng)計結果如表3所示，基于t分布馬爾科夫鏈的蒙特卡洛自動歷史擬合方法的擬合誤差統(tǒng)計結果如表4所示：表3基于正態(tài)分布馬爾科夫鏈的蒙特卡洛自動歷史擬合方法預測誤差表表4基于t分布馬爾科夫鏈的蒙特卡洛自動歷史擬合方法預測誤差表綜上所述，基于t分布的馬爾科夫鏈蒙特卡洛油藏自動歷史擬合法與基于正態(tài)分布的馬爾科夫鏈蒙特卡洛自動歷史擬合法相比，大部分井參數(shù)預測值與真實值間的誤差更小，減小了油藏模型認識的不確定性，提高模型的預測能力，效果更優(yōu)。實施例2、基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合系統(tǒng)。下面結合圖17至圖20對本實施例提供的系統(tǒng)進行詳細說明。參見圖17至圖20，基于t分布的馬爾科夫鏈蒙特卡洛自動歷史擬合系統(tǒng)，其特征在于，所述系統(tǒng)包括初始化模塊、構造模塊、優(yōu)化模塊、后驗估計模塊以及輸出模塊。初始化模塊，用于采用t分布隨機初始化得到初始的油藏靜態(tài)參數(shù)。構造模塊，用于根據貝葉斯公式構造油藏模型的目標函數(shù)。具體的，所述構造模塊具體包括第一構造單元、第二構造單元以及目標函數(shù)構造單元。第一構造單元，用于根據貝葉斯公式得到油藏靜態(tài)參數(shù)的后驗分布函數(shù)的正比公式，所述油藏靜態(tài)參數(shù)的后驗分布函數(shù)正比于油藏靜態(tài)參數(shù)的先驗t分布的概率函數(shù)與油藏生產動態(tài)數(shù)據的正態(tài)分布的似然函數(shù)的乘積。第二構造單元，用于根據t分布公式和正態(tài)分布公式，得到所述油藏靜態(tài)參數(shù)的后驗分布函數(shù)的等式公式，具體由先驗項函數(shù)和似然項函數(shù)相加得到。目標函數(shù)構造單元，用于將所述油藏靜態(tài)參數(shù)的后驗分布函數(shù)的公式作為油藏模型的目標函數(shù)。優(yōu)化模塊，用于采用馬爾科夫鏈蒙特卡洛歷史擬合法對所述油藏靜態(tài)參數(shù)進行迭代優(yōu)化，得到優(yōu)化油藏靜態(tài)參數(shù)。所述優(yōu)化模塊具體包括設置單元、分布函數(shù)計算單元、迭代單元、接受率計算單元、替換單元以及判斷單元。設置單元，用于設置馬爾科夫鏈的鏈長和優(yōu)化停止條件，將t分布得到的初始的油藏靜態(tài)參數(shù)作為當前狀態(tài)對應的優(yōu)化油藏靜態(tài)參數(shù)，并將所述當前狀態(tài)放入馬爾科夫鏈中。分布函數(shù)計算單元，用于根據當前狀態(tài)的對應的優(yōu)化油藏靜態(tài)參數(shù)，計算得到當前狀態(tài)對應的后驗分布函數(shù)值，以及用于根據下一個狀態(tài)對應的油藏靜態(tài)參數(shù)，計算得到下一個狀態(tài)對應的后驗分布函數(shù)值。所述分布函數(shù)計算單元具體包括先驗項計算單元、油藏生產動態(tài)數(shù)據計算單元、似然項計算單元以及函數(shù)值計算單元。先驗項計算單元，用于根據油藏靜態(tài)參數(shù)計算得到目標函數(shù)中先驗項函數(shù)對應的先驗項函數(shù)值。油藏生產動態(tài)數(shù)據計算單元，對所述油藏靜態(tài)參數(shù)采用油藏模擬器進行油藏擬合模擬計算，得到油藏生產動態(tài)數(shù)據。似然項計算單元，用于根據所述油藏靜態(tài)參數(shù)和所述油藏生產動態(tài)數(shù)據計算得到所述目標函數(shù)中似然項函數(shù)對應的似然項函數(shù)值。函數(shù)值計算單元，用于根據所述先驗項函數(shù)值和所述似然項函數(shù)值計算得到所述后驗分布函數(shù)值。所述迭代單元，用于迭代產生下一個狀態(tài)對應的油藏靜態(tài)參數(shù)。所述接受率計算單元，用于根據所述下一個狀態(tài)對應的后驗分布函數(shù)值和當前狀態(tài)對應的后驗分布函數(shù)值，計算得到接受率R。所述替換單元，用于從0～1的均勻分布中隨機取一個數(shù)y，如果y≤R,則接受所述下一個狀態(tài)，并替代當前狀態(tài)成為新的當前狀態(tài)放入馬爾科夫鏈中，所述下一個狀態(tài)對應的油藏靜態(tài)參數(shù)成為新的當前狀態(tài)對應的優(yōu)化油藏靜態(tài)參數(shù)；否則不接受所述下一個狀態(tài)，依然將當前狀態(tài)放入馬爾科夫鏈中。具體的，若所述下一個狀態(tài)不被接受，該狀態(tài)對應的油藏靜態(tài)參數(shù)無法成為優(yōu)化油藏靜態(tài)參數(shù)，而是被直接舍棄掉。所述判斷單元，用于判斷是否滿足所述優(yōu)化停止條件，若滿足則結束流程，否則轉至所述迭代單元。后驗估計模塊，用于對迭代優(yōu)化得到的所有優(yōu)化油藏靜態(tài)參數(shù)對應的所有目標函數(shù)值進行最大后驗估計，得到最優(yōu)目標函數(shù)值。具體的，所述目標函數(shù)值即為計算得到的后驗分布函數(shù)值，二者的的值相同，另外，只有優(yōu)化油藏靜態(tài)參數(shù)對應的目標函數(shù)值才會進行最大后驗估計。輸出模塊，用于輸出所述最優(yōu)目標函數(shù)值及其對應的最優(yōu)油藏靜態(tài)參數(shù)。以上所述僅為本發(fā)明的較佳實施例，并不用以限制本發(fā)明，凡在本發(fā)明的精神和原則之內，所作的任何修改、等同替換、改進等，均應包含在本發(fā)明的保護范圍之內。當前第1頁1 2 3