基于ASSA方法的等值性電測深數據反演及數據誤差相關性分析

付代光 廖錦芳 周黎明 肖國強 王法剛 羅 榮

(①長江科學院水利部巖土力學與工程重點實驗室，湖北武漢 430010; ②湖北省頁巖氣開發有限公司，湖北武漢 430071)

·非地震·

基于ASSA方法的等值性電測深數據反演及數據誤差相關性分析

付代光*①廖錦芳②周黎明①肖國強①王法剛①羅 榮①

(①長江科學院水利部巖土力學與工程重點實驗室，湖北武漢 430010; ②湖北省頁巖氣開發有限公司，湖北武漢 430071)

付代光，廖錦芳，周黎明，肖國強，王法剛，羅榮.基于ASSA方法的等值性電測深數據反演及數據誤差相關性分析.石油地球物理勘探，2017,52(5):1077-1084.

針對直流電測深數據等值性反演的病態性(等值性層參數及層參數間的狹長平坦的強相關性)，嘗試采用自適應單純形模擬退火(ASSA)方法反演具等值性電測深模型，通過與Lamarckian-marked-constraint(LMC)算法對比，發現ASSA方法能較好地解決電測深等值性反演問題，且反演精度較高。目前地球物理反演多基于數據誤差為高斯分布和空間不相關(非對角元素為零)的假設，對貝葉斯反演而言，后一假設如果運用不當往往會低估反演結果的不確定性，從而導致對反演結果可靠性的誤判。為了合理評價反演結果的不確定性，文中采用Runs-test非參數化判別方法，通過對等值性模型的歸一化數據殘差的判別，確定數據誤差是空間相關的。由此得出的協方差矩陣應用于貝葉斯反演，進而對反演結果的不確定性、相關性等進行更合理、準確的評價。

電測深等值性 貝葉斯反演 ASSA 數據誤差相關性 不確定性 相關性

1 引言

直流電測深方法是直流電阻率法的一種，通過控制極距達到探測不同深度目標體的目的。該方法應用范圍廣，不但可以應用于尋找金屬和非金屬礦產、圈定油藏邊界、探測天然氣水合物，而且在水文與工程地質、城市環境調查、建筑基礎勘查等方面也發揮著重要的作用[1-5]。

直流電測深具有等值性，而具有等值性的電測深數據反演是病態的[6]。電測深等值性是由其方法本身的特性所決定的，無法從理論上消除，因而常采用改進的反演算法提高反演精度，但相關的研究并不多見。Basokur等[7]提出了混合算法Lamarckian-marked-constraint(LMC)，該算法對傳統的遺傳算法進行了改進，將自然選擇方式變為人為選擇，即增加約束條件，對不滿足條件的模型進行淘汰，同時將最小二乘與標記約束遺傳算法相結合，利用最小二乘的梯度信息，提高算法對強相關性參數反演的能力。但該算法只能應用于有限個參數的反演中，另外對不符合實際情況的約束條件，反演算法將會失效；同時，該算法的反演過程消耗大量的計算時間，尤其二維反演更是如此，因為二維正演過程需要大量的計算時間，網格越小所需時間就越多，因而限制了在實際中的應用。為了實現基于LMC方法的二維電阻率成像， Ak?a等[8]采用Delaunay 三角剖分法極大地減少了地電模型的網格數，從而基于反演結果實現了巖性單元的劃分。Fallat等[9]和Dosso等[10]提出的單純形模擬退火(SSA)和自適應單純形模擬退火是一種局部與全局、線性與非線性相結合的確定性搜索算法，該方法對解決強參數相關性問題非常有效，而且無需求導計算。針對ASSA方法的特點，本文嘗試采用該方法對等值性電測深數據進行反演。

貝葉斯反演是一種把觀測數據和模型參數都看成隨機變量的反演方法。傳統的反演過程是一個尋找最優解的過程，它只提供一個確定的反演結果，但并未定量描述反演結果的可靠度，也未將含有誤差和噪聲的數據加入到反演過程中；而貝葉斯反演能對上述問題提供較滿意的解決方法[11-13]。貝葉斯反演是關于協方差矩陣最大似然函數的求解過程，如果假設協方差矩陣中非對角元素為零，那么將得到貝葉斯反演中常規的目標函數公式(非對角元素為零)。目前很多地球物理反演都是基于此假設的[14-18]，但該假設往往造成對反演結果不確定性的低估。為了對反演結果的不確定性進行合理的評價，Dosso等[19，20]提出了一種非參數化方法Runs-test(亦可稱為Runs檢驗)用于檢驗協方差矩陣元素是否空間相關。通過Runs檢驗得到的協方差矩陣應用于貝葉斯反演中，能對反演結果的不確定性、相關性等給出合理的評價。

本文首先從正演角度對等值性模型的一維和二維參數變化引起的誤差變化程度曲線的特征進行了分析，并應用ASSA算法對等值性電測深模型進行反演，通過與LMC方法反演結果的對比驗證了ASSA算法的有效性。為了分析數據誤差相關性對反演結果的不確定性的影響，本文對比了數據誤差空間相關與不相關反演結果的差異，并采用非參數化方法Runs檢驗數據誤差相關與否。將Runs檢驗判定出的數據誤差相關性結果應用到貝葉斯反演中，從而對反演結果的不確定性、參數間相關性等進行更加合理的評價。

2 貝葉斯反演

2.1 貝葉斯定理

假定d表示實測數據，m為模型參數，則對于任意的d和m，滿足如下的貝葉斯關系

(1)

式中：p(m|d)表示已知d時，所確定的m的函數，稱為后驗概率密度函數(probability density function，簡稱為后驗PDF)；同理，p(d|m)表示在給定m時，所得d的概率密度函數；p(m)稱為m的先驗PDF；p(d)是d的先驗概率密度，通常d是已知的，因此該項往往都看作常數，此時式(1)可寫為

p(m|d)∝p(d|m)p(m)

(2)

p(d|m)與似然函數L(m)是對等的，即L(m)∝p(d|m)，因此式(2)可進一步寫為

p(m|d)∝L(m)p(m)

(3)

式中L(m)∝exp[-E(m)]，E(m)為目標函數。廣義目標函數包含了數據誤差p(m)和先驗信息，即

φ(m)=E(m)-lnp(m)

(4)

聯立式(1)～式(4)，并將向量m的PDF歸一化得到

(5)

式中V表示模型的積分域。式(5)即所求問題的解。為了對反演結果進行更合理的解釋，需要對反演參數進行估計、分析不確定性并計算參數相關性，即對后驗概率進行積分求解，具體求解公式參見文獻[21]。

2.2 自適應單純形模擬退火

ASSA方法是將模擬退火與單純形方法相結合的一種非線性方法，該方法具備全局與局部搜索能力，與傳統的全局優化方法(如模擬退火和遺傳算法等)和局部優化方法(如擬牛頓法、最小二乘法等)相比具有以下特點：

(1)反演結果與初始值無關，克服了局部算法對初始估計值的較高要求；

(2)能夠保留搜索過程中得到的最優解，克服了模擬退火不能保存最優解的缺點；

(3)由于單純形具有梯度搜索能力(但不需要計算導數或靈敏度矩陣)，因此該方法反演效率高，尤其對相關性強或在收斂值附近的問題能夠高效解決。

自適應單純形模擬退火的自適應體現在對參數擾動的自動更新，其模型擾動公式為

(6)

ASSA方法將單純形嵌入到模擬退火中，因此每一次的降溫過程不再是對單一模型進行操作，而是對所有模型的單純形進行操作。通過模擬退火擾動產生的初始模型，經由單純形算法進行優化，再由模擬退火對優化參數進行隨機擾動，然后利用Metropolis準則判斷新生成的模型是否被接受。在上述過程中需要注意的是，單純形會產生M+1個模型(M表示模型m的元素個數)，且模型每次都會被保留，因此避免了模擬退火中好的模型被丟棄；同時，在模擬退火后期，隨著溫度降低(溫度為0°C時)，會變成純粹的單純形的搜索模式[8]。

3 協方差計算

對貝葉斯反演中數據協方差矩陣的計算采用Dosso[19]提出的計算方式。數據協方差矩陣可定義為

Cd=〈(d-〈d〉)T(d-〈d〉)〉

(7)

式中〈·〉表示總體平均。由于采集數據的誤差不能代表整個樣本的誤差，因此該公式很難應用于實際。

如果將協方差矩陣近似為最常用的形式——對角協方差矩陣，即非對角元素為零的矩陣Cd=σ2I(I是單位矩陣)，那么此時的非線性貝葉斯公式的似然函數(目標函數)可簡化為常用的二范數形式。對角協方差矩陣其實忽略了非對角元素的相關性，因而容易導致對參數不確定性的低估。為了合理估計參數的不確定性，本文采用優化模型的數據殘差計算全元素協方差矩陣

(8)

i=1,…,N;j=1,…，N

(9)

(10)

4 ASSA反演

為驗證ASSA算法的有效性，本文選擇具有等值性特征的四層地電模型[7](表1)。表中第一行給出了模型各層的電阻率和厚度。在驗證ASSA方法之前，先對該模型進行正演分析。

表1 模型參數值及LMC與ASSA反演結果

4.1 電測深正演

參數敏感性與參數相關性是影響反演結果的關鍵因素。單個參數變化引起的適值變化程度，稱為敏感性[10]。如果隨著參數變化，適值發生較明顯改變(曲線的曲率較大)，則說明參數敏感，較窄的參數敏感性有利于反演；反之，則不利于參數的反演。圖1給出了第二、第三層電阻率和厚度的參數適值范圍。從圖中可以看出各參數相對敏感，在不考慮參數相關性、等值性等其他因素的情況下，各參數的反演并不會太復雜。

從圖2部分參數間的誤差能量等值線圖可知，參數間存在強相關性，且這種強相關性表現為狹長型。其中，ρ2-h2(圖2g)、ρ3-h3(圖2h)表現得尤為明顯，因此對其反演時更難捕捉到有效的參數值。因此，單純采用線性反演方法需要有較好的初始猜測值，否則對這種狹長相關性的參數反演可能會失敗；而對非線性全局優化反演方法而言，則要求大量空間搜索(大量的計算時間)，或許可以得到較為合理的解，有時也可能得到局部最優解(對等值性電測深而言，可能為等值性參數對)。因此該模型的反演具有一定的挑戰性。

圖1 一維參數誤差能量斷面圖

圖2 誤差能量等值線圖(黑色區域為低適值區)

4.2 電測深反演

為驗證ASSA算法的有效性，本文選擇Basokur等[7]提出的LMC反演方法進行對比。LMC法是一種將線性最小二乘法與遺傳算法相結合的混合算法。為了方便對比，設定與文獻[7]相同的參數搜索范圍(見表1)。ASSA反演參數為：初始溫度為0.03°C，迭代1500次，每個溫度擾動5次。圖3為ASSA反演的適值變化與迭代次數的關系，從圖中可以看出適值降低速度較快，即算法的收斂速度較快。根據圖4所示的反演結果可知：視電阻率—極距和電阻率—深度曲線的反演結果精度較高，具體的反演值見表1。由表1可見，ASSA的反演精度與LMC方法相當甚至更高。但需要注意的是，單純形采用近似梯度搜索方法，而對于等值性強且窄的相關性參數而言，有時單次反演也未必一定能取得很好的結果，一般取幾次反演結果的均值為宜。

圖3 ASSA反演的適值與迭代次數的關系曲線

5 數據誤差分析

數據誤差的相關與否，決定了目標函數的形式，而目標函數的形式是否會對反演結果產生影響，為此，采用文獻[7]的模型，ASSA反演參數設定為：迭代次數為500次，初始溫度為0.3°C，每個溫度擾動20次，均取2次反演結果的平均值。圖5為電阻率反演結果。從圖5b看出，空間相關協方差矩陣對反演結果的擬合精度影響較大，其中對第三層電阻率和厚度的影響最為明顯； 但從圖5a發現，視電阻率—極距的反演擬合精度較好，究其原因，可能是由數據誤差的空間相關性及模型的等值性引起電阻率—深度剖面的反演擬合效果較差。在對電阻率剖面進行反演時，如何判斷數據誤差是否空間相關，至關重要。

本文采用非參數化Runs方法判斷數據誤差是否相關。從圖6a可以看出，對角協方差矩陣的自相關圖形主峰較寬，而且跳躍幅度較大，因此初步判定數據誤差相關，通過Runs定量計算可知p=0.0001，再次證明數據誤差相關。由圖6b也可看出，自相關的波峰較窄，且p=0.94，說明數據誤差是空間相關的。需要注意的是，對于不同輪次的反演結果，數據殘差自相關圖形和p值大小會有所變化，但這并不影響最終的結論，因此不同輪次的反演對殘差自相關結果的估計不敏感。

前文已證明數據誤差空間相關，因此在貝葉斯反演和后驗概率積分時都應將其考慮進去。圖7～圖9為數據誤差相關時所得的一維邊緣概率分布。由圖7可見，ρ3和h3的不確定性較大，呈雙峰狀，且峰值均不在理論值附近，ρ3和h3在其上、下邊界附近的采樣值較多；h1和ρ4主峰較窄，h1基本在9～10間變化，ρ4基本在30～50間變化，且與理論值基本一致，不確定性小；ρ2和h2的峰值較緩，ρ2在20～80間變化，h2在5～20間變化，置信區間相對較小，不確定性相對較大。根據圖8可知，ρ3-ρ2、ρ3-h2、h3-ρ2、h3-h2的二維邊緣概率分布的采樣點極其分散，這是由于ρ3和h3為多峰多模態，ρ2和h2的不確定性較大，置信區間較小，且ρ2、h2、ρ3和h3彼此間存在不同程度的相關性(圖9)，影響了各參數的采樣值；對于h1和ρ4這樣主峰較窄的參數，ρ3-h1、ρ3-ρ4、h3-h1、h3-ρ4的二維邊緣概率分布的采樣點相對集中。結合圖8和圖9綜合分析發現，ρ2-h2呈強線性正相關，而ρ3-h3呈強非線性負相關，此結果與圖2中的結果也是吻合的。這些相關性將會影響反演結果的精度，同時也會影響參數的不確定性。強的相關性會導致反演精度降低，增加反演結果的不確定性，這也是ρ3和h3在反演過程中很難得到較好結果的原因。如何有效降低參數間的相關性、提高反演精度，是亟待解決的問題。

圖4 采用文獻[7]參數的ASSA反演結果

圖5 電阻率反演結果

圖6 Runs-test法反演結果

圖7 部分參數一維邊緣概率分布

圖8 部分參數二維邊緣概率分布白色虛線表示參數真值

圖9 模型參數間相關系數

6 結論

(1)通過對等值性電測深數據一維和二維正演分析發現，單個參數敏感性較窄，利于參數反演。但從二維正演曲線發現，大多數參數間存在強相關性，且在真值附近存在許多等值性參數對，因此要獲得好的反演結果較困難。本文根據ASSA方法的近似梯度搜索能力和收斂速度快等特點，對等值性電測深數據進行了反演，其反演結果精度較高，驗證了ASSA算法的優越性。

(2)通過對比數據誤差空間相關與不相關的反演結果發現，數據誤差相關與否對反演精度具有較大的影響。為了對數據誤差的相關性進行判斷，本文采用Runs檢驗法定性和定量地判定了電測深模型數據殘差的相關性，最終證明其數據誤差是空間相關的。

(3)通過對后驗概率數值積分，獲得反演結果的一維和二維邊緣概率分布、相關系數等特征量，由此得到本文所述模型的第二層與第三層參數的不確定性大、且第三層的厚度和電阻率呈多模態性(雙峰狀)的結論。參數間存在不同程度的相關性，其中以第二和第三層相關性最為明顯，第三層厚度和電阻率呈強非線性負相關，第二層參數間呈強正相關。相關性影響了反演結果精度并增大了參數的不確定性。

[1] 王橋，萬漢平，王聞文等.綜合物探方法在鋁土礦勘查中的應用.地球物理學進展,2012,27(2):709-714. Wang Qiao，Wan Hanping，Wang Wenwen et al.The application of integrated geophysical exploration in Bauxite.Progress in Geophysics,2012,27(2):709-714.

[2] 翁愛華，祝嗣安.天然氣水合物的近海底直流電測深響應.石油地球物理勘探,2010,45(3)：458-461. Weng Aihua,Zhu Si’an.Near sea-bottom DC sounding response for gas hydrate.OGP,2010,45(3)：458-461.

[3] 張天倫.用直流電阻率法確定油氣藏邊界的初步實驗.石油地球物理勘探,1990，25(5)：584-593. Zhang Tianlun.An experiment in locating reservoir boundary using direct-current resistivity method.OGP，1990，25(5)：584-593.

[4] 王玉，郭秀軍，賈永剛等.土壤—地下水系統非水相液體污染區的電性特征及電阻率法探測.地球物理學進展,2009,24(6):2316-2323. Wang Yu，Guo Xiujun，Jia Yonggang et al.Overview of the electrical characters and the electrical survey methods for the soil-groundwater system in organic contaminated areas.Progress in Geophysics,2009,24(6):2316-2323.

[5] 賈志寬，安西峰，李振峰等.非常規電測深法在工程地質勘查中的應用.地球物理學進展,2006,21(1):296-299. Jia Zhikuan,An Xifeng,Li Zhenfeng et al.Application of unconventional electric detecting method in engineering geology survey.Progress in Geophysics,2006,21(1):296-299.

[6] Alvarez J P F,Martinez J L F,Perez C O M.Feasibility analysis of the use of binary genetic algorithms as importance samplers application to a 1-D DC resistivity inverse problem.International Association for Mathematical Geology,2008,40(4):375-408.

[7] Basokur A T,Akca I,Siyam N W A.Hybrid genetic algorithms in view of the evolution theories with application for the electrical sounding method.Geophy-sical Prospecting,2007，55(3)：393-406.

[8] Ak?a I,Basokur A T.Extraction of structure-based geoelectric models by hybrid genetic algorithms.Geophysics,2010，75(1):F15-F22.

[9] Fallat M R,Dosso S E.Geoacoustic inversion via local,global,and hybrid algorithms.Acoustical Society of America.1999，105 (6):3219- 3230.

[10] Dosso S E,Wilmut M J,Lapinski A L S.An adaptive-hybrid algorithm for geo-acoustic inversion.Acous-tical Society of America，2001，26(3):324-336.

[11] Guo R W,Dosso S E,Liu J X.Non-linearity in Baye-sian 1-D magnetotelluric inversion.Geophysical Journal International，2011，185(2):663-675.

[12] Dong H F,Dosso S E.Bayesian inversion of interface-wave dispersion for seabed shear-wave speed profiles.IEEE Journal of Oceanic Engineering,2011,36(1):1-11.

[13] 尹彬,胡祥云.基于并行回火技術優化的大地電磁數據貝葉斯變維反演.石油地球物理勘探,2016,51(6):1212-1218. Yin Bin,Hu Xiangyun.Magnetotelluric Bayesian trans-dimensional inversion improved by parallel tempering algorithm.OGP,2016,51(6):1212-1218.

[14] 宗兆云,印興耀,張繁昌.基于彈性阻抗貝葉斯反演的拉梅參數提取方法研究.石油地球物理勘探,2011,46(4):598-604,609. Zong Zhaoyun,Yin Xingyao and Zhang Fanchang.Elastic impedance Bayesian inversion for Lame parameters extracting.OGP,2011,46(4):598-604,609.

[15] 張豐麒，金之鈞，盛秀杰等．一種穩健的彈性阻抗反演方法．石油地球物理勘探，2017，52(2)：294-303. Zhang Fengqi,Jin Zhijun,Sheng Xiujie et al.A stable elastic impedance inversion approach.OGP,2017，52(2)：294-303

[16] 張廣智，潘新朋，孫昌路等．縱橫波聯合疊前自適應MCMC反演方法．石油地球物理勘探，2016，51(5)：938-946. Zhang Guangzhi,Pan Xinpeng,Sun Changlu et al.PP- & PS-wave prestack nonlinear inversion based on adaptive MCMC algorithm.OGP,2016，51(5)：938-946.

[17] 陳建江,印興耀.基于貝葉斯理論的AVO三參數波形反演.地球物理學報,2007,50(4):1251-1260. Chen Jianjiang,Yin Xingyao.Three parameter AVO waveform inversion based on Bayesian theorem.Chinese Journal of Geophysics,2007,50(4):1251-1260

[18] 侯棟甲，劉洋，胡國慶等．基于貝葉斯理論的疊前多波聯合反演彈性模量方法．地球物理學報，2014,57(4)：1251-1264. Hou Dongjia,Liu Yang,Hu Guoqing et al.Prestack multiwave joint inversion for elastic moduli based on Bayesian theory.Chinese Journal of Geophysics,2014,57(4):1251-1264.

[19] Dosso S E,Nielsen P L,Wilmut M J.Data error covariance in matched-field geoacoustic inversion.A-coustical Society of America，2006，119(1):208-219.

[20] Dosso S E,Holland C W.Geoacoustic uncertainties from viscoelastic inversion of seabed reflection data.IEEE Journal of Oceanic Engineering,2006,31(3):657-671.

[21] 付代光，劉江平，周黎明等.基于貝葉斯理論的軟夾層多模式瑞雷波頻散曲線反演研究.巖土工程學報，2015,37(2):321-329. Fu Daiguang,Liu Jiangping,Zhou Liming et al.Inversion of multimode Rayleigh-wave dispersion curves of soft interlayer based on Bayesian theory.Chinese Journal of Geotechnical Engineering,2015,37(2):321-329.

(本文編輯：劉海櫻)

付代光 1987年生;2011年畢業于中國地質大學(武漢)應用地球物理專業，獲學士學位；2013年畢業于中國地質大學(武漢)地質工程專業，獲碩士學位；目前在長江水利委員會長江科學院主要從事地球物理正、反演的相關研究。

1000-7210(2017)05-1077-08

P631

A

10.13810/j.cnki.issn.1000-7210.2017.05.022

*湖北省武漢市江岸區后九萬方長江科學院巖土重點實驗室，430010。Email：fudaiguang@163.com

本文于2016年12月1日收到，最終修改稿于2017年7月27日收到。

本項研究受國家自然科學基金青年基金資助項目(41202223、51409013)資助。