二維與三維大地電磁反演的敏感度研究

孟綠汀， 張慧茜， 黃清華*

1 北京大學(xué)地球與空間科學(xué)學(xué)院地球物理系， 北京 100871 2 河北紅山巨厚沉積與地震災(zāi)害國家野外科學(xué)觀測站(北京大學(xué))， 北京 100871

0 引言

大地電磁法以天然交變電磁場為場源，通過在地表觀測電場強(qiáng)度和磁場強(qiáng)度來研究地球內(nèi)部電性結(jié)構(gòu).由于其對低阻層的高敏感度以及深達(dá)軟流圈的探測深度，大地電磁法被廣泛應(yīng)用于資源勘探(Queralt et al., 2007)、巖石圈結(jié)構(gòu)探測(Zhang et al., 2016)、地震孕震環(huán)境研究(葉濤等，2021)等方面.該方法假設(shè)電磁波以平面波形式透射入地表，利用電磁波的探測深度隨頻率變化的特征來重構(gòu)地下結(jié)構(gòu)的電阻率模型(陳樂壽和王光鄂，1990).由于大地電磁反演為非線性反演且耗時較長， 大多數(shù)實(shí)際研究僅關(guān)注反演給出的物理參數(shù)模型并結(jié)合地質(zhì)資料對其進(jìn)行解釋，易于忽視反演模型的可靠性與分辨率等信息.事實(shí)上由于反演問題的多解性，反演模型可能與真實(shí)地質(zhì)結(jié)構(gòu)存在一定差異.此外，由于電磁擴(kuò)散場的物理本質(zhì)，大地電磁法重建的深部電性結(jié)構(gòu)往往還具有分辨率不足的問題.定量估計(jì)大地電磁反演模型的可靠性和分辨率，能為后續(xù)精細(xì)地質(zhì)解釋提供依據(jù).因此，如何構(gòu)建經(jīng)濟(jì)高效、易于實(shí)現(xiàn)的反演結(jié)果評價方法，成為受到廣泛關(guān)注的問題.

反演模型的可靠性可以從多種角度評價，評價方法各有優(yōu)劣.采用解析方法計(jì)算趨膚深度(陳小斌等, 2007, 2019; Borah and Patro, 2019)和最大探測深度(Spies, 1989; Huang, 2005; 肖調(diào)杰等, 2015)簡便易行，但只適用于一維層狀模型，難以擴(kuò)展到二維、三維場景.三維反演通常因?yàn)槟Ｐ蛥?shù)過多導(dǎo)致分辨率矩陣、協(xié)方差矩陣難以直接計(jì)算與存儲(Menke, 2015)，因此衍生了許多間接計(jì)算的方法.非線性敏感度測試(Dong et al., 2014; Zhang et al., 2016)、探測深度指數(shù)(Oldenburg and Li, 1999; Oldenborger et al., 2007; Carrière et al., 2017)和目標(biāo)函數(shù)極值分析(Jackson, 1976; Meju, 2009; Kalscheuer et al., 2010; de Wit et al., 2012)等間接方法常被應(yīng)用于評價電磁反演結(jié)果.此類方法需要額外設(shè)計(jì)模型或重新求解最優(yōu)化問題，所以很少應(yīng)用于三維場景.基于高斯過程(Astic et al., 2020; Olierook et al., 2021) 或貝葉斯反演(Minsley, 2011; Xiang et al., 2018; 周思杰和黃清華, 2018; Ray and Myer, 2019; Manassero et al., 2020; Seillé and Visser, 2020; Blatter et al., 2021; Peng et al., 2021)等概率類反演方法可同時獲得解的分布與不確定性.盡管可以采取高維奇異值分解、主成分分析等方法(Tompkins et al., 2011; Tompkins, 2012; Fernández-Martínez, 2015; Fernández-Martínez et al., 2019; Grana et al., 2019)對模型參數(shù)和數(shù)據(jù)參數(shù)進(jìn)行降維處理，但耗時的正演計(jì)算和低效的采樣算法依舊限制概率類方法與自抽樣方法(Schnaidt and Heinson, 2015)在三維場景的應(yīng)用.

隨著算法的完善與計(jì)算力的提升，三維大地電磁反演技術(shù)已成為主流，但受限于野外工作環(huán)境、施工成本等因素，往往無法得到陣列式分布的大地電磁臺站，臺陣分布為多條長剖面相互交叉，這種情況下二維反演技術(shù)依然具有重要的應(yīng)用價值 (Habibian Dehkordi et al., 2019; Nagarjuna et al., 2021).前人通過對三維數(shù)據(jù)進(jìn)行二維反演揭示了二維反演潛在的局限性(Ledo et al., 2002; Ledo, 2006)，同時也利用三維反演解釋二維剖面數(shù)據(jù)(Siripunvaraporn et al., 2005; 林昌洪等, 2011).盡管三維解釋技術(shù)具有一定的優(yōu)越性，但如何針對實(shí)際情況對比和分析二維與三維反演模型的差異依然具有較強(qiáng)的研究意義.本文嘗試發(fā)展一種基于敏感度矩陣的反演模型評價方法，該方法的運(yùn)算規(guī)模較小，且能夠定性分析二維與三維反演結(jié)果的可信度.此外，大地電磁阻抗張量的對角元素、非對角元素與傾子在反演中對模型重建具有不同的貢獻(xiàn)(Kiyan et al., 2014)，不同分量在數(shù)據(jù)采集中的信噪比及反演設(shè)定的誤差限也存在一定差異，從理論上探討各阻抗分量的敏感度對于實(shí)測資料解釋具有一定指導(dǎo)意義.本文從求解三維大地電磁正演模型的伴隨問題來獲取完整的敏感度矩陣，并基于敏感度矩陣定量分析大地電磁數(shù)據(jù)不同阻抗分量敏感度，形成了一套經(jīng)濟(jì)高效的反演結(jié)果評估方法.該方法能夠獲得異常體的邊界、垂向分辨率等信息.研究還利用合成數(shù)據(jù)定量討論非陣列式數(shù)據(jù)二維與三維反演的有效探測深度，為二維反演與三維反演結(jié)果的對比提供了理論依據(jù).

1 敏感度矩陣計(jì)算方法

如式(1)所示，敏感度矩陣表征數(shù)據(jù)對模型變化的響應(yīng)，能夠展示數(shù)據(jù)對模型的約束能力：

(1)

其中m為模型參數(shù)，f(m)為模型到數(shù)據(jù)的映射函數(shù).

本文基于三維交錯網(wǎng)格有限差分法(Egbert and Kelbert, 2012) 求解正演問題的伴隨問題，間接得到完整敏感度矩陣計(jì)算公式.離散形式的頻率域電磁偏微分方程一般表示為式(2)：

Sme=b,

(2)

其中b向量為邊界條件與場源項(xiàng)，Sm為稀疏系數(shù)矩陣，e為交錯網(wǎng)格中主網(wǎng)格上的電場值.由于磁場分量可由電場的旋度給出，此種情況下計(jì)算阻抗張量的函數(shù)與模型參數(shù)無直接關(guān)系，大地電磁正演過程可用式(3)描述：

f(m)=d(e(m)),

(3)

其中d為從離散電場值到阻抗張量、傾子的映射函數(shù).對式(3)關(guān)于m求偏導(dǎo)可得敏感度矩陣：

(4)

將式(4)寫作矩陣形式(5)：

J=LF,

(5)

Sm0F=P,

(6)

(7)

(8)

其中diag為提取矩陣對角線元素，Wd為數(shù)據(jù)協(xié)方差矩陣，表示對敏感度矩陣做關(guān)于數(shù)據(jù)誤差的歸一化，E為Nd×Nd(數(shù)據(jù)個數(shù))的對角系數(shù)矩陣.E的對角線元素取值區(qū)間為0到1，取值代表該數(shù)據(jù)分量對模型影響的權(quán)重.當(dāng)E為單位矩陣時，表示所有數(shù)據(jù)對模型參數(shù)影響的總和.通過至多Nd次計(jì)算可得到完整敏感度矩陣.

2 模擬結(jié)果

2.1 合成模型

合成數(shù)據(jù)測試使用的模型如圖1所示，在100 Ωm的均勻半空間內(nèi)交錯排列有6個1000 Ωm的高阻異常(藍(lán)色)和6個10 Ωm的低阻異常(紅色)，厚度分別為10 km、25 km、60 km，中心深度12.5 km、32 km、92.5 km.模型網(wǎng)格劃分為60×60×60個，中心計(jì)算區(qū)域網(wǎng)格為40×40×60個，水平網(wǎng)格尺寸5 km×5 km，邊界區(qū)域由10個以1.5倍遞增的網(wǎng)格組成，縱向網(wǎng)格第一層為20 m，以1.2倍向下遞增.

圖1 (a)合成模型示意圖藍(lán)色為1000 Ωm 高阻異常體，紅色為10 Ωm 低阻異常體.黑色實(shí)線(X=40 km)和白色實(shí)線(X=60 km)為本文研究涉及的垂向切片的位置，黑色三角表示數(shù)據(jù)接收點(diǎn)，白色三角表示模型中心處數(shù)據(jù)接收點(diǎn).(b)沿黑色實(shí)線的垂向切片，其中C1、C2和C3表示低阻異常體，R1表示高阻異常體.Fig.1 (a) Synthetic model composed of anomaliesBlue and red colors indicate 1000 Ωm high and 10 Ωm low electrical resistivity anomalies, respectively. The solid black line (X=40 km) and solid white line (X=60 km) indicate the location of the vertical slices in this study. The black triangles indicate the stations and the white triangles indicates the locations of the station in model center. (b) The vertical slice along the solid black line where C1—C3 represent the three conductive anomalies and R1 represents the resistivity anomalies.

2.2 不同分量的敏感度

使用2.1節(jié)的合成模型，計(jì)算所有臺站阻抗張量非對角元素、阻抗張量對角元素和傾子矢量的敏感度矩陣，周期范圍0.2～8000 s，等對數(shù)間隔頻點(diǎn)25個點(diǎn).圖2為部分代表性周期下阻抗張量非對角元素Zxy+Zyx的敏感度.在短周期時，高敏感度分布于淺層，集中在測點(diǎn)附近，這是因?yàn)榇蟮仉姶欧M足擴(kuò)散方程，短周期對應(yīng)的趨膚深度較小.長周期時模型深部敏感度增加，不同周期下敏感度峰值處于不同的深度，與大地電磁反演常規(guī)認(rèn)識一致.隨著周期增加，靈敏度高的主要分布區(qū)間的厚度增大，一方面是因?yàn)槟Ｐ途W(wǎng)格尺寸隨深度遞增，導(dǎo)致靈敏度高的主要分布區(qū)間增厚.另一方面即使使用等間距網(wǎng)格，也會存在逐漸增厚的高靈敏度主要分布區(qū)間.表明大地電磁法對深部結(jié)構(gòu)的分辨率有限.非對角元素的敏感度主要分布在橫向邊界附近，且低阻異常體附近的敏感度較大.

圖2 不同周期下所有臺站Zxy+Zyx敏感度沿X=40 km的垂向切片黑色虛線指示異常體位置.Fig.2 Vertical slices of all stations Zxy+Zyx sensitivity along X=40 km at different periodsThe dashed black lines indicate the locations of the anomalies.

如圖3所示，阻抗張量對角元素不僅在低阻異常體附近有高敏感度，在高阻異常體邊界處也有高敏感度，可以重建異常體邊界.圖4中傾子矢量的敏感度集中于異常體邊界處，在異常體內(nèi)部很小.值得注意的是，阻抗張量非對角元素的敏感度大于阻抗張量對角元素和傾子矢量的敏感度，表明權(quán)重相同時，非對角元素多用于重建研究區(qū)域內(nèi)大尺寸結(jié)構(gòu)，而對角元素和傾子矢量對部分邊界敏感度高，在重建異常體的邊界等精細(xì)結(jié)構(gòu)上具有一定貢獻(xiàn).

圖3 不同周期下所有臺站Zxx+Zyy敏感度沿X=40 km的垂向切片黑色虛線指示異常體位置.Fig.3 Vertical slices of all stations Zxx+Zyy sensitivity along X=40 km at different periodsThe dashed black lines indicate the locations of the anomalies.

圖4 不同周期下所有臺站Tx+Ty敏感度沿X=40 km的垂向切片黑色虛線指示異常體位置.Fig.4 Vertical slices of all stations Tx+Ty sensitivity along X=40 km at different periodsThe dashed black lines indicate the locations of the anomalies.

圖5為模型中心位置數(shù)據(jù)接收點(diǎn)(圖1中白色三角)位置的全分量敏感度，雖然該點(diǎn)遠(yuǎn)離異常體正上方，但同樣對異常體有敏感度.由于該點(diǎn)位于合成模型的中心對稱點(diǎn)，敏感度在淺部橫向?qū)ΨQ分布，圖6的深度5 km的橫向切片顯示敏感度呈現(xiàn)中心對稱的十字型，周期4.3 s處敏感度達(dá)到峰值，表明在5 km處該周期起主導(dǎo)作用.

圖5 不同周期下模型中心點(diǎn)位置全分量敏感度沿X=40 km的垂向切片黑色虛線指示異常體位置.Fig.5 Vertical slices of full component sensitivity of the model midpoint along X=40 km at different periodsThe dashed black lines indicate the locations of the anomalies.

圖6 不同周期下模型中心點(diǎn)位置全分量敏感度在深度5 km的水平切片黑色虛線指示異常體位置.Fig.6 Horizontal slices of full component sensitivity of the model midpoint along Z=5 km at different periodsThe dashed black lines indicate the locations of the anomalies.

3 二維與三維反演的對比

本文介紹的敏感度矩陣計(jì)算方法同樣適用于二維反演，通過計(jì)算二維大地電磁反演中的敏感度矩陣，定量對比理論二維、三維反演的敏感度.研究選用圖1中平行于YOZ平面的黑色實(shí)線剖面計(jì)算二維敏感度，周期范圍為0.2～8000 s，等對數(shù)間隔頻點(diǎn)25個.二維反演通常使用TM模式(蔡軍濤和陳小斌, 2010)，當(dāng)剖面平行于YOZ平面時，二維TM模式對應(yīng)三維YX模式，三維Zyx分量和二維TM分量具有一定的可比性.將剖面上所有網(wǎng)格求和并除以敏感度矩陣最大值歸一化，如圖7所示，雖然隨著周期增加，兩種反演的敏感度均有所下降，但二維反演敏感度始終高于三維反演，且在長周期部分差距更為明顯.需要指出的是，真實(shí)的地下結(jié)構(gòu)往往具有一定三維性，實(shí)測數(shù)據(jù)可能不滿足二維構(gòu)造假設(shè).為更好地指導(dǎo)實(shí)際應(yīng)用，有必要對比三維合成數(shù)據(jù)的二維與三維反演結(jié)果，并探究反演模型中異常體邊界的劃定方法.

圖7 二維與三維下反演歸一化敏感度對比Fig.7 Comparison of normalized sensitivity in two-dimensional and three-dimensional inversions

3.1 異常體邊界的劃定

二維與三維反演使用第2節(jié)中加入誤差的三維合成數(shù)據(jù)，其中傾子給定0.02的絕對誤差，阻抗張量主對角元素給定10%的相對誤差，非對角元素給定5%的相對誤差.三維反演使用全阻抗張量和傾子矢量，二維反演使用TM分量.反演的初始模型為100 Ωm均勻半空間，初始正則化因子100，二維和三維反演結(jié)果的均方根數(shù)據(jù)擬合誤差均收斂到1.05.本文選擇了兩條分別代表弱三維性(圖1黑實(shí)線)和強(qiáng)三維性(圖1白實(shí)線)剖面來進(jìn)行二維與三維反演結(jié)果的對比.

由于體積效應(yīng)的存在，大地電磁反演結(jié)果中，異常體與背景構(gòu)造之間存在電阻率值過渡帶.如圖8所示，統(tǒng)計(jì)三維反演中三個低阻異常體C1—C3及其過渡帶的電阻率分布，直方圖顯示電阻率分布存在峰值，在超過95%分位數(shù)時，電阻率值出現(xiàn)頻次迅速降低，表明該閾值可以視為異常體邊界.二維反演結(jié)果的邊界重建方法相同，選取95%分位數(shù)作為閾值.多次的合成數(shù)據(jù)測試顯示，95%分位數(shù)的閾值雖然是人為給定的數(shù)值，但具有一定的通用性.異常體中心測線的二維與三維反演結(jié)果如圖9所示，黑色虛線指示異常體的真實(shí)位置，黑色菱形為重建異常體內(nèi)部的極值，表示異常體的中心位置.三維反演和二維反演對淺部低阻異常體C1和C2的重建結(jié)果相似，三維反演得到的邊界與真實(shí)邊界更吻合.對于深部低阻異常體C3，三維反演敏感度在該深度范圍較小，難以恢復(fù)異常體電阻率值，但使用全阻抗張量進(jìn)行反演，重建的異常體邊界與真實(shí)邊界較為一致.雖然三維反演比二維反演更加吻合真實(shí)邊界，但三維反演得到的異常體中心位置偏淺，二維反演得到的異常體中心更接近真實(shí)位置，并且在此深度范圍內(nèi)依然有較大敏感度，能夠有效恢復(fù)電阻率值.圖10為遠(yuǎn)離異常體中心測線的二維和三維反演，結(jié)果顯示三維反演依然能夠重建異常體邊界，強(qiáng)三維性導(dǎo)致二維反演無法重建異常體，“拖尾”現(xiàn)象嚴(yán)重，出現(xiàn)虛假高阻異常.

圖8 三維反演中異常體(a) C1, (b) C2和(c) C3及其過渡帶的電阻率分布黑色虛線指示95%分位數(shù)的電阻率值.Fig.8 Electrical resistivity distribution of anomalies (a) C1, (b) C2 and (c) C3 and their transition zones in the three-dimensional inversionThe electrical resistivity values at the black dashed line are at 95% quantile.

圖9 合成數(shù)據(jù)(a)三維、(b)二維反演的電阻率模型垂直切片(X=40 km)黑色虛線指示異常體真實(shí)位置.黑色實(shí)線指示重建的異常體邊界.黑色三角為計(jì)算垂向分辨率矩陣的位置.黑色菱形指示重建的異常體中心.Fig.9 Vertical slices of synthetic data (a) three-dimensional inversion and (b) two-dimensional inversion electrical resistivity model along X=40 kmThe dashed black lines indicate the true locations of the anomalies. The solid black lines indicate the reconstruction boundary of anomalies. The black triangles show the location of vertical resolution matrix. The black diamonds indicate the reconstruction center of anomalies.

圖10 合成數(shù)據(jù)(a)三維、(b)二維反演的電阻率模型垂直切片(X=60 km)黑色虛線指示異常體真實(shí)位置.黑色實(shí)線指示重建的異常體邊界.黑色三角為計(jì)算垂向分辨率矩陣的位置.黑色菱形形指示重建的異常體中心.Fig.10 Vertical slices of synthetic data (a) three-dimensional inversion and (b) two-dimensional inversion electrical resistivity model along X=60 kmThe dashed black lines indicate the true locations of the anomalies. The solid black lines indicate the reconstruction boundary of anomalies. The black triangles show the location of vertical resolution matrix. The black diamonds indicate the reconstruction center of anomalies.

3.2 模型垂向分辨率

根據(jù)完整的敏感度矩陣，仿照線性反演中模型分辨率矩陣的定義式(9)(Egbert and Booker, 1992)計(jì)算了反演結(jié)果中不同位置的模型垂向分辨率其中P為二階正則化矩陣，λ為最終反演結(jié)果對應(yīng)的正則化因子.本文中垂向網(wǎng)格劃分為60個，頻點(diǎn)25個，阻抗張量和傾子矢量共6分量，考慮所有測點(diǎn)的影響總和，實(shí)際計(jì)算時只截取三維模型中的垂向一列，得到150×60的敏感度矩陣.理想狀態(tài)下的模型垂向分辨率矩陣應(yīng)該為一單位矩陣.如圖11a所示，三維反演在(X=40 km,Y=-40 km)處分辨率極值主要集中于對角線兩側(cè)，等值線隨深度增加逐漸擴(kuò)散，表明反演結(jié)果是真實(shí)模型均勻平滑后的結(jié)果.圖11b中二維反演的分辨率等值線在深度20～50 km向上偏離對角線，100 km后向下彎曲偏離對角線，說明二維反演結(jié)果在20～50 km深度上可能偏向于深部結(jié)構(gòu)，在100 km深度以下偏向于淺于100 km的真實(shí)結(jié)構(gòu)(即C3).圖11c和圖11d分別為三維、二維反演結(jié)果位于(X=60 km,Y=-40 km)處的垂向分辨率矩陣，由于該測點(diǎn)下方結(jié)構(gòu)的三維性更強(qiáng)，二維反演分辨率等值線偏離對角線程度更大.對比結(jié)果顯示當(dāng)結(jié)構(gòu)具有強(qiáng)三維性時，二維反演的垂向分辨率較差.

圖11 (a)三維反演和(b)二維反演位于(X=40 km, Y=-40 km)處的垂向分辨率矩陣.(c)三維反演和(d)二維反演位于(X=60 km, Y=-40 km)處的垂向分辨率矩陣黑色虛線指示分辨率矩陣的對角線.Fig.11 Vertical resolution matrix at the location of (X=40 km, Y=-40 km) for (a) three-dimensional inversion and (b) two-dimensional inversion. Vertical resolution matrix at the location of (X=60 km, Y=-40 km) for (c) three-dimensional inversion and (d) two-dimensional inversionThe dashed black line indicates the diagonal line of the resolution matrix.

R=P-1JT(JP-1JT+λI)-1J,

(9)

4 結(jié)論

本文采用求解大地電磁正演伴隨問題的方法，計(jì)算各個頻點(diǎn)、觀測分量和數(shù)據(jù)接收點(diǎn)在模型上的響應(yīng)，得到了完整敏感度矩陣.敏感度矩陣對網(wǎng)格大小、數(shù)據(jù)接收點(diǎn)位置和目標(biāo)函數(shù)的變化等均有響應(yīng)，因此分析完整敏感度矩陣，有助于優(yōu)化反演網(wǎng)格劃分、野外測點(diǎn)布設(shè)和反演參數(shù)選擇.本文在敏感度矩陣基礎(chǔ)上計(jì)算模型垂向分辨率，并結(jié)合反演結(jié)果的電阻率分布，提出了一種定量表征異常體邊界的方法.二維與三維反演合成模型測試表明，三維反演對淺部低阻異常體的重建優(yōu)于二維反演.當(dāng)電性結(jié)構(gòu)三維性不強(qiáng)時，二維反演在深部具有更高敏感度，可以有效恢復(fù)深部低阻異常體的電阻率值，為深部地質(zhì)解釋提供定量約束.三維反演雖然能有效重建異常體邊界，但無法恢復(fù)深部低阻異常體的電阻率值，并且重建的異常體中心偏淺.當(dāng)電性結(jié)構(gòu)三維性較強(qiáng)，三維反演重建的邊界與真實(shí)位置更加吻合.二維反演所出現(xiàn)的“拖尾”現(xiàn)象并非對應(yīng)真實(shí)結(jié)構(gòu)，虛假異常體會對精細(xì)地質(zhì)解釋產(chǎn)生阻礙.因此針對由多條長測線交叉組成的非陣列式數(shù)據(jù)，建議在三維性不強(qiáng)的研究區(qū)域，聯(lián)合使用二維和三維反演技術(shù)，提高深部地質(zhì)解釋的可靠性和準(zhǔn)確性.