何浩源， 李桐林， 張镕哲*， 朱威

1 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春 130012 2 中國地質(zhì)科學(xué)院地球物理地球化學(xué)勘查研究所自然資源部地球物理電磁法探測技術(shù)重點(diǎn)實(shí)驗(yàn)室， 廊坊 065000

0 引言

重力、重力梯度和大地電磁方法都是地球物理勘探領(lǐng)域中的重要技術(shù)手段，無論是在理論還是實(shí)際應(yīng)用方面均已得到國內(nèi)外學(xué)者的廣泛關(guān)注與研究.這些方法出于原理或?qū)嵺`上的差異，在分辨能力上表現(xiàn)出不同特點(diǎn).例如，重力和重力梯度數(shù)據(jù)雖具有較高的橫向分辨率，卻因信號隨深度增加的迅速衰減，使其縱向分辨率較低.重力梯度數(shù)據(jù)作為重力數(shù)據(jù)的一階導(dǎo)數(shù)，含有更多高頻信息，但低頻信息量少于重力數(shù)據(jù).因此，重力梯度數(shù)據(jù)相比重力數(shù)據(jù)對淺層區(qū)域的分辨能力更強(qiáng)，而對深層區(qū)域的分辨能力較弱.大地電磁數(shù)據(jù)則因包含大量頻率信息而更具測深能力，但通常大地電磁法的野外實(shí)測間距較大，測點(diǎn)較稀疏，使該方法的橫向分辨能力相對較弱.單一物探數(shù)據(jù)只能從單一物性角度對地下介質(zhì)作出評價(jià)，考慮到重力、重力梯度和大地電磁法各自的優(yōu)缺點(diǎn)，單一物探方法的解釋結(jié)果經(jīng)常有失準(zhǔn)確性.但如果采用多種數(shù)據(jù)從多個(gè)角度對同一區(qū)域的地下介質(zhì)進(jìn)行研究，則可以更加全面地預(yù)測地質(zhì)信息.因此，多種地球物理數(shù)據(jù)綜合解釋已成為地球物理勘探發(fā)展的重要趨勢.而聯(lián)合反演作為多種地球物理數(shù)據(jù)綜合解釋的有效方法之一，是降低反演多解性和改善單一方法局限性問題的重要手段.

20世紀(jì)70年代，Vozoff和Jupp(1975)首次提出地球物理聯(lián)合反演的思想，并將其應(yīng)用于直流電阻率法和大地電磁法的聯(lián)合解釋.隨后的45年中，聯(lián)合反演技術(shù)發(fā)展迅速，從一維聯(lián)合反演發(fā)展到三維聯(lián)合，從同一物性聯(lián)合反演發(fā)展到不同物性聯(lián)合，從物性關(guān)系聯(lián)合反演發(fā)展到結(jié)構(gòu)約束聯(lián)合.根據(jù)物性參數(shù)間構(gòu)建聯(lián)合關(guān)系方式的不同，可以將聯(lián)合反演分為兩大類.第一類為巖石物性關(guān)系聯(lián)合反演，該類方法需要提前構(gòu)建不同物性參數(shù)間的經(jīng)驗(yàn)關(guān)系式，通過經(jīng)驗(yàn)關(guān)系式進(jìn)行物性參數(shù)間的轉(zhuǎn)換(Nielsen and Jacobsen，2000；Afnimar et al.，2002；Colombo and De Stefano，2007).然而不同地區(qū)的物性關(guān)系存在差異，且復(fù)雜地區(qū)的經(jīng)驗(yàn)關(guān)系式難以建立，因此該方法具有一定局限性.第二類為結(jié)構(gòu)相似性聯(lián)合反演，該類方法假設(shè)同一區(qū)域內(nèi)，不同物性參數(shù)具有相同或部分相同的地下結(jié)構(gòu).目前，由Gallardo和Meju(2003)提出的基于交叉梯度約束的聯(lián)合反演已成為最受歡迎的結(jié)構(gòu)相似性聯(lián)合反演方法之一.該方法通過約束不同物性參數(shù)的梯度方向一致，來尋找結(jié)構(gòu)相似的模型，既不依賴經(jīng)驗(yàn)關(guān)系式，也不需要對初始模型做任何假設(shè)，是一種有效的地球物理綜合解釋手段(Gallardo and Meju，2004；彭淼等，2013；Zhang et al.，2019).

重力勘探和大地電磁勘探因具有覆蓋面廣、效率高和費(fèi)用低等特點(diǎn)而獲得廣泛應(yīng)用，兩種方法間的交叉梯度聯(lián)合反演也受到了地球物理學(xué)者們的關(guān)注和研究.其中，李桐林等(2016)實(shí)現(xiàn)了部分區(qū)域約束的重力、磁法和大地電磁的二維交叉梯度聯(lián)合反演；Pak等(2017)開展了二維數(shù)據(jù)空間下的重力、磁法和大地電磁的交叉梯度聯(lián)合反演；張镕哲等(2019)將交叉梯度約束應(yīng)用于重力、磁法、大地電磁和地震初至波走時(shí)四種方法的二維聯(lián)合反演中；閆政文等(2020)實(shí)現(xiàn)了三維重力、磁法和大地電磁法的交叉梯度聯(lián)合反演.上述研究成果均表明，基于交叉梯度約束的聯(lián)合反演對異常的分辨能力相比單獨(dú)反演更強(qiáng)，但也存在一定的局限性:1)現(xiàn)有的重力和大地電磁聯(lián)合反演更多集中在二維研究，三維研究較少.2)由于重力數(shù)據(jù)縱向分辨能力差，通過交叉梯度結(jié)構(gòu)約束傳遞的不準(zhǔn)確縱向密度信息通常會誤導(dǎo)電阻率模型，導(dǎo)致大地電磁反演結(jié)果偏離真實(shí)模型.3)平滑正則化方法在聯(lián)合反演中的應(yīng)用會使結(jié)果表現(xiàn)出異常發(fā)散的現(xiàn)象，影響聯(lián)合反演的分辨率.4)傳統(tǒng)的連續(xù)式聯(lián)合反演策略，常因交叉梯度的結(jié)構(gòu)相似約束使反演模型向真實(shí)模型迭代的趨勢受到限制，導(dǎo)致反演模型難以收斂.

為解決上述問題，本文提出了一種基于交叉梯度約束的三維重力、重力梯度和大地電磁數(shù)據(jù)聯(lián)合反演方法.首先，引入重力梯度數(shù)據(jù)能夠?yàn)榻徊嫣荻嚷?lián)合提供更為可靠的密度模型，有利于提高聯(lián)合反演的準(zhǔn)確性.其次，為了提高電阻率模型的自身分辨率，本文提出了一種平滑聚焦的正則化方法.該方法既能解決平滑反演結(jié)果異常體發(fā)散的問題，又能避免聚焦反演伴生的近地表假異常問題.同時(shí)，為聚焦和平滑聚焦正則化加入了聚焦強(qiáng)度閾值方法，以限制高物性異常區(qū)域的聚焦優(yōu)勢，避免過度聚焦問題.此外，借鑒Um等(2014)交替進(jìn)行結(jié)構(gòu)約束和數(shù)據(jù)約束的思想，本文提出了一套間斷式聯(lián)合的交叉梯度聯(lián)合反演流程.整個(gè)反演流程僅經(jīng)歷一次收斂過程，且避免了某一方法對聯(lián)合反演結(jié)果絕對主導(dǎo)的問題.在保證約束不同物性模型間結(jié)構(gòu)相似的同時(shí)，又能一定程度上減小結(jié)構(gòu)相似性約束對反演模型向真實(shí)模型迭代帶來的限制.最后，通過理論模型測試，驗(yàn)證了本文所提出方法的有效性和準(zhǔn)確性.

1 三維重力、重力梯度和大地電磁正反演理論

本文中，重力、重力梯度及大地電磁的正反演均在笛卡爾坐標(biāo)系下進(jìn)行計(jì)算，地下模型被劃分為若干塊物性均勻的規(guī)則長方體.

1.1 三維重力、重力梯度和大地電磁正演

根據(jù)Nagy(1966)，地下長方體在測點(diǎn)處產(chǎn)生的重力Vz及重力梯度Vzz的正演計(jì)算公式如下：

(1)

(2)

圖1 長方體與測點(diǎn)坐標(biāo)關(guān)系示意圖Fig.1 Coordinate relationship between cuboid and measuring point

其中，(ξ1，ξ2)、(η1，η2)和(ζ1，ζ2)依次為長方體沿x、y、z軸的坐標(biāo)范圍，(x,y,z)為觀測點(diǎn)P的坐標(biāo)，空間關(guān)系如圖1所示.G為萬有引力常數(shù)，ρ為長方體密度，r為長方體中心與觀測點(diǎn)P的距離.

假設(shè)地下模型被劃分成M=Nx×Ny×Nz塊小長方體，則有N=Nx×Ny個(gè)測點(diǎn)位于各長方體中心沿Z方向在水平地表的投影處.重力或重力梯度的觀測數(shù)據(jù)d與模型物性m的正演關(guān)系可寫成如下矩陣形式：

d=Jm,

(3)

式中，J是N×M的靈敏度矩陣.

大地電磁正演及計(jì)算靈敏度矩陣J的部分主要參考三維大地電磁反演程序WSINV3D(Siripunvaraporn et al.，2005).該程序采用FDE(Finite Difference Equations)方案(Yee，1966)進(jìn)行有限差分交錯(cuò)網(wǎng)格劃分，即電場分量定于網(wǎng)格棱邊的中點(diǎn)，方向沿著棱邊，磁場分量定于網(wǎng)格表面的中心，方向沿著面法向分量.該方案由麥克斯韋方程組出發(fā)，消去磁場，得到電場的波動方程.對波動方程離散化后使用QMR算法求解電場，再由電場計(jì)算磁場，最終計(jì)算出觀測點(diǎn)處的四個(gè)阻抗分量Zxx、Zxy、Zyx和Zyy.

1.2 三維重力、重力梯度和大地電磁單獨(dú)反演

本文中，重力、重力梯度及大地電磁均采用了模型空間下的正則化反演，其目標(biāo)函數(shù)如下：

(4)

對式(4)求導(dǎo)，令求導(dǎo)結(jié)果等于0，可以得到模型空間下的反演模型迭代表達(dá)式：

(5)

在沒有明確參考模型的情況下，每一次迭代時(shí)將當(dāng)前模型mk設(shè)置為參考模型mref.使用共軛梯度法求解線性方程組即可得到每次迭代的模型更新量.

2 三維重力、重力梯度和大地電磁數(shù)據(jù)交叉梯度約束聯(lián)合反演理論

重力、重力梯度和大地電磁數(shù)據(jù)聯(lián)合反演的主要聯(lián)合關(guān)系為：重力與重力梯度直接聯(lián)立方程組，聯(lián)合求解剩余密度模型(本文中為了方便理清關(guān)系，將此聯(lián)合稱為密度單獨(dú)反演)；剩余密度模型與大地電磁對應(yīng)的電阻率模型進(jìn)行交叉梯度約束聯(lián)合.聯(lián)合關(guān)系框架如圖2所示.

圖2 重力、重力梯度和大地電磁聯(lián)合關(guān)系示意圖Fig.2 Joint relation diagram of gravity， gravity gradient and magnetotelluric

2.1 重力和重力梯度數(shù)據(jù)直接聯(lián)合反演

重力和重力梯度的聯(lián)合反演(即密度單獨(dú)反演)，為密度模型不同觀測數(shù)據(jù)之間的聯(lián)合反演，直接聯(lián)立它們的方程組求解即可.綜合考慮反演的計(jì)算量和結(jié)果的準(zhǔn)確性，本文只采用Vz和Vzz分量進(jìn)行聯(lián)合反演(秦朋波，2016).令：

(6)

(7)

(8)

將式(6)—(8)代入式(5)中，得到重力和重力梯度聯(lián)合反演的模型迭代表達(dá)式：

×[dgra-Jgramk],

(9)

A·Δm=b，

(10)

(11)

求解出WmΔm，進(jìn)而求得模型更新量.

2.2 交叉梯度約束與聯(lián)合網(wǎng)格剖分

根據(jù)Gallardo和Meju(2003)定義的三維交叉梯度函數(shù)：

(12)

其中，mρ為剩余密度，mσ為對數(shù)電阻率.交叉梯度的三個(gè)分量分別為

(13)

(14)

(15)

采用向前差分的方式對式(13)—(15)進(jìn)行離散(Fregoso and Gallardo，2009)：

-mρ d(mσ f-mσ c)],

(16)

-mρ r(mσ d-mσ c)],

(17)

-mρ f(mσ r-mσ c)],

(18)

其中，mc為當(dāng)前塊物性，mr為X方向下一塊的物性，mf為Y方向下一塊的物性，md為Z方向下一塊的物性.

將交叉梯度項(xiàng)加入式(4)，得到交叉梯度聯(lián)合反演的目標(biāo)函數(shù)：

+βρtT(mρ,mσ)t(mρ,mσ),

(19)

+βσtT(mρ,mσ)t(mρ,mσ),

(20)

其中，β為交叉梯度約束項(xiàng)權(quán)重系數(shù).令式(19)和式(20)的導(dǎo)數(shù)等于0，得到模型迭代表達(dá)式：

(21)

(22)

式(21)和式(22)的等號右端項(xiàng)對應(yīng)第k次迭代的數(shù)據(jù).B為交叉梯度的離散導(dǎo)數(shù)，其求取方式可參考Fregoso和Gallardo(2009)的文章，本文限于篇幅不作展開.

圖3 聯(lián)合網(wǎng)格劃分示意圖Fig.3 Diagram of joint grid generation

大地電磁反演通常要求縱向網(wǎng)格間距隨著深度增加而逐漸增加，而重力和重力梯度對網(wǎng)格劃分的要求則較為寬松.因此，如圖3所示，聯(lián)合網(wǎng)格的公共部分(中心區(qū))采用X、Y方向均勻剖分，Z方向網(wǎng)格距隨深度增加而增加的剖分方式.為滿足大地電磁邊界條件，還需要在中心區(qū)外設(shè)置擴(kuò)邊區(qū).擴(kuò)邊區(qū)內(nèi)，越遠(yuǎn)離中心區(qū)的網(wǎng)格，X、Y、Z方向的間距越大.

考慮到計(jì)算交叉梯度要求兩種物性模型網(wǎng)格一致，在為大地電磁聯(lián)合反演計(jì)算交叉梯度之前，需要先將剩余密度模型擴(kuò)邊，擴(kuò)邊區(qū)剩余密度按照背景值設(shè)定.而在為重力和重力梯度聯(lián)合反演計(jì)算交叉梯度之前，則需要先將電阻率模型縮邊，只保留中心區(qū).

2.3 平滑聚焦的模型正則化方法

平滑約束和聚焦約束的有效性無論在理論還是應(yīng)用上均已得到充分證明(Constable et al.，1987；Zhdanov，2009；劉小軍，2007).本文結(jié)合兩種傳統(tǒng)約束方法，為大地電磁反演提出了一種平滑聚焦的模型正則化方法，并將其應(yīng)用于聯(lián)合反演中.該正則化方法既能避免平滑反演存在的異常體發(fā)散問題，又能解決聚焦反演伴生的近地表假異常問題，具體方法如下.

平滑聚焦正則化的平滑部分如下：

(23)

(24)

(25)

(26)

其中，Rx、Ry和Rz分別用于約束模型X、Y、Z方向的一階導(dǎo)數(shù).對于Rx，若當(dāng)前塊非X方向的最后一塊，則用X方向的下一塊減當(dāng)前塊作為一階導(dǎo)數(shù)；若當(dāng)前塊是X方向的最后一塊，則用當(dāng)前塊減X方向的上一塊作為一階導(dǎo)數(shù).Ry和Rz同理.

平滑聚焦正則化的聚焦部分如下(Zhdanov，2015；孫思源，2019)：

(27)

Re=diag(1/me),

(28)

其中，常數(shù)e為聚焦因子，e越小，迭代模型越易聚焦.在沒有明確參考模型的情況下，通常令此處的mref等于0，即均勻半空間.me是M×1的矩陣，Re為模型聚焦項(xiàng)，是以1/me為主對角線的M×M矩陣.

但該聚焦方式經(jīng)常存在過度聚焦的問題.由式(27)和式(28)，當(dāng)前模型物性參數(shù)與均勻半空間相差越大的塊，越容易聚焦，越容易在下次迭代中產(chǎn)生更大的物性差異.這種聚焦優(yōu)勢經(jīng)常會隨著模型迭代逐漸增大，最終導(dǎo)致反演結(jié)果的物性異常聚集在比真實(shí)模型邊界更小的范圍內(nèi).僅根據(jù)實(shí)際資料在迭代過程中對物性參數(shù)值的上下限進(jìn)行約束尚不足以充分解決此問題.因此，本文對模型聚焦項(xiàng)進(jìn)行了改進(jìn)，通過設(shè)置閾值et對聚焦強(qiáng)度進(jìn)行限制.由式(27)得到me后，將me中所有高于et的值降至et，再將新得到的me代入后續(xù)計(jì)算.根據(jù)模型測試經(jīng)驗(yàn)，et通?？扇≡璵e最大值的1/2至2/3.

模型正則化的本質(zhì)是約束各模型塊物性參數(shù)間的數(shù)學(xué)關(guān)系，平滑約束令相鄰模型塊間的物性參數(shù)變化極小，而聚焦約束對物性參數(shù)與均勻半空間相差大的塊約束小，對相差小的塊約束大.同時(shí)利用兩種約束關(guān)系，豐富模型限制條件，有利于減少反演的多解性.于是為平滑聚焦正則化方法設(shè)計(jì)如下的模型協(xié)方差矩陣：

(29)

在迭代初期，由于物性起伏較小，Re的聚焦作用并不明顯，平滑聚焦約束整體偏向平滑.隨著迭代繼續(xù)，物性起伏愈發(fā)明顯，Re的聚焦作用也愈發(fā)突出，平滑聚焦約束逐漸偏向聚焦.考慮到如果任一物性在某處的梯度值為0，則該處兩種物性的交叉梯度也為0.那么在交叉梯度聯(lián)合反演的迭代初期，通過偏向平滑的正則化使物性異常發(fā)散，能夠得到更大范圍的有效交叉梯度約束，有利于提高兩種物性模型結(jié)構(gòu)變化的同步性.

(30)

(31)

其中，Rs為基于靈敏度矩陣的深度加權(quán)項(xiàng)(Mehanee et al.，1998)，是一個(gè)M×M的對角陣.

2.4 重力、重力梯度和大地電磁數(shù)據(jù)聯(lián)合反演流程

圖4 重力、重力梯度和大地電磁聯(lián)合反演流程圖Fig.4 Flow chart of gravity，gravity gradient and magnetotelluric joint inversion

如圖4所示，在交叉梯度聯(lián)合之前，先對兩種物性的初始模型進(jìn)行幾次單獨(dú)反演迭代，使物性模型初顯結(jié)構(gòu).考慮到交叉梯度聯(lián)合反演具有約束兩種物性結(jié)構(gòu)相似的特性，而迭代過程中常常存在兩種物性結(jié)構(gòu)已近似一致，但均與真實(shí)模型存在偏差的情況.在傳統(tǒng)的連續(xù)式聯(lián)合反演中，由于兩種物性間的結(jié)構(gòu)約束是分步、交替且連續(xù)進(jìn)行的，反而使反演模型向真實(shí)模型迭代的趨勢受到限制，導(dǎo)致聯(lián)合反演收斂困難，反演結(jié)果發(fā)散.

Um等(2014)為了克服地震和電磁聯(lián)合反演收斂困難的問題，提出了一種交替進(jìn)行結(jié)構(gòu)約束和數(shù)據(jù)約束的反演策略.其反演流程中，先分別令地震和電磁數(shù)據(jù)單獨(dú)反演至收斂，然后固定所得速度模型對所得電阻率模型進(jìn)行僅含模型約束項(xiàng)和交叉梯度約束項(xiàng)的純結(jié)構(gòu)約束反演，再將結(jié)構(gòu)約束后的電阻率結(jié)果作為新的初始模型進(jìn)行電磁單獨(dú)反演直至收斂.一次純結(jié)構(gòu)約束反演與一次電磁單獨(dú)反演被合稱為一次電磁精化過程(EM refinement process).重復(fù)數(shù)次電磁精化過程直至連續(xù)兩次精化的電阻率結(jié)果基本相近，再固定電阻率模型對地震施行類似的精化過程，如此往復(fù)，最終得到聯(lián)合反演結(jié)果.該策略雖有效改善了聯(lián)合反演的收斂性，但也明顯存在結(jié)構(gòu)約束過強(qiáng)的問題.尤其是在純結(jié)構(gòu)約束反演中缺少數(shù)據(jù)項(xiàng)約束，而多次電磁精化后才反輪到電阻率結(jié)構(gòu)約束地震反演，導(dǎo)致聯(lián)合反演的結(jié)果幾乎完全取決于地震反演所得結(jié)構(gòu).此外，由于每次精化都意味著一個(gè)完整的反演收斂過程，反復(fù)地精化使聯(lián)合反演流程過于冗長.

借鑒Um等交替進(jìn)行結(jié)構(gòu)約束和數(shù)據(jù)約束的思路，本文設(shè)計(jì)了一套間斷式聯(lián)合的反演流程，即根據(jù)當(dāng)前迭代次數(shù)間隔式地令交叉梯度約束項(xiàng)權(quán)重系數(shù)β=0，相當(dāng)于讓聯(lián)合反演和單獨(dú)反演交替進(jìn)行.根據(jù)參與交叉梯度聯(lián)合的物性方法在可信度上的差異，又可為流程分出兩種間斷方案：若兩種物性方法可信度相近，則采用如圖5a的方案，在奇數(shù)次迭代中對方法一進(jìn)行交叉梯度約束，令方法二單獨(dú)反演，在偶數(shù)次迭代中對方法二進(jìn)行交叉梯度約束，令方法一單獨(dú)反演，整體上兩種方法的物性結(jié)構(gòu)互相影響；若方法二的可信度明顯高于方法一，則采用如圖5b的方案，在奇數(shù)次迭代中對兩種方法都進(jìn)行交叉梯度約束，在偶數(shù)次迭代中都不約束，整體上偏向于以方法二的物性結(jié)構(gòu)為主導(dǎo).對于第二種方案，雖然迭代過程中總是先由方法二的物性結(jié)構(gòu)約束方法一，但由于約束策略的變動精細(xì)到每次迭代，被約束的物性結(jié)構(gòu)不會立即與約束方高度一致，使得方法一的物性結(jié)構(gòu)仍對方法二的物性結(jié)構(gòu)有一定的約束作用，避免了某一方法因結(jié)構(gòu)約束過強(qiáng)而對聯(lián)合反演結(jié)果絕對主導(dǎo)的問題.具體到本文的重力、重力梯度和大地電磁聯(lián)合反演，可通過單獨(dú)反演測試來選擇合適的間斷方案.若大地電磁單獨(dú)反演結(jié)果呈現(xiàn)出明顯的縱向電阻率結(jié)構(gòu)，則選擇偏向以大地電磁為主導(dǎo)的方案.若物性無明顯縱向結(jié)構(gòu)，則可以選擇互相約束的方案，抑或?yàn)橥怀鰴M向物性變化而選擇偏向以重力和重力梯度為主導(dǎo)的方案.這種間斷式聯(lián)合的交叉梯度聯(lián)合反演策略在整個(gè)反演流程中僅經(jīng)歷一次收斂過程，相比于現(xiàn)有的連續(xù)式聯(lián)合，既能保證不同物性參數(shù)間的結(jié)構(gòu)相似性約束，又能讓迭代模型的響應(yīng)向觀測數(shù)據(jù)充分?jǐn)M合.

2.5 聯(lián)合反演參數(shù)取值方案

本節(jié)對模型約束項(xiàng)權(quán)重系數(shù)λ、交叉梯度約束項(xiàng)權(quán)重系數(shù)β和交叉梯度聯(lián)合前單獨(dú)反演次數(shù)的取值方案進(jìn)行說明.λ和β的計(jì)算方法如下：

圖5 間斷式聯(lián)合方法示意圖 (a) 方案一； (b) 方案二.Fig.5 Diagram of discontinuous joint method (a) Scheme 1; (b) Scheme 2.

(32)

(33)

通過比值方法為系數(shù)λ和β賦值，能夠保證數(shù)據(jù)約束項(xiàng)、模型約束項(xiàng)和交叉梯度約束項(xiàng)在數(shù)量級上相近.對于重力和重力梯度反演，由式(32)為λ賦初值，并在迭代過程中采用與張镕哲(2020)相同的阻尼因子方法對λ進(jìn)行修正.而大地電磁反演的λ及所有反演的β則均采用隨迭代次數(shù)線性減小λ0和β0的方法進(jìn)行取值：

(34)

(35)

其中，k為當(dāng)前迭代次數(shù)，Nit為最大迭代次數(shù)，a和b為常數(shù).

3 模型試驗(yàn)

采用雙長方體理論模型對提出的算法進(jìn)行模型試驗(yàn).如圖6a和圖6b，左側(cè)長方體大小為2000 m×1500 m×2000 m，頂面埋深850 m，剩余密度為1 g·cm-3，電阻率為1000 Ωm；右側(cè)長方體大小為2000 m×1500 m×1900 m，頂面埋深1500 m，剩余密度為1 g·cm-3，電阻率為10 Ωm；背景剩余密度為0，電阻率為100 Ωm.中心區(qū)采用22×22×14的三維剖分，X、Y方向?yàn)榫鶆蚱史?，網(wǎng)格間距為500 m.Z方向網(wǎng)格間距隨著深度增加逐漸從200 m增加到2000 m.中心區(qū)外，在X、Y的正負(fù)方向及Z的正方向各擴(kuò)邊4個(gè)網(wǎng)格，間距依次為4000、8000、1600和32000 m，地下網(wǎng)格共30×30×18塊.重力和重力梯度的測點(diǎn)數(shù)均為22×22個(gè)，分布于中心區(qū)地表各網(wǎng)格的中心處.大地電磁共8×8個(gè)測點(diǎn)，如圖6c所示布設(shè)于地表.大地電磁采用9個(gè)頻率觀測，依次為100、10、5、3 Hz、2、1.5、1、0.5 Hz和0.25 Hz.為重力、重力梯度和大地電磁的合成數(shù)據(jù)分別加入3%隨機(jī)噪聲.

圖6 雙長方體模型設(shè)計(jì)示意圖 (a) 剩余密度模型三維透視圖； (b) 對數(shù)電阻率模型三維透視圖； (c) 大地電磁測點(diǎn)分布圖.Fig.6 Design diagram of double cuboid model (a) 3D perspective of residual density model； (b) 3D perspective of log resistivity model； (c) Distribution of magnetotelluric measuring points.

各組模型試驗(yàn)均使用均勻半空間作為初始模型，并保證充分迭代.使用數(shù)據(jù)擬合差RMS和模型擬合差RMSE對反演結(jié)果進(jìn)行評價(jià)：

(36)

(37)

其中，dobs為觀測數(shù)據(jù)，dinv為反演結(jié)果的正演響應(yīng)，mtrue為真實(shí)模型，minv為反演結(jié)果模型，err為觀測誤差.RMS越小表示數(shù)據(jù)擬合程度更高，RMSE越小則表示模型還原程度更高.

模型試驗(yàn)3.1—3.3節(jié)中的所有二維模型切片均在x=5500 m處，為方便對比，使用黑色粗線框在反演結(jié)果切片中標(biāo)記了異常體的真實(shí)邊界位置.試驗(yàn)中提及的所有“密度”均指剩余密度，“電阻率”均指對數(shù)電阻率.

3.1 正則化方法對比試算

3.2 不同數(shù)據(jù)的交叉梯度聯(lián)合反演對比試算

圖7 (a) 真實(shí)電阻率模型切片； (b) 大地電磁平滑反演結(jié)果切片； (c) 大地電磁聚焦反演結(jié)果切片； (d) 大地電磁平滑聚焦反演結(jié)果切片F(xiàn)ig.7 (a) Real resistivity slice； (b) Slice of magnetotelluric smooth inversion； (c) Slice of magnetotelluric focusing inversion； (d) Slice of magnetotelluric smooth-focusing inversion

圖8 不同正則化方法大地電磁單獨(dú)反演RMS曲線Fig.8 RMS curve of magnetotelluric inversion by different regularization methods

表1 不同正則化方法大地電磁單獨(dú)反演結(jié)果的 RMS和RMSE值Table 1 RMS and RMSE of magnetotelluric single inversion results by different regularization methods

圖9a和圖9b分別為真實(shí)密度模型切片和真實(shí)電阻率模型切片；圖9c和圖9d分別為重力單獨(dú)反演結(jié)果切片和大地電磁單獨(dú)反演結(jié)果切片；圖9e和圖9f分別為重力和大地電磁交叉梯度聯(lián)合反演結(jié)果的密度切片和電阻率切片；圖9g和圖9h分別為重力、重力梯度和大地電磁交叉梯度聯(lián)合反演結(jié)果的密度切片和電阻率切片.圖10為上述反演迭代的RMS曲線，各組反演結(jié)果的RMS和RMSE如表2所示.由圖9(a—f)發(fā)現(xiàn)，重力縱向分辨能力弱，其單獨(dú)反演結(jié)果的縱向密度結(jié)構(gòu)與真實(shí)模型相差較大，無法為交叉梯度聯(lián)合反演提供可靠的密度模型.由于交叉梯度聯(lián)合方法的結(jié)構(gòu)相似性約束，不可靠的縱向密度模型影響了本來更為可靠的縱向電阻率模型，導(dǎo)致聯(lián)合反演效果差.而如圖9g和圖9h，異常體的位置、幾何形態(tài)和物性參數(shù)幅值明顯更還原真實(shí)模型，可見重力梯度信息的引入極大地提高了密度模型的可靠性，從而改善了聯(lián)合反演效果.圖11對比了重力、重力梯度和大地電磁聯(lián)合反演中，剩余密度和對數(shù)電阻率間交叉梯度的模的變化.圖11a為單獨(dú)反演階段結(jié)束時(shí)的交叉梯度結(jié)果切片，圖11b則為聯(lián)合反演最終結(jié)果的交叉梯度切片.交叉梯度模的減小說明了交叉梯度聯(lián)合能夠使兩種物性反演結(jié)果具有更高的結(jié)構(gòu)一致性.

表2 不同數(shù)據(jù)交叉梯度聯(lián)合反演結(jié)果的RMS和RMSE值

3.3 聚焦強(qiáng)度閾值和間斷式聯(lián)合算法試算

沿用上一節(jié)的反演參數(shù)，基于重力、重力梯度和大地電磁交叉梯度聯(lián)合反演對聚焦強(qiáng)度閾值和間斷式聯(lián)合算法進(jìn)行模型試算分析.圖12a和圖12b為無聚焦強(qiáng)度閾值的交叉梯度聯(lián)合反演結(jié)果的密度模型切片和電阻率模型切片；圖12c和圖12d為連續(xù)式聯(lián)合的交叉梯度聯(lián)合反演結(jié)果的密度模型切片和電阻率模型切片；圖12e和圖12f(同圖9g和圖9h)為加入聚焦強(qiáng)度閾值且間斷式聯(lián)合的交叉梯度聯(lián)合反演結(jié)果的密度切片和電阻率切片，其中聚焦強(qiáng)度閾值設(shè)置為原最大值的0.6倍.各組反演迭代的RMS曲線如圖13所示，反演結(jié)果的RMS和RMSE如表3所示.

對比圖12a和圖12e，發(fā)現(xiàn)圖12a的反演結(jié)果存在過度聚焦問題，高密度塊聚集在比真實(shí)模型邊界更小的范圍內(nèi)，而加入了聚焦強(qiáng)度閾值的圖12e物性參數(shù)聚焦地更為均勻，更吻合真實(shí)模型.對比圖12b和圖12f也能得到同樣的結(jié)論，即聚焦強(qiáng)度閾值能有效改善聚焦反演和平滑聚焦反演中存在的過度聚焦問題.如圖12c和圖12d，連續(xù)式聯(lián)合的每次迭代都進(jìn)行交叉梯度約束，當(dāng)密度和電阻率的結(jié)構(gòu)已近似一致時(shí)，兩種交替迭代的物性模型反而受結(jié)構(gòu)相似性約束的限制，難以進(jìn)一步向真實(shí)模型收斂.而間斷式聯(lián)合的方法，如圖12e和圖12f，明顯改善了這一問題.

圖9 (a) 真實(shí)密度模型切片； (b) 真實(shí)電阻率模型切片； (c) 重力單獨(dú)反演結(jié)果切片； (d) 大地電磁單獨(dú)反演結(jié)果切片； (e) 重力和大地電磁聯(lián)合反演結(jié)果密度切片； (f) 重力和大地電磁聯(lián)合反演結(jié)果電阻率切片； (g) 重力、重力梯度和 大地電磁聯(lián)合反演結(jié)果密度切片； (h) 重力、重力梯度和大地電磁聯(lián)合反演結(jié)果電阻率切片F(xiàn)ig.9 (a) Real density slice； (b) Real resistivity slice； (c) Slice of gravity single inversion； (d) Slice of magnetotelluric single inversion； (e) Density slice of gravity and magnetotelluric joint inversion； (f) Resistivity slice of gravity and magnetotelluric joint inversion； (g) Density slice of gravity， gravity gradient and magnetotelluric joint inversion； (h) Resistivity slice of gravity， gravity gradient and magnetotelluric joint inversion

3.4 復(fù)雜理論模型試算

為進(jìn)一步驗(yàn)證本文聯(lián)合反演算法的有效性，本節(jié)對復(fù)雜理論模型進(jìn)行試算.沿用與前一模型相同的網(wǎng)格劃分、測點(diǎn)布設(shè)及大地電磁觀測頻率組，同樣為重力、重力梯度和大地電磁的合成數(shù)據(jù)分別加入3%隨機(jī)噪聲.如圖14所示，模型由深淺不一的兩塊低阻高密度長方體和一塊高阻高密度長方體組成，另有一條高阻高密度臺階斜搭在其中一塊長方體上.

圖10 不同數(shù)據(jù)交叉梯度聯(lián)合反演RMS曲線 (a) 重力數(shù)據(jù)RMS； (b) 重力梯度數(shù)據(jù)RMS； (c) 大地電磁數(shù)據(jù)RMS.Fig.10 RMS curve of cross-gradient joint inversion with different data (a) RMS of gravity data; (b) RMS of gravity gradient data; (c) RMS of magnetotelluric data.

圖11 重力、重力梯度和大地電磁聯(lián)合反演交叉梯度變化對比圖 (a) 單獨(dú)反演階段結(jié)果交叉梯度模切片； (b) 聯(lián)合反演最終結(jié)果交叉梯度模切片.Fig.11 Contrast of cross-gradient for gravity，gravity gradient and magnetotelluric joint inversion (a) Modulus slice of cross-gradient for results in single inversion phase； (b) Modulus slice of cross-gradient for joint inversion final result.

圖12 (a) 無聚焦強(qiáng)度閾值的聯(lián)合反演結(jié)果密度切片； (b) 無聚焦強(qiáng)度閾值的聯(lián)合反演結(jié)果電阻率切片； (c) 連續(xù)式聯(lián)合反演結(jié)果密度切片； (d) 連續(xù)式聯(lián)合反演結(jié)果電阻率切片； (e) 加入聚焦強(qiáng)度閾值且間斷式聯(lián)合的反演結(jié)果密度切片； (f) 加入聚焦強(qiáng)度閾值且間斷式聯(lián)合的反演結(jié)果電阻率切片F(xiàn)ig.12 (a) Density slice of joint inversion without focusing intensity threshold； (b) Resistivity slice of joint inversion without focusing intensity threshold； (c) Density slice of continuous joint inversion； (d) Resistivity slice of continuous joint inversion； (e) Density slice of discontinuous joint inversion with focusing intensity threshold； (f) Resistivity slice of discontinuous joint inversion with focusing intensity threshold

圖13 不同方式交叉梯度聯(lián)合反演RMS曲線 (a) 重力數(shù)據(jù)RMS； (b) 重力梯度數(shù)據(jù)RMS； (c) 大地電磁數(shù)據(jù)RMS.Fig.13 RMS curve of cross-gradient joint inversion by different methods (a) RMS of gravity data; (b) RMS of gravity gradient data; (c) RMS of magnetotelluric data.

圖14 復(fù)雜模型設(shè)計(jì)示意圖 (a) 剩余密度模型三維透視圖； (b) 對數(shù)電阻率模型三維透視圖.Fig.14 Design diagram of complex model (a) 3D perspective of residual density model； (b) 3D perspective of log resistivity model.

圖15 復(fù)雜模型的重力、重力梯度和大地電磁交叉梯度聯(lián)合反演RMS曲線 (a) 重力數(shù)據(jù)RMS； (b) 重力梯度數(shù)據(jù)RMS； (c) 大地電磁數(shù)據(jù)RMS.Fig.15 RMS curve of gravity, gravity gradient and magnetotelluric cross-gradient joint inversion for complex model (a) RMS of gravity data; (b) RMS of gravity gradient data; (c) RMS of magnetotelluric data.

圖16 (a) 真實(shí)密度模型切片1(x=3000 m)； (b) 真實(shí)密度模型切片2(y=5500 m)； (c) 真實(shí)電阻率模型切片1(x=3000 m)； (d) 真實(shí)電阻率模型切片2(y=5500 m)； (e) 聯(lián)合反演結(jié)果密度切片1(x=3000 m)； (f) 聯(lián)合反演結(jié)果密度 切片2(y=5500 m)； (g) 聯(lián)合反演結(jié)果電阻率切片1(x=3000 m)； (h) 聯(lián)合反演結(jié)果電阻率切片2(y=5500 m)Fig.16 (a) Real density slice 1 (x=3000 m)； (b) Real density slice 2 (y=5500 m)； (c) Real resistivity slice 1 (x=3000 m)； (d) Real resistivity slice 2 (y=5500 m)； (e) Density slice 1 of joint inversion (x=3000 m)； (f) Density slice 2 of joint inversion (y=5500 m)； (g) Resistivity slice 1 of joint inversion (x=3000 m)； (h) Resistivity slice 2 of joint inversion (y=5500 m)

表3 不同方式交叉梯度聯(lián)合反演結(jié)果的RMS和RMSE值Table 3 RMS and RMSE of cross-gradient joint inversion results by different methods

4 結(jié)論

本文提出了一種基于平滑聚焦正則化法的三維重力、重力梯度和大地電磁數(shù)據(jù)交叉梯度聯(lián)合反演算法，并通過理論模型試算對比分析得出以下結(jié)論：

(1)將交叉梯度結(jié)構(gòu)相似約束應(yīng)用于聯(lián)合反演的重要前提是，參與聯(lián)合的兩種方法單獨(dú)反演均能得到具有一定可靠性的物性模型.由于重力反演縱向分辨能力弱，其不可靠的縱向密度模型會因?yàn)榻Y(jié)構(gòu)相似性約束影響本來更可靠的縱向電阻率模型，導(dǎo)致聯(lián)合反演結(jié)果準(zhǔn)確性降低.本文通過引入重力梯度數(shù)據(jù)有效地提高了密度模型的可靠性，從而提高了聯(lián)合反演的分辨率.

(2)本文提出的平滑聚焦正則化方法可以為大地電磁反演恢復(fù)出更真實(shí)的異常體幾何形態(tài)和物性參數(shù)幅值，既能避免平滑反演異常體發(fā)散、邊界模糊的問題，又能有效解決聚焦反演伴生的近地表假異常問題.同時(shí)，結(jié)合提出的聚焦強(qiáng)度閾值方法，可以顯著地改善聚焦反演或平滑聚焦反演中易出現(xiàn)的過度聚焦問題，更好地恢復(fù)異常體特征.

(3)在交叉梯度聯(lián)合反演中，若兩種物性結(jié)構(gòu)已近似一致，但均與真實(shí)模型有所偏差.此時(shí)，由于兩種物性反演的分步性和交替性，連續(xù)的交叉梯度結(jié)構(gòu)相似約束反而會限制反演模型向真實(shí)模型迭代的趨勢.而本文提出的間斷式聯(lián)合方法，通過交替進(jìn)行單獨(dú)反演和交叉梯度聯(lián)合反演，既能保證約束物性結(jié)構(gòu)相似，又能允許迭代模型響應(yīng)向觀測數(shù)據(jù)充分?jǐn)M合，可以獲得更加收斂的聯(lián)合反演結(jié)果.

(4)由于重力、重力梯度和大地電磁方法存在分辨率差異，交叉梯度聯(lián)合之前，應(yīng)使單獨(dú)反演得到的密度異常和電阻率異常在形態(tài)上盡量相近.通過靈活調(diào)整反演參數(shù)來控制單獨(dú)反演結(jié)果的異常發(fā)散程度，得到盡量大范圍的有效交叉梯度，有利于改善最終聯(lián)合反演結(jié)果.

本文限于條件僅采用了理論模型對算法進(jìn)行研究，目前尚未將其應(yīng)用于實(shí)測數(shù)據(jù)反演，接下來將進(jìn)一步對此進(jìn)行研究.

致謝感謝評審專家在百忙之中為本文修改提出的寶貴意見.