楊 雯， 李予國,2??， 段雙敏， 韓 波， 羅 鳴,2

(1. 中國海洋大學海洋地球科學學院， 山東 青島 266100； 2. 青島海洋科學與技術試點國家實驗室 海洋礦產(chǎn)資源評價與探測技術功能實驗室， 山東 青島 266237； 3. 中國地質(zhì)大學(武漢) 地質(zhì)調(diào)查研究院，湖北 武漢 430074)

大地電磁法(Magnetotelluric,MT)是利用天然交變電磁信號研究地球內(nèi)部電性結(jié)構(gòu)的一種地球物理方法，現(xiàn)已廣泛應用于油氣資源勘探和俯沖帶、地殼與地幔電性結(jié)構(gòu)研究中[1-3]。大地電磁反演是獲得地下介質(zhì)電阻率信息的重要方法。然而，MT反演存在多解性，如何減少MT反演的多解性一直是重要的研究方向。自從Karl Pearson提出卡方檢驗以來，在地球物理反演中便一直使用均方根擬合差(RMS misfit)衡量反演模型響應與實測數(shù)據(jù)的擬合程度[4]。

1980年代末以來，為了改善大地電磁反演問題的不適定性，吉洪諾夫(Tikhonov)正則化方法開始被應用于MT反演中，并取得了迅速發(fā)展[5-7]。目前，常用的大地電磁反演方法包括OCCAM反演方法、非線性共軛梯度法(NLCG)、擬牛頓法和高斯-牛頓法等[7-10]。這些反演方法關于反演模型響應與觀測數(shù)據(jù)擬合程度的描述大都依賴于目標函數(shù)中的均方根擬合差，卻沒有考慮反演模型響應與觀測數(shù)據(jù)之差(常稱為殘差)的統(tǒng)計特征，這可能導致反演模型響應形態(tài)與觀測數(shù)據(jù)形態(tài)出現(xiàn)偏離較大的情況，尤其是當觀測數(shù)據(jù)存在較大誤差時，甚至可能得到錯誤的反演模型[11]。于是用于衡量實測數(shù)據(jù)與反演模型響應之間關系的其它統(tǒng)計量開始被應用于大地電磁反演中。Smith 和 Booker將斯皮爾曼(Spearman)相關系數(shù)應用于大地電磁(MT)反演結(jié)果評價中[12]。Jones利用DW統(tǒng)計量衡量觀測數(shù)據(jù)與模型響應之間殘差的自相關性[11]，并將其應用于加拿大Great Slave Lake地區(qū)大地電磁測深數(shù)據(jù)反演中[13]。將統(tǒng)計學度量應用于MT反演中，一方面可以評價反演結(jié)果的合理性，另一方面可以為反演算法的改進提供依據(jù)，以便進一步優(yōu)化反演方法。但是目前關于這方面的研究還不多見。本文將DW統(tǒng)計量應用于高斯-牛頓大地電磁一維反演中，通過修改目標函數(shù)，實現(xiàn)了基于DW統(tǒng)計量的大地電磁一維反演方法(以下簡稱為DW反演方法)。通過模型算例分析了該方法的有效性，并與傳統(tǒng)的高斯-牛頓反演結(jié)果進行了對比。

1 反演方法

1.1 DW統(tǒng)計量

杜賓-沃森(Durbin-Watson，簡稱DW)統(tǒng)計方法是用于檢驗回歸分析中殘差一階自相關性的一種方法，

該方法常用于檢驗計量經(jīng)濟學中的自相關問題[14-15]。一階自相關DW統(tǒng)計量的定義式為[11,16-17]：

(1)

圖1 序列自相關性與DW統(tǒng)計量的關系

表1 部分DW統(tǒng)計量臨界值分布表(k=2)

通常地，DW統(tǒng)計量與自相關系數(shù)R具有如下關系[16,18]：

DW≈2(1-R)。

(2)

式中，R為一階自相關系數(shù)，其表達式為：

(3)

一階自相關系數(shù)R的取值范圍為(-1，1)。通常認為，當|R|<0.2時，殘差序列的自相關性極弱。

1.2 DW反演方法

大地電磁反演問題實質(zhì)上是目標函數(shù)的最優(yōu)化問題。基于吉洪諾夫正則化反演思想，反演目標函數(shù)中通常包含數(shù)據(jù)擬合差項和正則化項[19]：

φ(m)=φd(m)+λ1φm(m)。

(4)

式中：m=[m1,m2,···,mM]T為模型參數(shù)向量，由反演模型中各地層的電阻率和厚度構(gòu)成；φd(m)為觀測數(shù)據(jù)與理論模型響應的擬合差，可表示為：

(5)

(6)

式中：Wm為拉普拉斯算子的差分近似；mref為參考模型。

在傳統(tǒng)的反演方法中，反演目標函數(shù)僅考慮了反演模型響應與觀測數(shù)據(jù)的擬合差，但沒有考慮其殘差的統(tǒng)計學特征。本文將DW統(tǒng)計量加入反演目標函數(shù)中，通過反演傾向于尋找在擬合差和殘差自相關性兩個方面都可以接受的模型。經(jīng)過修改后的反演目標函數(shù)具有如下形式[11]：

φ(m)=φd(m)+λ1φm(m)+λ2(DW-2)2。

(7)

上式中等號右端第三項是將DW統(tǒng)計量歸一化為零。在反演中，正則化因子λ1從一個較大的初值開始，隨迭代次數(shù)的增加逐漸減?。欢c此相反，對于權重系數(shù)λ2，則先給定一個較小的初值，在迭代過程中緩慢遞增。λ1遞減和λ2遞增的系數(shù)通常位于(1，2)之間。采用這種策略的基本思想是，在迭代反演的前期，注重數(shù)據(jù)擬合和模型約束，而在反演的后期通過增加DW項的權重，減少實測數(shù)據(jù)與反演模型響應之間殘差的自相關性。

本文采用高斯-牛頓最優(yōu)化方法求解目標函數(shù)式(7)的極小值。迭代反演過程中，模型更新量為：

(8)

其中gk和Hk分別為第k次迭代時目標函數(shù)φ(m)的梯度向量和Hessian矩陣，其表達式為：

(9)

(10)

在傳統(tǒng)的高斯-牛頓反演方法中，gk和Hk分別僅由式(9)和式(10)右端的前兩項構(gòu)成，而除該兩項以外的其它項反映了DW約束項對梯度向量和Hessian矩陣的貢獻。其中：

(11)

Hdw=2QJTKTKJ-2PQ2LLT+8Q2JTKTKe(LE)T+
8PQ3LE(LE)T。

(12)

1.3 合成數(shù)據(jù)添加誤差

大地電磁阻抗Z可以寫成下列指數(shù)形式

Z=reiφ。

(13)

由式(13)可以得到阻抗微小變化dZ的表達式

dZ=dreiφ+reiφ·idφ。

(14)

式(14)除以式(13)，得

(15)

上式表明，大地電磁阻抗Z的相對變化相當于其幅值的相對變化和相位的絕對變化。大多數(shù)時間序列分析表明，大地電磁阻抗的實部和虛部通常具有相等的誤差。于是，有

(16)

這意味著，對于大地電磁阻抗Z，其相位的絕對誤差與幅值的相對誤差相對應。也就是說，如果給大地電磁阻抗幅值添加相對誤差er, 則對應的應該給阻抗相位添加er弧度或er·180/π角度的絕對誤差。

(17)

大地電磁理論模型合成數(shù)據(jù)反演中，通常會給正演合成數(shù)據(jù)加入一定的隨機噪聲，并將其結(jié)果作為反演數(shù)據(jù)。然而，若直接給視電阻率添加相對誤差，可能會使得所添加的誤差具有自相關性，以2.1節(jié)將要討論的三層高阻薄層模型為例，對此進行說明。

首先，在頻率范圍10-3～103Hz之間選取21個對數(shù)等間隔的頻點，正演計算得到三層高阻薄層模型21個頻點的視電阻率和阻抗相位。然后，在視電阻率數(shù)據(jù)中加入1%的隨機相對誤差噪聲，相應地在阻抗相位數(shù)據(jù)中加入0.286 5°的隨機絕對誤差。用這樣的方法得到的視電阻率誤差和相位誤差隨周期的變化關系如圖2所示。

圖2 視電阻率誤差(a)和阻抗相位誤差(b)分布特征

從圖2中可以看到，阻抗相位誤差在周期范圍(10-3～103s)內(nèi)是隨機分布的。阻抗相位誤差的DW統(tǒng)計量DWpha=1.74和自相關系數(shù)Rpha=0.036，這表明阻抗相位誤差不具有自相關性。長周期視電阻率的誤差是在零線附近上下隨機分布的，而短周期視電阻率的誤差并非是隨機分布的。視電阻率誤差的DW統(tǒng)計量DWrho為 0.81，自相關系數(shù)Rrho為0.568，這表明視電阻率誤差具有一定的自相關性。圖3分別繪出了視電阻率誤差序列和相位誤差序列與其滯后一階誤差序列的散點圖，從圖3中也可以看到，視電阻率誤差沒有呈現(xiàn)出隨機分布的狀態(tài)，但相位誤差呈現(xiàn)出了明顯的隨機分布特征?；贒W統(tǒng)計量的大地電磁反演方法需要滿足噪聲誤差隨機性的條件，因此需要對視電阻率進行轉(zhuǎn)換。

圖3 視電阻率誤差序列(a)和阻抗相位誤差序列(b)與其滯后一階誤差序列散點圖

在大地電磁反演中，通常會將視電阻率的對數(shù)值作為反演數(shù)據(jù)。若采用給視電阻率對數(shù)值添加絕對誤差的方式，則所加誤差便不具有自相關性。依據(jù)前面的討論，如果給大地電磁阻抗添加相對誤差er, 則應給視電阻率加入2er的相對誤差，這相當于給視電阻率對數(shù)添加絕對誤差log10(1+2er)。其原因解釋如下。

假如給視電阻率ρa添加相對誤差2er，則絕對誤差Δρa=2ρa·er。若視電阻率取對數(shù)，y=log10ρa，則其絕對誤差為：

(18)

還是以高阻薄層模型為例，在視電阻率對數(shù)數(shù)據(jù)中加入log10(1+1%)的隨機絕對誤差。 圖4給出了視電阻率對數(shù)的絕對誤差隨周期的變化關系以及散點圖。由圖4可見，視電阻率對數(shù)的誤差呈現(xiàn)出了隨機分布狀態(tài)。視電阻率對數(shù)誤差的DW統(tǒng)計量DWrho= 1.74和自相關系數(shù)Rrho=0.036，這表明各個頻點視電阻率對數(shù)的絕對誤差之間不存在自相關性。

圖4 視電阻率對數(shù)值的絕對誤差分布情況(a)和視電阻率對數(shù)值的誤差與其滯后一階誤差序列散點圖(b)

考慮到理想情況下，各個頻點大地電磁數(shù)據(jù)的誤差之間應各自獨立，不存在自相關。因此，在大地電磁反演中應將視電阻率對數(shù)值和阻抗相位作為反演數(shù)據(jù)，在添加誤差時應加入相應的絕對誤差。

2 反演算例

2.1 高阻薄層模型

首先考慮一個含有高阻薄層的三層模型(如圖5中黑色虛線)。該模型是根據(jù)文獻[11,20]建立的頁巖油氣高阻薄層模型。第一層的電阻率和厚度分別為30 Ω·m和120 m；中間層為高阻薄層，其電阻率和厚度分別為120 Ω·m和10 m，高阻層下方是電阻率為2.5 Ω·m的均勻半空間。假設21個頻點對數(shù)等間隔地均勻分布在頻率范圍10-3～103Hz之間。在各個頻點的視電阻率對數(shù)和阻抗相位數(shù)據(jù)中，分別加入0.004 3(log10(1+0.01)=0.004 3，相當于給視電阻率添加1%的相對誤差)和0.005弧度的隨機噪聲，從而構(gòu)成反演數(shù)據(jù)。

圖5 反演結(jié)果對比

反演初始模型為100 Ω·m的均勻半空間。正則化因子和DW權重因子的初值分別為λ1=105和λ2=10-4。在迭代過程中，λ1逐漸減小，而λ2逐漸增大。在本例中，λ1的遞減系數(shù)為1.23，λ2的遞增系數(shù)為1.60。圖5展示了經(jīng)過第50次迭代后獲得的反演結(jié)果。從反演結(jié)果中可以看到，高阻薄層厚度及其電阻率均得到了很好的重構(gòu)，其厚度非常接近真實值，其電阻率為151.30 Ω·m，比其真實值偏高。圖6展示了第47次至第50次迭代的反演結(jié)果。從圖6可以看出DW反演方法重構(gòu)高阻薄層的過程。為了比較起見，圖5也給出了傳統(tǒng)高斯-牛頓法反演結(jié)果。高斯-牛頓法反演方法很好地恢復了第一層和第三層的電阻率,但是對高阻薄層幾乎完全沒有反映。

由于大地電磁法存在等值性現(xiàn)象，含有高阻薄層的正演響應曲線與不含高阻薄層的正演響應曲線十分接近，使得大地電磁法對高阻薄層的分辨率極低。海洋大地電磁方法對高阻薄層分辨率低的特性制約了其在高阻薄層探測方面的發(fā)展，而本例則表明含DW項的反演方法可以通過減小反演擬合殘差的自相關性，在一定程度上提高海洋大地電磁方法對高阻薄層的探測能力。

((a)第47次迭代，(b)第48次迭代，(c)第49次迭代，(d)第50次迭代。(a) The 47th iteration， (b) The 48th iteration， (c) The 49th iteration， (d) The 50th iteration.)

圖7給出了DW反演迭代過程中均方根RMS、視電阻率對數(shù)的DW統(tǒng)計量(DWrho)及相位DW統(tǒng)計量(DWpha)的變化情況。從圖7中可以看出，均方根RMS隨著迭代次數(shù)的增加不斷下降，而DWrho和DWpha隨迭代次數(shù)逐漸增大，并逐漸趨近于它們的期望值。

圖7 DW反演迭代過程中RMS與DW值的變化

圖8繪出了DW方法第50次迭代反演模型與真實模型的視電阻率對數(shù)和相位之間的殘差。從圖8可以看到，視電阻率對數(shù)殘差和相位殘差在各頻段上基本呈現(xiàn)出隨機波動的狀態(tài)。由表 1可知，在數(shù)據(jù)個數(shù)為21的情形下，DW統(tǒng)計量位于(1.54, 2.46)內(nèi)時，反演數(shù)據(jù)和反演模型響應之間的殘差將不存在自相關性。圖5所示DW反演模型與真實模型視電阻率對數(shù)殘差和相位殘差的DW統(tǒng)計量，分別為DWrho= 1.638和DWpha= 1.915，它們均位于1.54～2.46之間，這表明DW反演模型的視電阻率對數(shù)殘差和相位殘差之間都不存在自相關性。

圖8 DW反演模型與真實模型的視電阻率殘差(a)和相位殘差(b)

2.2 八層模型反演

圖9中灰色實線所示一個八層地電模型。該模型常被用于驗證大地電磁反演方法的效果[21-22]。先在0.000 1～1 000 Hz之間呈對數(shù)等間距選取40個頻點，經(jīng)正演計算得到40個頻點的大地電磁響應，對視電阻率對數(shù)和相位分別添加0.004 3和0.005弧度的隨機誤差，將其作為反演數(shù)據(jù)。然后分別使用傳統(tǒng)高斯-牛頓反演方法(GN)、OCCAM反演方法和DW反演方法進行反演計算。反演時，初始模型均為100 Ω·m均勻半空間?；谌N反演方法的反演結(jié)果如圖9所示。從反演結(jié)果中可以看出，相較于傳統(tǒng)高斯-牛頓反演方法和OCCAM反演方法，DW方法反演結(jié)果更接近真實模型，對高阻層和低阻層的識別更加準確。DW反演模型視電阻率對數(shù)和相位的DW統(tǒng)計量分別為DWrho=1.208和DWpha= 1.978。由表 1可知，對于40個樣本點，當DW統(tǒng)計量位于(1.60, 2.40)區(qū)間內(nèi)時，反演數(shù)據(jù)和反演模型響應之間的殘差將不存在自相關性。而當DW統(tǒng)計量小于1.39時， 則認為殘差具有正一階自相關。據(jù)此認為，DW反演模型的視電阻率殘差之間具有一定的正自相關性，而相位殘差不具有自相關性。

圖9 DW反演結(jié)果與傳統(tǒng)GN和OCCAM反演結(jié)果對比

圖10 DW反演模型與真實模型的視電阻率殘差(a)與相位殘差(b)

2.3 實測數(shù)據(jù)反演

在本節(jié)中，研究將DW反演方法應用于南黃海實測海洋大地電磁數(shù)據(jù)反演中。南黃海中部隆起區(qū)由于高速海相碳酸鹽巖的屏蔽作用，導致地震勘探方法應用不理想，成為制約海相油氣勘探的關鍵[23]。海洋大地電磁方法由于不受高阻層屏蔽、對低阻層敏感度高的特點，已應用于該海域，初步探明海相碳酸鹽巖層的厚度[24]。由于研究區(qū)水深僅21 m，海水運動(如波浪，潮汐等)使得大地電磁視電阻率和相位在1～10 s周期受到較大干擾，一定程度降低了數(shù)據(jù)質(zhì)量和分辨率，因此，在數(shù)據(jù)處理時需對受海水運動影響較大頻點的大地電磁數(shù)據(jù)進行剔除。下面，給出利用S4站位27個頻點的視電阻率數(shù)據(jù)和19個頻點的相位數(shù)據(jù)的反演結(jié)果。

對S4站位實測MT數(shù)據(jù)進行了DW 反演，并將反演結(jié)果與測點附近的測井電阻率數(shù)據(jù)進行對比(該站位位于CSDP-2井東南方向4 km處[25])。反演時，初始模型為10 Ω·m均勻半空間。圖11為該測點實測MT數(shù)據(jù)的反演結(jié)果。由圖11可以看出，DW反演結(jié)果與傳統(tǒng)高斯-牛頓法反演結(jié)果[24]相比，對六百多米深處高阻界面的識別更為清晰。圖12為DW反演迭代過程中均方根RMS與DW統(tǒng)計量的變化情況。隨著迭代次數(shù)的增加，均方根誤差RMS逐漸下降，而視電阻率與相位的DW值逐漸增大，最后趨于穩(wěn)定。圖13為反演模型的視電阻率和相位與實測數(shù)據(jù)的對比。圖14為它們之間的殘差。由表 1可知，在樣本數(shù)量分別為27與19時，如果DW統(tǒng)計量分別位于區(qū)間(1.56, 2.44)與(1.53，2.47)內(nèi)，則反演數(shù)據(jù)和反演模型響應之間的殘差不存在自相關性。通過對比分析DW值以及擬合曲線可以得到，DWrho= 0.81,相位的DWpha= 0.61, 反演后殘差具有正自相關性，這是因為實測數(shù)據(jù)經(jīng)過處理后的誤差無法滿足隨機性的特點，由于受海水運動、儀器噪聲等影響，實測數(shù)據(jù)的誤差具有較強的正自相關性。但是，較于傳統(tǒng)的高斯-牛頓法反演，DW反演結(jié)果對淺層高阻層與620 m處的高阻層界面識別更清晰。

(測井和GN反演結(jié)果來自文獻[24]。The GN inversion results are from reference[24].)

圖12 DW反演迭代過程中RMS與DW值的變化

圖13 反演模型響應(實線)與實測數(shù)據(jù)的視電阻率(a)和相位曲線(b)對比

圖14 DW反演模型響應的視電阻率殘差(a)和相位殘差(b)

3 結(jié)語

在MT反演的目標函數(shù)中加入Durbin-Watson(DW)統(tǒng)計量，可以降低反演模型響應與觀測數(shù)據(jù)之間殘差的自相關性，在一定程度上避免反演結(jié)果中出現(xiàn)的擬合差較小而反演模型響應偏離觀測數(shù)據(jù)的情況，從而優(yōu)化反演結(jié)果，降低反演的多解性。

本文通過兩個層狀模型模擬數(shù)據(jù)反演，驗證了DW反演方法的有效性，并與傳統(tǒng)高斯-牛頓法或OCCAM反演方法進行了對比，結(jié)果表明基于DW統(tǒng)計量的反演方法提高了對高阻薄層的反演分辨能力。南黃海實測大地電磁數(shù)據(jù)反演結(jié)果表明，該方法能夠較準確地重構(gòu)高阻層界面，具有良好的應用前景和效果。

針對目標函數(shù)中各項權重系數(shù)的選取策略以及如何將DW反演方法推廣到二維和三維大地電磁反演中，有待下一步繼續(xù)研究。