楊 雯, 李予國,2??, 段雙敏, 韓 波, 羅 鳴,2
(1. 中國海洋大學(xué)海洋地球科學(xué)學(xué)院, 山東 青島 266100; 2. 青島海洋科學(xué)與技術(shù)試點(diǎn)國家實(shí)驗(yàn)室 海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室, 山東 青島 266237; 3. 中國地質(zhì)大學(xué)(武漢) 地質(zhì)調(diào)查研究院,湖北 武漢 430074)
大地電磁法(Magnetotelluric,MT)是利用天然交變電磁信號(hào)研究地球內(nèi)部電性結(jié)構(gòu)的一種地球物理方法,現(xiàn)已廣泛應(yīng)用于油氣資源勘探和俯沖帶、地殼與地幔電性結(jié)構(gòu)研究中[1-3]。大地電磁反演是獲得地下介質(zhì)電阻率信息的重要方法。然而,MT反演存在多解性,如何減少M(fèi)T反演的多解性一直是重要的研究方向。自從Karl Pearson提出卡方檢驗(yàn)以來,在地球物理反演中便一直使用均方根擬合差(RMS misfit)衡量反演模型響應(yīng)與實(shí)測(cè)數(shù)據(jù)的擬合程度[4]。
1980年代末以來,為了改善大地電磁反演問題的不適定性,吉洪諾夫(Tikhonov)正則化方法開始被應(yīng)用于MT反演中,并取得了迅速發(fā)展[5-7]。目前,常用的大地電磁反演方法包括OCCAM反演方法、非線性共軛梯度法(NLCG)、擬牛頓法和高斯-牛頓法等[7-10]。這些反演方法關(guān)于反演模型響應(yīng)與觀測(cè)數(shù)據(jù)擬合程度的描述大都依賴于目標(biāo)函數(shù)中的均方根擬合差,卻沒有考慮反演模型響應(yīng)與觀測(cè)數(shù)據(jù)之差(常稱為殘差)的統(tǒng)計(jì)特征,這可能導(dǎo)致反演模型響應(yīng)形態(tài)與觀測(cè)數(shù)據(jù)形態(tài)出現(xiàn)偏離較大的情況,尤其是當(dāng)觀測(cè)數(shù)據(jù)存在較大誤差時(shí),甚至可能得到錯(cuò)誤的反演模型[11]。于是用于衡量實(shí)測(cè)數(shù)據(jù)與反演模型響應(yīng)之間關(guān)系的其它統(tǒng)計(jì)量開始被應(yīng)用于大地電磁反演中。Smith 和 Booker將斯皮爾曼(Spearman)相關(guān)系數(shù)應(yīng)用于大地電磁(MT)反演結(jié)果評(píng)價(jià)中[12]。Jones利用DW統(tǒng)計(jì)量衡量觀測(cè)數(shù)據(jù)與模型響應(yīng)之間殘差的自相關(guān)性[11],并將其應(yīng)用于加拿大Great Slave Lake地區(qū)大地電磁測(cè)深數(shù)據(jù)反演中[13]。將統(tǒng)計(jì)學(xué)度量應(yīng)用于MT反演中,一方面可以評(píng)價(jià)反演結(jié)果的合理性,另一方面可以為反演算法的改進(jìn)提供依據(jù),以便進(jìn)一步優(yōu)化反演方法。但是目前關(guān)于這方面的研究還不多見。本文將DW統(tǒng)計(jì)量應(yīng)用于高斯-牛頓大地電磁一維反演中,通過修改目標(biāo)函數(shù),實(shí)現(xiàn)了基于DW統(tǒng)計(jì)量的大地電磁一維反演方法(以下簡稱為DW反演方法)。通過模型算例分析了該方法的有效性,并與傳統(tǒng)的高斯-牛頓反演結(jié)果進(jìn)行了對(duì)比。
杜賓-沃森(Durbin-Watson,簡稱DW)統(tǒng)計(jì)方法是用于檢驗(yàn)回歸分析中殘差一階自相關(guān)性的一種方法,
該方法常用于檢驗(yàn)計(jì)量經(jīng)濟(jì)學(xué)中的自相關(guān)問題[14-15]。一階自相關(guān)DW統(tǒng)計(jì)量的定義式為[11,16-17]:
(1)
圖1 序列自相關(guān)性與DW統(tǒng)計(jì)量的關(guān)系
表1 部分DW統(tǒng)計(jì)量臨界值分布表(k=2)
通常地,DW統(tǒng)計(jì)量與自相關(guān)系數(shù)R具有如下關(guān)系[16,18]:
DW≈2(1-R)。
(2)
式中,R為一階自相關(guān)系數(shù),其表達(dá)式為:
(3)
一階自相關(guān)系數(shù)R的取值范圍為(-1,1)。通常認(rèn)為,當(dāng)|R|<0.2時(shí),殘差序列的自相關(guān)性極弱。
大地電磁反演問題實(shí)質(zhì)上是目標(biāo)函數(shù)的最優(yōu)化問題?;诩橹Z夫正則化反演思想,反演目標(biāo)函數(shù)中通常包含數(shù)據(jù)擬合差項(xiàng)和正則化項(xiàng)[19]:
φ(m)=φd(m)+λ1φm(m)。
(4)
式中:m=[m1,m2,···,mM]T為模型參數(shù)向量,由反演模型中各地層的電阻率和厚度構(gòu)成;φd(m)為觀測(cè)數(shù)據(jù)與理論模型響應(yīng)的擬合差,可表示為:
(5)
(6)
式中:Wm為拉普拉斯算子的差分近似;mref為參考模型。
在傳統(tǒng)的反演方法中,反演目標(biāo)函數(shù)僅考慮了反演模型響應(yīng)與觀測(cè)數(shù)據(jù)的擬合差,但沒有考慮其殘差的統(tǒng)計(jì)學(xué)特征。本文將DW統(tǒng)計(jì)量加入反演目標(biāo)函數(shù)中,通過反演傾向于尋找在擬合差和殘差自相關(guān)性兩個(gè)方面都可以接受的模型。經(jīng)過修改后的反演目標(biāo)函數(shù)具有如下形式[11]:
φ(m)=φd(m)+λ1φm(m)+λ2(DW-2)2。
(7)
上式中等號(hào)右端第三項(xiàng)是將DW統(tǒng)計(jì)量歸一化為零。在反演中,正則化因子λ1從一個(gè)較大的初值開始,隨迭代次數(shù)的增加逐漸減小;而與此相反,對(duì)于權(quán)重系數(shù)λ2,則先給定一個(gè)較小的初值,在迭代過程中緩慢遞增。λ1遞減和λ2遞增的系數(shù)通常位于(1,2)之間。采用這種策略的基本思想是,在迭代反演的前期,注重?cái)?shù)據(jù)擬合和模型約束,而在反演的后期通過增加DW項(xiàng)的權(quán)重,減少實(shí)測(cè)數(shù)據(jù)與反演模型響應(yīng)之間殘差的自相關(guān)性。
本文采用高斯-牛頓最優(yōu)化方法求解目標(biāo)函數(shù)式(7)的極小值。迭代反演過程中,模型更新量為:
(8)
其中g(shù)k和Hk分別為第k次迭代時(shí)目標(biāo)函數(shù)φ(m)的梯度向量和Hessian矩陣,其表達(dá)式為:
(9)
(10)
在傳統(tǒng)的高斯-牛頓反演方法中,gk和Hk分別僅由式(9)和式(10)右端的前兩項(xiàng)構(gòu)成,而除該兩項(xiàng)以外的其它項(xiàng)反映了DW約束項(xiàng)對(duì)梯度向量和Hessian矩陣的貢獻(xiàn)。其中:
(11)
Hdw=2QJTKTKJ-2PQ2LLT+8Q2JTKTKe(LE)T+
8PQ3LE(LE)T。
(12)
大地電磁阻抗Z可以寫成下列指數(shù)形式
Z=reiφ。
(13)
由式(13)可以得到阻抗微小變化dZ的表達(dá)式
dZ=dreiφ+reiφ·idφ。
(14)
式(14)除以式(13),得
(15)
上式表明,大地電磁阻抗Z的相對(duì)變化相當(dāng)于其幅值的相對(duì)變化和相位的絕對(duì)變化。大多數(shù)時(shí)間序列分析表明,大地電磁阻抗的實(shí)部和虛部通常具有相等的誤差。于是,有
(16)
這意味著,對(duì)于大地電磁阻抗Z,其相位的絕對(duì)誤差與幅值的相對(duì)誤差相對(duì)應(yīng)。也就是說,如果給大地電磁阻抗幅值添加相對(duì)誤差er, 則對(duì)應(yīng)的應(yīng)該給阻抗相位添加er弧度或er·180/π角度的絕對(duì)誤差。
(17)
大地電磁理論模型合成數(shù)據(jù)反演中,通常會(huì)給正演合成數(shù)據(jù)加入一定的隨機(jī)噪聲,并將其結(jié)果作為反演數(shù)據(jù)。然而,若直接給視電阻率添加相對(duì)誤差,可能會(huì)使得所添加的誤差具有自相關(guān)性,以2.1節(jié)將要討論的三層高阻薄層模型為例,對(duì)此進(jìn)行說明。
首先,在頻率范圍10-3~103Hz之間選取21個(gè)對(duì)數(shù)等間隔的頻點(diǎn),正演計(jì)算得到三層高阻薄層模型21個(gè)頻點(diǎn)的視電阻率和阻抗相位。然后,在視電阻率數(shù)據(jù)中加入1%的隨機(jī)相對(duì)誤差噪聲,相應(yīng)地在阻抗相位數(shù)據(jù)中加入0.286 5°的隨機(jī)絕對(duì)誤差。用這樣的方法得到的視電阻率誤差和相位誤差隨周期的變化關(guān)系如圖2所示。
圖2 視電阻率誤差(a)和阻抗相位誤差(b)分布特征
從圖2中可以看到,阻抗相位誤差在周期范圍(10-3~103s)內(nèi)是隨機(jī)分布的。阻抗相位誤差的DW統(tǒng)計(jì)量DWpha=1.74和自相關(guān)系數(shù)Rpha=0.036,這表明阻抗相位誤差不具有自相關(guān)性。長周期視電阻率的誤差是在零線附近上下隨機(jī)分布的,而短周期視電阻率的誤差并非是隨機(jī)分布的。視電阻率誤差的DW統(tǒng)計(jì)量DWrho為 0.81,自相關(guān)系數(shù)Rrho為0.568,這表明視電阻率誤差具有一定的自相關(guān)性。圖3分別繪出了視電阻率誤差序列和相位誤差序列與其滯后一階誤差序列的散點(diǎn)圖,從圖3中也可以看到,視電阻率誤差沒有呈現(xiàn)出隨機(jī)分布的狀態(tài),但相位誤差呈現(xiàn)出了明顯的隨機(jī)分布特征?;贒W統(tǒng)計(jì)量的大地電磁反演方法需要滿足噪聲誤差隨機(jī)性的條件,因此需要對(duì)視電阻率進(jìn)行轉(zhuǎn)換。
圖3 視電阻率誤差序列(a)和阻抗相位誤差序列(b)與其滯后一階誤差序列散點(diǎn)圖
在大地電磁反演中,通常會(huì)將視電阻率的對(duì)數(shù)值作為反演數(shù)據(jù)。若采用給視電阻率對(duì)數(shù)值添加絕對(duì)誤差的方式,則所加誤差便不具有自相關(guān)性。依據(jù)前面的討論,如果給大地電磁阻抗添加相對(duì)誤差er, 則應(yīng)給視電阻率加入2er的相對(duì)誤差,這相當(dāng)于給視電阻率對(duì)數(shù)添加絕對(duì)誤差log10(1+2er)。其原因解釋如下。
假如給視電阻率ρa(bǔ)添加相對(duì)誤差2er,則絕對(duì)誤差Δρa(bǔ)=2ρa(bǔ)·er。若視電阻率取對(duì)數(shù),y=log10ρa(bǔ),則其絕對(duì)誤差為:
(18)
還是以高阻薄層模型為例,在視電阻率對(duì)數(shù)數(shù)據(jù)中加入log10(1+1%)的隨機(jī)絕對(duì)誤差。 圖4給出了視電阻率對(duì)數(shù)的絕對(duì)誤差隨周期的變化關(guān)系以及散點(diǎn)圖。由圖4可見,視電阻率對(duì)數(shù)的誤差呈現(xiàn)出了隨機(jī)分布狀態(tài)。視電阻率對(duì)數(shù)誤差的DW統(tǒng)計(jì)量DWrho= 1.74和自相關(guān)系數(shù)Rrho=0.036,這表明各個(gè)頻點(diǎn)視電阻率對(duì)數(shù)的絕對(duì)誤差之間不存在自相關(guān)性。
圖4 視電阻率對(duì)數(shù)值的絕對(duì)誤差分布情況(a)和視電阻率對(duì)數(shù)值的誤差與其滯后一階誤差序列散點(diǎn)圖(b)
考慮到理想情況下,各個(gè)頻點(diǎn)大地電磁數(shù)據(jù)的誤差之間應(yīng)各自獨(dú)立,不存在自相關(guān)。因此,在大地電磁反演中應(yīng)將視電阻率對(duì)數(shù)值和阻抗相位作為反演數(shù)據(jù),在添加誤差時(shí)應(yīng)加入相應(yīng)的絕對(duì)誤差。
首先考慮一個(gè)含有高阻薄層的三層模型(如圖5中黑色虛線)。該模型是根據(jù)文獻(xiàn)[11,20]建立的頁巖油氣高阻薄層模型。第一層的電阻率和厚度分別為30 Ω·m和120 m;中間層為高阻薄層,其電阻率和厚度分別為120 Ω·m和10 m,高阻層下方是電阻率為2.5 Ω·m的均勻半空間。假設(shè)21個(gè)頻點(diǎn)對(duì)數(shù)等間隔地均勻分布在頻率范圍10-3~103Hz之間。在各個(gè)頻點(diǎn)的視電阻率對(duì)數(shù)和阻抗相位數(shù)據(jù)中,分別加入0.004 3(log10(1+0.01)=0.004 3,相當(dāng)于給視電阻率添加1%的相對(duì)誤差)和0.005弧度的隨機(jī)噪聲,從而構(gòu)成反演數(shù)據(jù)。
圖5 反演結(jié)果對(duì)比
反演初始模型為100 Ω·m的均勻半空間。正則化因子和DW權(quán)重因子的初值分別為λ1=105和λ2=10-4。在迭代過程中,λ1逐漸減小,而λ2逐漸增大。在本例中,λ1的遞減系數(shù)為1.23,λ2的遞增系數(shù)為1.60。圖5展示了經(jīng)過第50次迭代后獲得的反演結(jié)果。從反演結(jié)果中可以看到,高阻薄層厚度及其電阻率均得到了很好的重構(gòu),其厚度非常接近真實(shí)值,其電阻率為151.30 Ω·m,比其真實(shí)值偏高。圖6展示了第47次至第50次迭代的反演結(jié)果。從圖6可以看出DW反演方法重構(gòu)高阻薄層的過程。為了比較起見,圖5也給出了傳統(tǒng)高斯-牛頓法反演結(jié)果。高斯-牛頓法反演方法很好地恢復(fù)了第一層和第三層的電阻率,但是對(duì)高阻薄層幾乎完全沒有反映。
由于大地電磁法存在等值性現(xiàn)象,含有高阻薄層的正演響應(yīng)曲線與不含高阻薄層的正演響應(yīng)曲線十分接近,使得大地電磁法對(duì)高阻薄層的分辨率極低。海洋大地電磁方法對(duì)高阻薄層分辨率低的特性制約了其在高阻薄層探測(cè)方面的發(fā)展,而本例則表明含DW項(xiàng)的反演方法可以通過減小反演擬合殘差的自相關(guān)性,在一定程度上提高海洋大地電磁方法對(duì)高阻薄層的探測(cè)能力。
((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、視電阻率對(duì)數(shù)的DW統(tǒng)計(jì)量(DWrho)及相位DW統(tǒng)計(jì)量(DWpha)的變化情況。從圖7中可以看出,均方根RMS隨著迭代次數(shù)的增加不斷下降,而DWrho和DWpha隨迭代次數(shù)逐漸增大,并逐漸趨近于它們的期望值。
圖7 DW反演迭代過程中RMS與DW值的變化
圖8繪出了DW方法第50次迭代反演模型與真實(shí)模型的視電阻率對(duì)數(shù)和相位之間的殘差。從圖8可以看到,視電阻率對(duì)數(shù)殘差和相位殘差在各頻段上基本呈現(xiàn)出隨機(jī)波動(dòng)的狀態(tài)。由表 1可知,在數(shù)據(jù)個(gè)數(shù)為21的情形下,DW統(tǒng)計(jì)量位于(1.54, 2.46)內(nèi)時(shí),反演數(shù)據(jù)和反演模型響應(yīng)之間的殘差將不存在自相關(guān)性。圖5所示DW反演模型與真實(shí)模型視電阻率對(duì)數(shù)殘差和相位殘差的DW統(tǒng)計(jì)量,分別為DWrho= 1.638和DWpha= 1.915,它們均位于1.54~2.46之間,這表明DW反演模型的視電阻率對(duì)數(shù)殘差和相位殘差之間都不存在自相關(guān)性。
圖8 DW反演模型與真實(shí)模型的視電阻率殘差(a)和相位殘差(b)
圖9中灰色實(shí)線所示一個(gè)八層地電模型。該模型常被用于驗(yàn)證大地電磁反演方法的效果[21-22]。先在0.000 1~1 000 Hz之間呈對(duì)數(shù)等間距選取40個(gè)頻點(diǎn),經(jīng)正演計(jì)算得到40個(gè)頻點(diǎn)的大地電磁響應(yīng),對(duì)視電阻率對(duì)數(shù)和相位分別添加0.004 3和0.005弧度的隨機(jī)誤差,將其作為反演數(shù)據(jù)。然后分別使用傳統(tǒng)高斯-牛頓反演方法(GN)、OCCAM反演方法和DW反演方法進(jìn)行反演計(jì)算。反演時(shí),初始模型均為100 Ω·m均勻半空間?;谌N反演方法的反演結(jié)果如圖9所示。從反演結(jié)果中可以看出,相較于傳統(tǒng)高斯-牛頓反演方法和OCCAM反演方法,DW方法反演結(jié)果更接近真實(shí)模型,對(duì)高阻層和低阻層的識(shí)別更加準(zhǔn)確。DW反演模型視電阻率對(duì)數(shù)和相位的DW統(tǒng)計(jì)量分別為DWrho=1.208和DWpha= 1.978。由表 1可知,對(duì)于40個(gè)樣本點(diǎn),當(dāng)DW統(tǒng)計(jì)量位于(1.60, 2.40)區(qū)間內(nèi)時(shí),反演數(shù)據(jù)和反演模型響應(yīng)之間的殘差將不存在自相關(guān)性。而當(dāng)DW統(tǒng)計(jì)量小于1.39時(shí), 則認(rèn)為殘差具有正一階自相關(guān)。據(jù)此認(rèn)為,DW反演模型的視電阻率殘差之間具有一定的正自相關(guān)性,而相位殘差不具有自相關(guān)性。
圖9 DW反演結(jié)果與傳統(tǒng)GN和OCCAM反演結(jié)果對(duì)比
圖10 DW反演模型與真實(shí)模型的視電阻率殘差(a)與相位殘差(b)
在本節(jié)中,研究將DW反演方法應(yīng)用于南黃海實(shí)測(cè)海洋大地電磁數(shù)據(jù)反演中。南黃海中部隆起區(qū)由于高速海相碳酸鹽巖的屏蔽作用,導(dǎo)致地震勘探方法應(yīng)用不理想,成為制約海相油氣勘探的關(guān)鍵[23]。海洋大地電磁方法由于不受高阻層屏蔽、對(duì)低阻層敏感度高的特點(diǎn),已應(yīng)用于該海域,初步探明海相碳酸鹽巖層的厚度[24]。由于研究區(qū)水深僅21 m,海水運(yùn)動(dòng)(如波浪,潮汐等)使得大地電磁視電阻率和相位在1~10 s周期受到較大干擾,一定程度降低了數(shù)據(jù)質(zhì)量和分辨率,因此,在數(shù)據(jù)處理時(shí)需對(duì)受海水運(yùn)動(dòng)影響較大頻點(diǎn)的大地電磁數(shù)據(jù)進(jìn)行剔除。下面,給出利用S4站位27個(gè)頻點(diǎn)的視電阻率數(shù)據(jù)和19個(gè)頻點(diǎn)的相位數(shù)據(jù)的反演結(jié)果。
對(duì)S4站位實(shí)測(cè)MT數(shù)據(jù)進(jìn)行了DW 反演,并將反演結(jié)果與測(cè)點(diǎn)附近的測(cè)井電阻率數(shù)據(jù)進(jìn)行對(duì)比(該站位位于CSDP-2井東南方向4 km處[25])。反演時(shí),初始模型為10 Ω·m均勻半空間。圖11為該測(cè)點(diǎn)實(shí)測(cè)MT數(shù)據(jù)的反演結(jié)果。由圖11可以看出,DW反演結(jié)果與傳統(tǒng)高斯-牛頓法反演結(jié)果[24]相比,對(duì)六百多米深處高阻界面的識(shí)別更為清晰。圖12為DW反演迭代過程中均方根RMS與DW統(tǒng)計(jì)量的變化情況。隨著迭代次數(shù)的增加,均方根誤差RMS逐漸下降,而視電阻率與相位的DW值逐漸增大,最后趨于穩(wěn)定。圖13為反演模型的視電阻率和相位與實(shí)測(cè)數(shù)據(jù)的對(duì)比。圖14為它們之間的殘差。由表 1可知,在樣本數(shù)量分別為27與19時(shí),如果DW統(tǒng)計(jì)量分別位于區(qū)間(1.56, 2.44)與(1.53,2.47)內(nèi),則反演數(shù)據(jù)和反演模型響應(yīng)之間的殘差不存在自相關(guān)性。通過對(duì)比分析DW值以及擬合曲線可以得到,DWrho= 0.81,相位的DWpha= 0.61, 反演后殘差具有正自相關(guān)性,這是因?yàn)閷?shí)測(cè)數(shù)據(jù)經(jīng)過處理后的誤差無法滿足隨機(jī)性的特點(diǎn),由于受海水運(yùn)動(dòng)、儀器噪聲等影響,實(shí)測(cè)數(shù)據(jù)的誤差具有較強(qiáng)的正自相關(guān)性。但是,較于傳統(tǒng)的高斯-牛頓法反演,DW反演結(jié)果對(duì)淺層高阻層與620 m處的高阻層界面識(shí)別更清晰。
(測(cè)井和GN反演結(jié)果來自文獻(xiàn)[24]。The GN inversion results are from reference[24].)
圖12 DW反演迭代過程中RMS與DW值的變化
圖13 反演模型響應(yīng)(實(shí)線)與實(shí)測(cè)數(shù)據(jù)的視電阻率(a)和相位曲線(b)對(duì)比
圖14 DW反演模型響應(yīng)的視電阻率殘差(a)和相位殘差(b)
在MT反演的目標(biāo)函數(shù)中加入Durbin-Watson(DW)統(tǒng)計(jì)量,可以降低反演模型響應(yīng)與觀測(cè)數(shù)據(jù)之間殘差的自相關(guān)性,在一定程度上避免反演結(jié)果中出現(xiàn)的擬合差較小而反演模型響應(yīng)偏離觀測(cè)數(shù)據(jù)的情況,從而優(yōu)化反演結(jié)果,降低反演的多解性。
本文通過兩個(gè)層狀模型模擬數(shù)據(jù)反演,驗(yàn)證了DW反演方法的有效性,并與傳統(tǒng)高斯-牛頓法或OCCAM反演方法進(jìn)行了對(duì)比,結(jié)果表明基于DW統(tǒng)計(jì)量的反演方法提高了對(duì)高阻薄層的反演分辨能力。南黃海實(shí)測(cè)大地電磁數(shù)據(jù)反演結(jié)果表明,該方法能夠較準(zhǔn)確地重構(gòu)高阻層界面,具有良好的應(yīng)用前景和效果。
針對(duì)目標(biāo)函數(shù)中各項(xiàng)權(quán)重系數(shù)的選取策略以及如何將DW反演方法推廣到二維和三維大地電磁反演中,有待下一步繼續(xù)研究。