基于馬爾科夫隨機場的巖性識別方法

田玉昆，周 輝，袁三一

中國石油大學油氣資源與探測國家重點實驗室，CNPC物探重點實驗室，北京 102249

1 引 言

近年來，中國陸上含油氣盆地勘探開發(fā)已經(jīng)日趨成熟，在勘探思路上已經(jīng)由構(gòu)造油氣藏向巖性油氣藏轉(zhuǎn)變.在這種條件下，準確識別巖性、流體，對油氣勘探開發(fā)具有重要意義.

自1984年Ostrander［1］提出了從地震資料中估算巖石彈性特征的AVO技術(shù)之后，通過地震資料進行巖性識別的方法得到了飛速發(fā)展.目前，通過地震數(shù)據(jù)反演得到彈性參數(shù)，再通過巖石物理模型的轉(zhuǎn)換得到含流體飽和度、孔隙度等儲層物性參數(shù)是進行巖性識別、儲層評價的常規(guī)手段.直接通過彈性參數(shù)進行巖性識別的方法主要有兩種：一是建立彈性參數(shù)量板，即統(tǒng)計各種巖性的彈性參數(shù)范圍，然后對反演得到的彈性參數(shù)進行判斷，落于某種巖性量板范圍內(nèi)巖性的彈性參數(shù)歸于該種巖性.由于不同巖性的彈性參數(shù)常常存在一定范圍的重疊，所以這種方法的識別精度不高，而且，由于地區(qū)間的差異，每次識別都要建立不同的量板，工作量大［2］；二是通過交會圖的方式，采用幾種對巖性敏感的彈性參數(shù)進行交會，選取在交會圖上能夠使巖性區(qū)分開的彈性參數(shù)或者巖石物理參數(shù)進行巖性的識別，這種方法通常是選取測井資料中已有的彈性參數(shù)進行交會，然后用交會的結(jié)果識別地震反演剖面上的巖性，同樣的，這種方法也會受限于彈性參數(shù)重疊的情況，如果通過交會圖，找不到合適的參數(shù)用于巖性識別，那么對反演剖面進行巖性識別就會陷入困境.此外，在交會圖進行區(qū)域劃分方面，一般都采取粗略地描述，或者手動勾繪，這種方法存在很大的不確定性及人為因素［3］.但是相對于直接通過地震剖面解釋的結(jié)果判別巖性，這種方法有一定的數(shù)據(jù)分析作為基礎，無疑要比僅定性分析的地震剖面識別方法更具有可信度.

基于概率統(tǒng)計理論進行巖性識別的方法在國外已趨成熟，這種方法通過統(tǒng)計巖石物理模型，從疊前地震數(shù)據(jù)反演得到儲層物性參數(shù)，進而進行巖性、流體等的預測.Mukerji等［4］給出了統(tǒng)計巖石物理方法應用流程，考慮實際數(shù)據(jù)誤差，并給出模型精確度描述.Mukerji等［4］和 Eidsvik等［5］研究了利用統(tǒng)計巖石物理模型預測儲層物性參數(shù)的方法，并給出了不確定性分析.Bachrach［6］基于 Gassmann方程建立統(tǒng)計巖石物理模型，計算孔隙度和含水飽和度用于巖性、流體預測.Larsen［7］、Gunning和 Glinsky［8］、Buland［9］、Eidsvik［10］以及 Marit等［11－14］使用疊前反演得到的彈性參數(shù)，經(jīng)由巖石物理模型，通過馬爾科夫鏈、馬爾科夫隨機場的條件概率形式作為巖性流體的垂向先驗分布，成功利用該方法進行巖性流體預測.鄧繼新和王尚旭［15］實現(xiàn)了基于統(tǒng)計巖石物理模型計算孔隙度和含水飽和度；胡華鋒［16］以地震資料和測井資料為反演基礎資料，統(tǒng)計巖石物理模型，反演得到孔隙度、含水飽和度和泥質(zhì)含量信息，利用這些信息來進行儲層的預測.

本文試圖直接利用疊前反演得到的彈性參數(shù)，通過對測井資料進行統(tǒng)計，運用概率統(tǒng)計的方法對巖性進行分析.這種方法在早期是通過巖性量板或者交會圖，建立起巖性同彈性參數(shù)之間的關系來進行巖性識別.但是由于彈性參數(shù)精度不高，并且在參數(shù)范圍上存在大量重疊，導致識別準確度不理想，同時，不論巖性量板法還是交會圖分析法都需要進行大量的準備工作，識別效率較低.在這種情況下，本文通過統(tǒng)計測井資料，建立起描述不同巖性的似然模型，然后用馬爾科夫隨機場建立先驗模型.在貝葉斯的框架下進行巖性識別，有利于綜合各種不同的先驗信息.該方法不僅僅限于本文提到的三種彈性參數(shù)，可以根據(jù)實際情況，采用對巖性敏感的彈性參數(shù)或者巖石物理參數(shù)約束巖性識別過程.本方法的優(yōu)勢在于，能夠通過馬爾科夫隨機場先驗模型，建立起點與點之間的相互關系作為空間約束，通過這種鄰域內(nèi)的相互關系，增強識別的準確性和空間上的連續(xù)性.該方法的意義在于，能夠克服參數(shù)范圍存在重疊情況下導致的識別不準問題，這就意味著，即使彈性參數(shù)反演的結(jié)果存在一定誤差，通過鄰域系統(tǒng)地調(diào)整，也能得到較準的識別結(jié)果.簡而言之，該方法不僅僅得到一個點的巖性，而且能夠描述相鄰位置的相互關系和相互作用，得到更符合實際情況的巖性剖面；同時，該方法加入了測井的指導，只要測井中出現(xiàn)的巖性，均能在剖面上識別.所以，通過對多口井的統(tǒng)計，能夠全面得到各種巖性的彈性參數(shù)分布范圍，從而準確得到巖性識別的結(jié)果.

2 方法原理

Marit和 Henning［12－13］介紹了利用馬爾科夫隨機場的條件概率形式作為先驗模型進行巖性/流體反演，受此啟發(fā)，結(jié)合圖像分割中的一些方法，本文利用反演的彈性參數(shù)作為特征場，以馬爾科夫模型作為先驗的巖性標號場對巖性識別進行了研究.考慮到地下儲層邊緣幾何特征，彈性參數(shù)在空間上的分布等因素，用馬爾科夫隨機場模型描述待識別的地下巖性的先驗分布符合沉積規(guī)律和地下構(gòu)造實際情況.馬爾科夫隨機場方法是建立在馬爾科夫隨機場模型和Bayesian理論的基礎上，馬爾科夫隨機場模型提供了不確定性描述與先驗知識聯(lián)系的紐帶，并利用疊前反演得到的彈性參數(shù)，根據(jù)統(tǒng)計決策和估計理論中的最優(yōu)準則確定巖性分類問題的目標函數(shù)，求解滿足這些條件的最大可能分布，從而將分類問題轉(zhuǎn)化為最優(yōu)化問題.

馬爾科夫隨機場模型在儲層屬性隨機模擬中應用已久，它利用變量狀態(tài)的轉(zhuǎn)移概率矩陣對隨機過程的未來狀態(tài)進行分析預測，是一種半定量的研究方法［17－20］.簡而言之，某點的性質(zhì)，只與其鄰域系統(tǒng)有關.基于馬爾科夫隨機場的地質(zhì)統(tǒng)計模型可以很好地反映復雜空間的連續(xù)性，顯然這種結(jié)果對巖性識別是有利的.

2.1 馬爾科夫隨機場

馬爾科夫隨機場是對一維馬爾科夫過程的擴展.對于二維平面上的巖性剖面，可以將其視為一個馬爾科夫隨機場.

圖1 不同階鄰域系統(tǒng)及其勢團Fig.1 Different order neighborhood system and its clique

設關于巖性π的隨機場是關于鄰域系統(tǒng)δ的馬爾科夫隨機場，ω為關于巖性π的分類標號.馬爾科夫隨機場的性質(zhì)可以描述為

式（2）稱為馬爾科夫隨機場的條件概率形式，事實上，它很難求取、Marit和Henning利用從測井資料統(tǒng)計得到的計數(shù)矩陣對其進行了描述，反映了巖性剖面的局部性質(zhì)，但是這種方法無法對整個剖面進行定義.Hammersley－Clifford定理揭示了馬爾科夫隨機場與Gibbs分布之間的等價性，把馬爾科夫隨機場與Gibbs分布對應起來，通過單個位置點及其鄰域的簡單特性得到整個剖面全局特性.

Hammersley－Clifford定理：設S是一個鄰域系統(tǒng)，當且僅當p（π）是一個關于S的Gibbs分布時，π是關于S的馬爾科夫隨機場.Gibbs的分布形式為［21］

其中，Z是歸一化常數(shù)，T表示溫度，U（ω）表示能量函數(shù).

2.2 MRF—MAP框架下的巖性識別

本文巖性識別的方法是在Bayesian框架下，利用馬爾科夫隨機場對巖性的先驗信息進行描述，反演目標為巖性類別，已知參數(shù)為疊前反演得到的彈性參數(shù).這里 MAP是指最大后驗概率，它是Bayesian框架下反演最常用的最優(yōu)化準則.

如之前定義，S為定義在巖性剖面上的空間位置集合，假設反演得到的彈性參數(shù)場為R，那么R的任意一次實現(xiàn)可以記為r＝ （r1，r2，…，rM×N），巖性反演問題需要求解的是滿足最大后驗概率準則條件下對于每一個位置點巖性的分類，這個分類過程通過對巖性進行標號完成，將標號結(jié)果稱為標號場，記為ωk，表示共有k個類別.根據(jù)貝葉斯公式［22］，有

其中p（r）為彈性參數(shù)場的先驗分布，通常情況下是一個常數(shù)，不參與計算，一般不予考慮.p（ω）是巖性標號場的先驗Gibbs分布，滿足馬爾科夫性.在本研究中，采用馬爾科夫隨機場的potts模型［23］，其能量函數(shù)具體表達式為

其中，δ（ωi，ωj）為 Kronecker delta，U（ω）為式（3）中的能量函數(shù)，通過其能量來描述某點屬于某種巖性的概率，最大概率的狀態(tài)對應于鄰域點為同一巖性的狀態(tài)，通過這種描述，使得在一定范圍內(nèi)，相近的彈性參數(shù)被劃歸為同一種巖性，保證空間上的連續(xù)性.p（r｜ω）是似然函數(shù)，一般情況下，假定彈性參數(shù)之間彼此獨立，即

其中，n為參與計算的彈性參數(shù)的種類數(shù)目.假設p（r｜ω）服從高斯分布，即每個類都可以由彈性參數(shù)的期望μi和標準差σi來唯一確定其高斯分布：

給定巖性類別的似然函數(shù)與先驗函數(shù)以后，通過先驗函數(shù)能夠保證空間上的連續(xù)性，同時，通過似然函數(shù)的調(diào)整，能夠使每一點的巖性在多種彈性參數(shù)共同作用的情況下盡可能準確地分類.

綜上所述，在MAP準則下求解分類結(jié)果，可以歸結(jié)為求使式（4）最大時候ω的估計，記為

整個建模過程具體如下：

（1）根據(jù)問題的性質(zhì)，定義S上的鄰域系統(tǒng)，選擇合適的階數(shù)得到相應的勢團；

（2）定義勢函數(shù)的具體形式并寫出Gibbs分布；

（3）定義似然函數(shù)，選擇合適的分布形式；

（4）寫出后驗概率最簡化的表示形式；

（5）估計模型參數(shù)并優(yōu)化尋找MAP準則下的解.

2.3 算法流程

針對本文研究的巖性反演目標函數(shù)，采用條件迭代（Iterated Conditional Mode簡寫為ICM）的算法求解該反問題.假定通過疊前反演得到的彈性參數(shù)場為R，巖性標號場為ω，具體算法流程如下：

（1）通過統(tǒng)計測井資料，得到不同巖性下參數(shù)與巖性反演的各彈性參數(shù)的均值μi與標準差σi；

（2）根據(jù)似然函數(shù)最大準則求取最大似然分類，以此作為初始分類0，即對剖面上的每一點取其最大似然概率下的標號，遍歷整個剖面得到初始分類；

（4）比較第k次與第k－1次的分類結(jié)果，如果滿足迭代終止條件則將第k次結(jié)果作為最終分類結(jié)果，否則取k＝k＋1，轉(zhuǎn)到第3步繼續(xù)計算，直到滿足迭代終止條件為止.

將上述過程表示成流程圖的形式如圖2所示.

圖2 ICM算法流程圖Fig.2 ICM algorithm flow chart

3 模型試算

為了測試上述方法，本文設計了兩個模型進行檢驗：（1）水平和楔形地層模型，用以檢驗本文方法在彈性參數(shù)范圍有重疊情況下識別巖性的能力和描述薄層的能力；（2）Marmousi II模型，用以檢驗本文方法描述地下復雜構(gòu)造空間形態(tài)的能力.

3.1 水平和楔形地層模型

圖3為一個垂向上為500個樣點，橫向上為100道的巖性剖面，整個剖面上一共分為4種巖性.巖性為塊狀分布，在水平層位置，沒有橫向上彈性參數(shù)的變化.在垂直位置101～200的范圍內(nèi)，有一楔形模型，用來驗證對薄層形狀的刻畫能力.

假設在第50道有一口井，并有縱波速度、橫波速度和密度曲線.通過自上而下統(tǒng)計各巖性的彈性參數(shù)范圍，求取高斯分布的均值和標準差.各巖性的彈性參數(shù)范圍如圖4所示：

通過圖4可以看到，不同巖性彈性參數(shù)范圍之間存在一定的重疊，這樣，僅僅通過觀察彈性參數(shù)剖面很難區(qū)分出不同的巖性.在實際資料中，也存在這樣的情況，所以，通過彈性參數(shù)直接識別巖性存在一定的困難.

圖5是用本文方法得到的識別結(jié)果，從圖中可以看到，楔形模型得到了完美的識別，表明本文方法能很好地描述地下的構(gòu)造形態(tài)，并盡可能地保留反演結(jié)果的分辨率.這對于刻畫砂體展布、河道形態(tài)和儲層范圍具有重要意義.

但是在圖5中也可以看到，在垂向坐標值51～55中間有一個薄層未被識別，這是因為在測井資料統(tǒng)計的過程中，是按照點數(shù)統(tǒng)計，該巖性的這個薄層所占點數(shù)很少，所以對于最終的概率分布函數(shù)的形態(tài)影響不大，另外，該部分與其上覆的巖性在彈性參數(shù)數(shù)值上有較大的重疊，所以在識別結(jié)果中歸為上覆的地層.在實際中，這種情況并不多見，對于一種巖性，通過統(tǒng)計其彈性參數(shù)得到的概率分布函數(shù)，能夠全面均勻地包含各點信息，因此，能夠準確地描述各點屬于該種巖性的概率.

圖5中，除去垂向坐標值51～55中的一層外，每種巖性都得到了準確分類，表明本文提出的巖性識別方法合理有效.當?shù)貙訌椥詤?shù)在某范圍內(nèi)大致均勻分布時，即使與其它巖性的彈性參數(shù)有較大的重疊，在多參數(shù)聯(lián)合參與識別以及先驗函數(shù)共同調(diào)整下，這些巖性也能得到很好地識別.

3.2 Marmousi II模型

為了進一步檢驗該方法對于復雜地下構(gòu)造形態(tài)的識別能力，以下采用Marmousi II模型進行檢驗.

Marmousi II模型是一個垂向、橫向上彈性參數(shù)變化比較大的模型.在如圖6所示橫向上第100道、第700道的位置處，假設有兩口井，對100道位置的井根據(jù)其彈性參數(shù)范圍，人為的分為5種巖性，并使用同樣的方法流程統(tǒng)計彈性參數(shù)范圍和計算未知參數(shù).700道位置處的井不參與統(tǒng)計，用于后面同識別結(jié)果進行對比.最終得到基于兩種不同鄰域系統(tǒng)的分類結(jié)果（如圖6）.圖6a為基于二階鄰域系統(tǒng)的分類結(jié)果，圖6b為基于一階鄰域系統(tǒng)的分類結(jié)果.

從圖6a中可以發(fā)現(xiàn)，整個剖面上分類層次清楚，5種巖性都得到了準確地識別，結(jié)果的形態(tài)符合Marmousi II模型.每種巖性層內(nèi)平滑，不存在分類錯誤的噪聲點.垂向上的薄層得到了較好辨別，橫向上的結(jié)果比較連續(xù)，并且圖中背斜、斷層的形態(tài)都得到了準確地刻畫，這對于油氣的勘探是很有利的.表明本文方法在地下彈性參數(shù)連續(xù)且有較大重疊的情況下，也能比較準確地獲得理想的分類結(jié)果.

圖6b是一階鄰域系統(tǒng)對應的分類結(jié)果，與二階鄰域系統(tǒng)相比較差別不大，主要區(qū)別是在空間上的連續(xù)性不如二階.

圖7為橫向第100道、700道位置處的真實巖性與識別巖性結(jié)果的對比.在兩個位置，兩條曲線都基本重合，在巖性發(fā)生劇烈變化的位置處，識別結(jié)果有一些偏差，但是這種偏差很小，不影響最終整個剖面上的分類結(jié)果.這表明本文巖性識別的方法精度還是很高的.

圖8描述該方法的收斂速度，取誤差的對數(shù)顯示.從圖中可以看到，本文方法收斂很快，經(jīng)過10次迭代以后，基本得到收斂，因此在計算效率上是可以接受的.

為了討論反演處理不確定性對該方法的影響，在Marmousi II模型中，分別加入了1%和5%的隨機誤差（5%的誤差分布如圖9所示），以檢驗處理過程中產(chǎn)生的不確定性對巖性識別的影響.

圖8 收斂速度Fig.8 Convergence rate

圖10為加入1%和5%的誤差后巖性識別的結(jié)果.從圖中可以看出，隨著誤差的增大，各種巖性大致上還是比較清楚地區(qū)分開，同時，巖層的空間結(jié)構(gòu)也得到較好地刻畫.這表明該方法在反演結(jié)果不太準確的條件下，也有比較準確的識別結(jié)果.但同時，對比圖10a和圖10b可以發(fā)現(xiàn)，隨著誤差的增大，巖性識別的結(jié)果會受到一定影響，一些巖性層變得不連續(xù)，產(chǎn)生模糊的邊界.這表明反演的不確定性會對巖性識別產(chǎn)生一定程度的影響，并且隨著誤差比例的增加，對巖性識別結(jié)果的影響也會增強.

4 結(jié) 論

本文提出了一種利用反演的彈性參數(shù)識別巖性的方法.該方法主要的優(yōu)勢在于，對于有一定范圍重疊或者誤差的彈性參數(shù)，本方法能夠在最大程度上避免識別錯誤，且高效快捷.該方法的巖性識別主要通過馬爾科夫隨機場鄰域系統(tǒng)的調(diào)整來實現(xiàn).事實上馬爾科夫隨機場本身并不是一種巖性識別方法，但它可以作為一個先驗模型嵌入到某種特定的算法中，以達到影響巖性識別結(jié)果的目的.具體而言，其通過調(diào)整相鄰點的相互關系和相互作用，從全局角度描述巖性，將各點的相互關系加以傳播，從而得到空間上有一定連續(xù)性的識別結(jié)果.同時，該方法的識別過程不損害反演結(jié)果已有的分辨率，能夠最有效地保護反演分辨率.

從本文中的例子可以看到，馬爾科夫隨機場能夠很好的描述剖面上彈性參數(shù)表示出的輪廓和結(jié)構(gòu)，空間上有一定的連續(xù)性，能夠細致刻畫地下構(gòu)造空間上的形態(tài)，得到合適的分類結(jié)果.馬爾科夫隨機場模型在反映巖性剖面隨機性的同時，又能反映剖面的一些潛在固有的空間結(jié)構(gòu)，可以更有效地描述地下構(gòu)造.本方法采用縱橫波速度和密度參數(shù)進行巖性識別，事實上，在具體使用中，并不限于這三種彈性參數(shù)，采用對巖性更為敏感的參數(shù)，能夠使本方法的識別結(jié)果更精確.此外，該方法也存在一定的局限性，馬爾科夫隨機場先驗模型產(chǎn)生的是一個塊的識別結(jié)果，如果識別的地層需要強調(diào)一種連續(xù)的變化，該方法會有一定的局限.

（References）

［1］ Ostrander W J.Plane－wave reflection coefficients for gas sands at nonnormal angles of incidence.Geophysics，1984，49（10）：1637－1648.

［2］ 何又雄，姚姚.基于參考道的巖性識別與巖性剖面非線性反演.石油勘探與開發(fā)，2005，32（3）：61－63.He Y X，Yao Y.Lithological profile non－linear inversion by reference trace－based lithological identification.Petroleum Exploration and Development （in Chinese），2005，32（3）：61－63.

［3］ 李國福.多參數(shù)儲層預測及流體識別方法研究［碩士論文］.成都：成都理工大學，2011.Li G F. Multi－parameter reservoir prediction and fluid identification method research（in Chinese）［Master′s thesis］.Chengdu：Chengdu University of Technology，2011.

［4］ Mukerji T，Avseth T，Mavko G，et al.Statistical rock physics：Combining rock physics，information theory，and geostatistics to reduce uncertainty in seismic reservoir characterization.The Leading Edge，2001，20（3）：313－319.

［5］ Eidsvik J，Avseth P，Omre H，et al.Stochastic reservoir characterization using prestack seismic data.Geophysics，2004，69（4）：978－993.

［6］ Bachrach R.Joint estimation of porosity and saturation using stochastic rock－physics modeling.Geophysics，2006，71（5）：O53－O63.

［7］ Larsen A L，Ulvmoen M，Omre H，et al.Bayesian lithology/fluid prediction and simulation on the basis of a Markov－chain prior model.Geophysics，2006，71（5）：R69－R78.

［8］ Gunning J，Glinsky M E.Detection of reservoir quality using Bayesian seismic inversion.Geophysics，2007，72（3）：R37－R49.

［9］ Buland A，Kolbj?rnsen O， Hauge R，et al.Bayesian lithology and fluid prediction from seismic prestack data.Geophysics，2008，73（3）：C13－C21.

［10］ Eidsvik J，Omre H，Mukerji T，et al.Seismic reservoir prediction using Bayesian integration of rock physics and markov random fields：A North Sea example.The Leading Edge，2002，21（3）：290－294.

［11］ Spikes K，Mukerji T，Dvorkin J，et al.Probabilistic seismic inversion based on rock－physics models.Geophysics，2007，72（5）：R87－R97.

［12］ Ulvmoen M，More H.Improved resolution in Bayesian lithology/fluid inversion from prestack seismic data and well observations：Part 1－Methodology.Geophysics，2010，75（2）：R21－R35.

［13］ Ulvmoen M，More H，Buland A.Improved resolution in Bayesian lithology/fluid inversion from prestack seismic data and well observations：Part 2－Real case study.Geophysics，2010，75（2）：B73－B82.

［14］ Ulvmoen M，Hammer H.Bayesian lithology/fluid inversioncomparison of two algorithms.Computational Geosciences，2010，14（2）：357－367.

［15］ 鄧繼新，王尚旭.基于統(tǒng)計巖石物理的含氣儲層飽和度與孔隙度聯(lián)合反演.石油天然氣學報，2009，31（1）：48－53.Deng J X，Wang S X.Joint inversion of saturation and porosity in gas reservoirs based on statistical rock physics.Journal of Oil and Gas Technology （in Chinese），2009，31（1）：48－53.

［16］ 胡華鋒.基于疊前道集的儲層參數(shù)反演方法研究［碩士論文］.青島：中國石油大學（華東），2009.Hu H F.The Study of Petrophysical－Properties Inversion Base on Pre－stack Seismic Data （in Chinese）［Master′s thesis］.Qingdao：China University of Petroleum，2009.

［17］ Elfeki A，Dekking M.A Markov chain model for subsurface characterization：Theory and applications. Mathematical Geology，2001，33（5）：569－589.

［18］ Carle S F，F(xiàn)ogg G E.Modeling spatial variability with one and multidimensional continuous－Lag Markov chains.Mathematical Geology，1997，29（7）：891－918.

［19］ Weissmann G S， Fogg G E. Multi－scale alluvial fan heterogeneity modeled with transition probability geostatistics in a sequence stratigraphic framework.Journal of Hydrology，1999，226（1－2）：48－56.

［20］ Norberg T，Rosén L，Baran A，et al.On modeling discrete geological structures as Markov random fields.Mathematical Geology，2002，34（1）：63－77.

［21］ 李旭超，朱善安.圖像分割中的馬爾可夫隨機場方法綜述.中國圖像圖形學報，2007，12（5）：789－798.Li X C，Zhu S A.A survey of the Markov random field method for image segmentation.Journal of Image and Graphics （in Chinese），2007，12（5）：789－798.

［22］ Yuan S Y，Wang S X，Li G F.Random noise reduction using Bayesian inversion.Journal of Geophysics and Engineering，2012，9（1）：60－68.

［23］ Pérez P.Markov random fields and images.CWI Quarterly，1998，11（4）：413－437.