基于分形高頻初始模型和低頻先驗(yàn)信息的物性參數(shù)隨機(jī)反演

印興耀，孫瑞瑩，張廣智，王保麗

(中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島266580)

在勘探地球物理中，大尺度(低頻)成分和小尺度(高頻)成分都是儲層和油藏描述的重要信息。地震數(shù)據(jù)的橫向連續(xù)性好，但是具有帶限性質(zhì)；測井?dāng)?shù)據(jù)可提供鉆孔附近的高分辨率的彈性參數(shù)和儲層物性參數(shù)，但是缺乏區(qū)域覆蓋。貝葉斯理論框架下的隨機(jī)反演方法可以更好地融合地震數(shù)據(jù)和測井信息。常規(guī)的隨機(jī)反演通常假定先驗(yàn)概率密度服從高斯分布，再利用測井?dāng)?shù)據(jù)建立變差函數(shù)來約束零空間。然而，服從高斯分布的假設(shè)有時(shí)是不恰當(dāng)?shù)?，而且高斯分布僅由均值和方差表示。分形高斯概率密度函數(shù)由3個(gè)參數(shù)描述，即均值、方差和Hurst系數(shù)[1]。該方法可以產(chǎn)生高頻的先驗(yàn)信息，并且基于分形理論產(chǎn)生的初始模型可以減少超出測井頻帶的假頻出現(xiàn)[2]?；诜中卫碚摰南闰?yàn)信息主要表征小尺度變化，大尺度的趨勢可以由測井?dāng)?shù)據(jù)進(jìn)行克里金插值并通過低通濾波得到，該初始模型常被用于確定性反演中，以獲取地層趨勢。這樣，混合先驗(yàn)信息就包含了高頻成分和低頻成分。

儲層物性參數(shù)是進(jìn)行儲層預(yù)測與流體識別的重要依據(jù)?；诘卣鸱囱莸玫降膹椥詤?shù)，通過巖石物理模型估計(jì)儲層物性參數(shù)，會存在不確定性傳遞，因此有必要開展物性參數(shù)的直接反演。而巖石物理是建立地震彈性參數(shù)和儲層物性參數(shù)之間關(guān)系的橋梁，運(yùn)用巖石物理理論可以進(jìn)行巖石骨架和孔隙流體性質(zhì)定性評價(jià)。Mukerji等[3-4]為了利用疊前AVO資料預(yù)測儲層參數(shù)并對相應(yīng)的不確定性進(jìn)行評價(jià)，引入了統(tǒng)計(jì)巖石物理模型。隨后Eidsvik等[5]、Spikes等[6]、Grana等[7-8]、胡華鋒等[9]和夏麗娜等[10]也都進(jìn)行了儲層物性參數(shù)的反演研究。我們采用貝葉斯理論進(jìn)行反演以獲取后驗(yàn)概率密度的樣本[11]，通過貝葉斯理論框架結(jié)合分形高斯模型算法、克里金插值、SA-PSO優(yōu)化算法和統(tǒng)計(jì)巖石物理理論，克服了地震資料的帶限性質(zhì)，可以得到寬頻帶的反演結(jié)果。

1 方法原理

首先，基于分形理論得到高頻初始模型，利用克里金插值得到低頻初始模型，兩者混合就是本文的初始先驗(yàn)實(shí)現(xiàn)；然后，利用SA-PSO優(yōu)化算法進(jìn)行優(yōu)化，得到后驗(yàn)概率密度的樣本。其中，優(yōu)化算法中目標(biāo)函數(shù)的建立包含了統(tǒng)計(jì)巖石物理理論。圖1為基于分形高頻初始模型和低頻先驗(yàn)信息的儲層物性參數(shù)隨機(jī)反演流程。

圖1 基于分形高頻初始模型和低頻先驗(yàn)信息的儲層物性參數(shù)隨機(jī)反演流程

1.1 巖石物理模型

巖石物理模型可以是測井曲線的簡單回歸擬合，也可以是更為復(fù)雜的物理模型。測井曲線的簡單回歸擬合一般不適用于整個(gè)工區(qū)，我們引入微分等效介質(zhì)(DEM)模型和Gassmann方程，采用變化的孔隙縱橫比逐點(diǎn)擬合測井曲線的方式來建立確定性巖石物理模型。首先，利用DEM解析模型和Gassmann方程建立巖石的縱、橫波速度，密度與孔隙度,飽和度和礦物組分等各參數(shù)之間的關(guān)系;其次，將巖石孔隙等效為具有單一縱橫比的理想橢球孔，應(yīng)用非線性全局尋優(yōu)算法來尋找最佳的等效孔隙縱橫比，使得理論預(yù)測與實(shí)際測量的彈性模量之間的誤差最小;最后，將反演得到的等效孔隙縱橫比代入到Gassmann方程和DEM模型中構(gòu)建縱波速度、橫波速度和密度[12]。但是常規(guī)的確定性巖石物理模型都有一定的適用范圍和局限性，因此我們在此基礎(chǔ)上引入統(tǒng)計(jì)巖石物理模型。

統(tǒng)計(jì)巖石物理模型由確定性巖石物理模型和隨機(jī)誤差項(xiàng)構(gòu)成，其中確定性巖石物理模型表征儲層物性參數(shù)與彈性參數(shù)之間的巖石物理關(guān)系，隨機(jī)誤差項(xiàng)則用于彌補(bǔ)確定性巖石物理模型的局限性。統(tǒng)計(jì)巖石物理模型的數(shù)學(xué)表達(dá)式為

(1)

式中：m表示彈性參數(shù)體，m=[vP,vS,ρ];R代表儲層物性參數(shù)，R=[φ,Vsh,Sw];vP,vS,ρ,φ,Vsh,Sw分別表示縱波速度、橫波速度、密度、孔隙度、泥質(zhì)含量和含水飽和度;fRPM(R)表示確定性巖石物理模型;ε為誤差項(xiàng)，可以通過確定性巖石物理模型與實(shí)際測井資料之間的相對差異估算，通常取零均值截?cái)喔咚拐`差。

1.2 貝葉斯反演方法

AVO反演的理論基礎(chǔ)是Zoeppritz方程及其近似式。我們提出的非線性反演方法在計(jì)算反射系數(shù)時(shí)，直接采用精確的Zoeppritz方程計(jì)算縱波反射系數(shù)，表達(dá)式為

(2)

式中：RPP表示P波反射系數(shù);RPS表示S波反射系數(shù);TPP表示P波透射系數(shù);TPS表示S波透射系數(shù);θ1和θ2分別為P波反射角和透射角;φ1和φ2分別為S波反射角和透射角;Δα，Δβ和Δρ分別表示界面兩側(cè)的縱、橫波速度差及密度差，Δα=α2-α1，Δβ=β2-β1，Δρ=ρ2-ρ1;α，β和ρ分別表示縱波速度、橫波速度及密度;下標(biāo)1和2分別表示界面上、下兩側(cè)的參數(shù)信息。

假設(shè)地震褶積模型為

(3)

式中：d為觀測地震數(shù)據(jù);r為反射系數(shù)序列;e為服從高斯分布的隨機(jī)噪聲;G為子波褶積矩陣?？梢越⒎瓷湎禂?shù)序列r和反演的儲層物性參數(shù)R之間的關(guān)系，即得到

(4)

根據(jù)貝葉斯理論可知[13]

(5)

式中：ρM(R)表示物性參數(shù)的先驗(yàn)信息(基于分形理論和克里金插值建立的先驗(yàn)?zāi)Ｐ?;L(d/R)代表似然函數(shù)，表示模型與數(shù)據(jù)的匹配程度;k是概率歸一化常數(shù);σM(R)是貝葉斯后驗(yàn)分布。

1.3 先驗(yàn)信息建立

1.3.1 基于分形高斯模型的高頻先驗(yàn)信息

測井資料分析表明，測井曲線的功率譜、變差函數(shù)和協(xié)方差通常呈冪律分布，冪指數(shù)為Hurst系數(shù)[14]。這樣可以應(yīng)用均值為0，具有特定自協(xié)方差結(jié)構(gòu)的分形高斯分布來描述測井曲線形態(tài)[1]。自協(xié)方差函數(shù)為

(6)

式中：R是自協(xié)方差;σ是測井曲線的標(biāo)準(zhǔn)差;t是時(shí)間間隔;H是通過R/S分析方法求取的孔隙度、泥質(zhì)含量和含水飽和度的Hurst系數(shù)[15]。

我們采用Srivastava等[1]的方法生成基于分形理論的高頻初始模型。為了生成分形高斯模型，首先計(jì)算給定測井曲線的均值μ，方差σ2和Hurst系數(shù)H。然后利用公式(6)計(jì)算自協(xié)方差。根據(jù)維納-辛欽定理，協(xié)方差函數(shù)的傅里葉變換是功率譜。該方法利用自協(xié)方差函數(shù)的傅里葉變換與和標(biāo)準(zhǔn)化的測井曲線同分布的隨機(jī)數(shù)進(jìn)行褶積，構(gòu)建精確的分形高斯模型。

首先，由序列R，即R(0),R(1),…,R(M/2-1),R(M/2),R(M/2-1),…,R(1)的離散傅里葉變換計(jì)算精確的功率譜S(k)：

(7)

對于所有k，所有的功率譜S(k)≥0。

然后，產(chǎn)生均值為μ的獨(dú)立同分布的高斯隨機(jī)數(shù)W(k),k=0,1,…,M-1，利用測井?dāng)?shù)據(jù)的方差σ2進(jìn)行標(biāo)準(zhǔn)化處理。

(8)

式中：“*”表示V(k)和V(M-k)是復(fù)共軛。

最后，用V(k)的離散傅里葉變換的前N個(gè)元素計(jì)算模擬序列Y(k)：

(9)

c=0,1,…,N-1

1.3.2 基于克里金插值的低頻先驗(yàn)信息

克里金插值是在考慮了信息樣品的幾何特征和變量的空間結(jié)構(gòu)信息后，為了達(dá)到線性無偏和最小估計(jì)方差的條件，對每個(gè)樣品值賦予一定權(quán)值，利用加權(quán)平均值法對待估計(jì)區(qū)域的未知量進(jìn)行估計(jì)的方法，也就是說是一種待定的滑動(dòng)加權(quán)平均方法。圖2為克里金估計(jì)示意圖，即根據(jù)周圍位置的數(shù)據(jù)來估計(jì)中間點(diǎn)位置的數(shù)值。以簡單克里金為例，表達(dá)式可簡要表示為

(10)

式中：C(ui,uj)表示采樣點(diǎn)ui,uj間的協(xié)方差;C(uj,u0)表示未采樣點(diǎn)u0與采樣點(diǎn)uj之間的協(xié)方差;λi為待求的克里金系數(shù)。然后將其代入(11)式中，即可得到簡單克里金的估計(jì)量和估計(jì)方差。

(11)

圖2 克里金估計(jì)示意圖解

1.3.3 混合初始模型建立

低頻模型是通過測井?dāng)?shù)據(jù)進(jìn)行克里金插值，然后進(jìn)行低通濾波得到的。這個(gè)初始模型常用于確定性反演中獲取地層趨勢。在獲得低頻模型m1和基于分形理論的高頻模型m2后，就可以進(jìn)行簡單的線性加權(quán)得到混合初始模型m：

(12)

式中：α為加權(quán)值，根據(jù)實(shí)際情況選取。如果原始測井曲線的儲層參數(shù)變化較劇烈，可適當(dāng)減小α;當(dāng)儲層參數(shù)變化較平緩時(shí)，可以適當(dāng)增大α。如圖3 和圖4，由于選取的數(shù)據(jù)波動(dòng)性不大，所以取α=0.7。

圖3為分形初始模型反演結(jié)果與真實(shí)測井曲線的對比。圖4為基于分形理論和克里金插值的混合初始模型反演結(jié)果與真實(shí)值的對比。由圖3可以看出，基于分形理論初始模型的反演結(jié)果波動(dòng)性較大，與真實(shí)值差異較大；而圖4紅色曲線所示的混合初始模型在一定程度上克服了抖動(dòng)，更好地給出了曲線的趨勢，與真實(shí)測井?dāng)?shù)據(jù)值更加接近。進(jìn)行反演時(shí)，采用混合初始模型進(jìn)行反演所需的計(jì)算時(shí)間更短，計(jì)算效率更高。

圖3 分形初始模型與真實(shí)值的對比

圖4 混合初始模型與真實(shí)值的對比

1.4 似然函數(shù)構(gòu)建

本文方法的似然函數(shù)可表示為

(13)

式中：s，μ，σ和H分別代表模擬地震波形的振幅、已知測井?dāng)?shù)據(jù)的物性參數(shù)的均值、標(biāo)準(zhǔn)差和Hurst系數(shù);sobs，μobs，σobs和Hobs分別是觀測地震數(shù)據(jù)、反演算法得到的物性參數(shù)的均值、標(biāo)準(zhǔn)差和Hurst系數(shù);σ是期望數(shù)據(jù)不確定性的標(biāo)準(zhǔn)差;R代表反演的儲層物性參數(shù);R0是確定性反演中的低頻約束信息;α，α1，α2和α3是可調(diào)參數(shù)，可通過不斷測試得到所需的參數(shù);s為由孔隙度、飽和度和泥質(zhì)含量計(jì)算得到的合成記錄，其中彈性參數(shù)和物性參數(shù)的關(guān)系由統(tǒng)計(jì)巖石物理模型確立。

(13)式似然函數(shù)中的第1項(xiàng)為反演算法中目標(biāo)函數(shù)的基本結(jié)構(gòu)(觀測合成記錄與反演合成記錄差的L2范數(shù))，第2項(xiàng)為確定性反演中常用的低頻約束信息，第3，4，5項(xiàng)為基于分形理論的約束信息(均值、標(biāo)準(zhǔn)差和Hurst系數(shù))。目標(biāo)函數(shù)的后4項(xiàng)依賴于井資料，所以在離井近的位置，相應(yīng)的權(quán)重可以取較大值;而遠(yuǎn)離井的位置，權(quán)重適當(dāng)?shù)販p小。

1.5 SA-PSO混合優(yōu)化算法

粒子群優(yōu)化算法(Particle Swarm Optimization,PSO)是一種基于群體的隨機(jī)優(yōu)化算法，最先由美國學(xué)者Kennedy等[16]通過模擬鳥群覓食過程中的遷徙和群聚行為提出。粒子群優(yōu)化算法通過個(gè)體之間的協(xié)作來尋找最優(yōu)解。PSO中粒子向量的速度和位置變化公式為

(14)

式中：V為粒子i的速度向量更新;x為粒子i的位置更新;ω是慣性加權(quán)值，表示舊速度對新速度的影響;R1和R2為兩個(gè)隨機(jī)向量，是在[0,1]上均勻分布的隨機(jī)數(shù);c1和c2是學(xué)習(xí)因子;xpbest是粒子的局部最優(yōu)位置;xgbest是粒子的全局最優(yōu)位置。

然而，當(dāng)粒子的飛行速度較大時(shí)，有利于全局搜索，但有可能飛過最優(yōu)解;當(dāng)粒子的飛行速度較小時(shí)，容易陷入局部最優(yōu)[17]。鑒于以上缺陷，高鷹等[18]建議應(yīng)用基于模擬退火的粒子群優(yōu)化算法。該算法不僅基本保持了粒子群優(yōu)化算法容易實(shí)現(xiàn)的特點(diǎn)，還改善了粒子群優(yōu)化算法擺脫局部極值點(diǎn)的能力，提高了算法的收斂速度和精度。

模擬退火算法(Simulated Annealing,SA)的思想最早由Metropolis等[19]提出，Rothman[20-21]首先將模擬退火算法應(yīng)用于地球物理問題?；谀M退火改進(jìn)的粒子群優(yōu)化算法的(SA-PSO)流程如圖5所示，其中SA選用快速模擬退火方法(Very Fast Simulated Annealing，VFSA)。

圖5 SA-PSO優(yōu)化算法流程

因此，首先根據(jù)分形理論和克里金插值得到具有高頻和低頻成分的初始模型，然后利用統(tǒng)計(jì)巖石物理模型建立彈性參數(shù)和物性參數(shù)之間的關(guān)系，最后依據(jù)SA-PSO優(yōu)化算法進(jìn)行優(yōu)化，就可以得到后驗(yàn)概率密度的樣本。

2 反演實(shí)例分析

2.1 一維數(shù)據(jù)測試

選取A井?dāng)?shù)據(jù)進(jìn)行反演算法測試與分析。首先求取與A井相鄰4口井的孔隙度、泥質(zhì)含量以及含水飽和度的均值、方差和Hurst系數(shù)，從而建立基于分形高頻模型和低頻信息的混合初始模型(圖4)，同時(shí)反演A井的孔隙度、泥質(zhì)含量和含水飽和度，然后與真實(shí)數(shù)據(jù)值進(jìn)行對比。圖6是實(shí)際數(shù)據(jù)與基于分形高頻初始模型和低頻先驗(yàn)信息的儲層物性參數(shù)反演方法得到的孔隙度、泥質(zhì)含量和含水飽和度的對比結(jié)果，可以看出反演結(jié)果和真實(shí)測井曲線匹配得很好。

圖6 反演結(jié)果與實(shí)際數(shù)據(jù)的對比

圖7和圖8分別是資料信噪比為4和1的反演結(jié)果，可以看出，只有部分層段與測井?dāng)?shù)據(jù)有一些差別，即使是在信噪比為1時(shí)，仍可很好地反演出儲層物性參數(shù)的趨勢，可以用于區(qū)分巖性和流體。通過加噪分析仍然可以得到可信的反演結(jié)果。

圖7 資料信噪比為4時(shí)反演結(jié)果與實(shí)際數(shù)據(jù)的對比

圖8 資料信噪比為1時(shí)反演結(jié)果與實(shí)際數(shù)據(jù)的對比

2.2 二維資料分析

二維地震資料來自于國內(nèi)某油田實(shí)際工區(qū)。首先從該實(shí)際資料中提取正相位子波，然后對其進(jìn)行基于分形高頻初始模型和低頻先驗(yàn)信息的儲層物性參數(shù)隨機(jī)反演。地震資料的縱向采樣率為2ms，時(shí)間范圍為1.10～1.30s，共有101道數(shù)據(jù)，圖9顯示的是疊后地震剖面。取該工區(qū)的5口井，基于除A井外的4口井進(jìn)行分段交會分析，建立統(tǒng)計(jì)巖石物理模型和先驗(yàn)信息，用A井(盲井)驗(yàn)證反演方法的可行性。

圖10，圖11和圖12分別是隨機(jī)反演得到的孔隙度、泥質(zhì)含量和含水飽和度，圖中標(biāo)示的曲線為A井的井曲線?？梢钥闯?，反演結(jié)果和A井的真實(shí)測井曲線匹配得很好，分辨率較高。對于孔隙度和泥質(zhì)含量來說，即使是3.5m左右的薄層都可以分辨得比較清晰。由于測井解釋的含水飽和度曲線分辨率稍差，所以反演結(jié)果的分辨率不是很高，但還是與解釋的測井曲線吻合較好。該二維實(shí)際資料的應(yīng)用結(jié)果驗(yàn)證了本文反演方法的可行性。

圖9 疊后地震剖面

圖10 反演的孔隙度剖面

圖11 反演的泥質(zhì)含量剖面

圖12 反演的含水飽和度剖面

3 結(jié)論與認(rèn)識

我們引入統(tǒng)計(jì)巖石物理模型反演儲層物性參數(shù)，彌補(bǔ)了確定性巖石物理模型的局限性。該反演方法直接采用精確的Zoeppritz方程計(jì)算正演合成地震記錄，可以減少誤差。基于分形方法產(chǎn)生的初始模型具有高頻信息，同時(shí)可以減少假頻的出現(xiàn)?？死锝鸩逯到Y(jié)果具有低頻趨勢，所以二者的混合初始模型同時(shí)具有高頻和低頻信息。利用基于模擬退火改進(jìn)的粒子群優(yōu)化算法進(jìn)行優(yōu)化，避免了粒子群優(yōu)化算法易于陷入局部極值的缺點(diǎn)，提高了算法的精度和收斂速度?；诜中胃哳l初始模型和低頻先驗(yàn)信息的物性參數(shù)反演通過貝葉斯理論框架結(jié)合了分形高斯模型算法、克里金插值、SA-PSO優(yōu)化算法和統(tǒng)計(jì)巖石物理理論。即使在信噪比相對較低時(shí)，仍然可以反演出合理的儲層物性參數(shù)信息，從而為孔隙度、泥質(zhì)含量以及含水飽和度的預(yù)測提供了一種可行的方法。

參 考 文 獻(xiàn)

[1] Srivastava R P,Sen M K.Stochastic inversion of prestack seismic data using fractal-based initial models[J].Geophysics,2010,75(3): R47-R59

[2] Tao Y,Spikes K,Sen M K.Stochastic seismic inversion using both fractal and low-frequency priors[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011,2732-2735

[3] Mukerji T,J?rstad A,Avseth P,et al.Mapping lithofacies and pore-fluid probabilities in a North Sea reservoir:seismic inversions and statistical rock physics[J].Geophysics,2001,66(4):988-1001

[4] Mukerji T,Avseth P,Mavko G,et al.Statistical rock physics:combining rock physics,information theory,and geostatistics to reduce uncertainty in seismic reservoir characterization[J].The Leading Edge,2001,20(3):313-319

[5] Eidsvik J,Avseth P,Omre H,et al.Stochastic reservoir characterization using prestack seismic data[J].Geophysics,2004,69(4):978-993

[6] Spikes K,Mukerji T,Dvorkin J,et al.Probabilistic seismic inversion based on rock-physics models[J].Geophysics,2007,72(5):R87-R97

[7] Grana D,Della Rossa E.Probabilistic petrophysical-properties estimation integrating statistical rock physics with seismic inversion[J].Geophysics,2010,75(3):O21-O37

[8] Grana D,Dvorkin J.The link between seismic inversion,rock physics,and geostatistical simulations in seismic reservoir characterization studies[J].The Leading Edge,2011,30(1):54-61

[9] 胡華鋒,印興耀,吳國忱.基于貝葉斯分類的儲層物性參數(shù)聯(lián)合反演方法[J].石油物探,2012,51(3):225-232

Hu H F,Yin X Y,Wu G C.Joint inversion of petrophysical parameters based on Bayesian classification[J].Geophysical Prospecting for Petroleum,2012,51(3):225-232

[10] 夏麗娜.基于貝葉斯理論框架的儲層參數(shù)地震反演方法研究[D].青島:中國石油大學(xué)(華東),2013

Xia L N.The study of petro-physical properties seismic inversion based on Bayesian theory[D].Qingdao:China University of Petroleum (Huadong),2013

[11] Kj?nsberg H,Hauge R,Kolbj?rnsen O,et al.Bayesian Monte Carlo method for seismic predrill prospect assessment[J].Geophysics,2010,75(2):O9-O19

[12] 李宏兵,張佳佳,姚逢昌.巖石的等效孔隙縱橫比反演及其應(yīng)用[J].地球物理學(xué)報(bào),2013,56(2):608-615

Li H B,Zhang J J,Yao F C.Inversion of effective pore aspect ratios for porous and its applications[J].Chinese Journal of Geophysics,2013,56(2):608-615

[13] Scales J A,Smith M L,Treitel S.Introductory geophysical inverse theory[M].Colorado:Samizdat Press,2001:1-193

[14] Hewett T A.Fractal distributions of reservoir heterogeneity and their influence on fluid transport[J].Expanded Abstracts of 61stSPE Annual Conference,1986,15386

[15] 韓忠玲.分形高斯噪聲Hurst參數(shù)估計(jì)的實(shí)驗(yàn)評價(jià)[D].上海:華東師范大學(xué),2008

Han Z L.Experimental evaluation on methods of the Hurst parameter estimation using fractional Gausssian noise[D].Shanghai:East China Normal University,2008

[16] Kennedy J,Eberhart R.Particle swarm optimization[C]∥Proceedings of IEEE international conference on neural networks.Neural Networks:Institute of Electrical and Electronics Engineers,1995:1942-1948

[17] 朱麗莉,楊志鵬,袁華.粒子群優(yōu)化算法分析及研究進(jìn)展[J].計(jì)算機(jī)工程與應(yīng)用,2007,43(5):24-27.

Zhu L L,Yang Z P,Yuan H.Analysis and development of particle swarm optimization[J].Computer Engineering and Applications,2007,43(5):24-27

[18] 高鷹,謝勝利.基于模擬退火的粒子群優(yōu)化算法[J].計(jì)算機(jī)工程與應(yīng)用,2004,40(1):47-50

Gau Y,Xie S L.Particle swarm optimization algorithms based on simulated annealing[J].Computer Engineering and Applications,2004,40(1):47-50

[19] Metropolis N,Rosenbluth A,Rosenbluth M,et a1.Equation of state calculations by fast computing machines[J].The Journal of Chemical Physics,1953,21(6):1087-1092

[20] Rothman D H.Nonlinear inversion statistical mechanics and residual statics estimation[J].Geophysics,1985,50(12):2784-2796

[21] Rothman D H.Automatic estimation of large residual statics correction[J].Geophysics,1986,51(2):337-346