李桐林，張镕哲，樸英哲,*，金鋼燮，金長鮮，催邴敏

1 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春　130026，中國 2 金策工業(yè)綜合大學(xué)資源探測工學(xué)系，平壤　999093，朝鮮

?

部分區(qū)域約束下的交叉梯度多重地球物理數(shù)據(jù)聯(lián)合反演

李桐林1，張镕哲1，樸英哲1,2*，金鋼燮2，金長鮮2，催邴敏2

1 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春130026，中國 2 金策工業(yè)綜合大學(xué)資源探測工學(xué)系，平壤999093，朝鮮

交叉梯度聯(lián)合反演方法通過對多種地球物理模型實現(xiàn)結(jié)構(gòu)耦合，在巖石物性關(guān)系不確定的情況下，既能提高反演結(jié)果的可靠性，又能減少反演的多解性，還能減少不同方法解釋結(jié)果之間的矛盾．當不同的模型觀測數(shù)據(jù)覆蓋范圍不一致時，交叉梯度聯(lián)合反演通常需要取出重疊區(qū)域數(shù)據(jù)進行聯(lián)合反演，并且建模時還要擴展一些模型范圍．本文首先提出并實現(xiàn)了部分區(qū)域約束下的交叉梯度多重地球物理數(shù)據(jù)聯(lián)合反演算法；接著進行了算法的模型試算；最后，我們將該反演算法用于本溪—集安深部地質(zhì)調(diào)查重磁電綜合地質(zhì)地球物理解釋中．結(jié)果表明：該算法不但能在重疊區(qū)域內(nèi)很好地恢復(fù)結(jié)構(gòu)相似的模型，而且在非重疊區(qū)域與重疊區(qū)域的邊界處仍然可以得到平滑變化的模型；在本溪—集安10號剖面所獲得的結(jié)構(gòu)上相似的電阻率、密度及磁化率模型較好地反映了該區(qū)的深部地質(zhì)結(jié)構(gòu)，對于確定深部地質(zhì)體的性質(zhì)提供了有力的證據(jù)．

交叉梯度；聯(lián)合反演；大地電磁測深；初至波走時；重磁法

1　引言

在地球物理反演中，由于在空間和時間上的觀測局限性等原因，始終存在多解性問題(Tarantola，2004)．減少多解性的通常方法是應(yīng)用模型光滑約束(Constable et al．，1987)，其能夠有效地消除不光滑構(gòu)造信息，并得到可靠的反演結(jié)果(吳小平和徐果明，2000)．為了進一步減少反演結(jié)果不確定性和提高解釋精度，需要結(jié)合多種地球物理方法．不同地球物理數(shù)據(jù)的結(jié)合是基于相應(yīng)的物性之間的某些耦合關(guān)系．最簡單的形式是對同一種物性的不同探測方法數(shù)據(jù)的聯(lián)合解釋(Vozoff and Jupp，1975；Sasaki，1989)．不同物性的耦合關(guān)系可以分成巖石物性關(guān)系和結(jié)構(gòu)相似性關(guān)系．巖石物性關(guān)系因為建立了不同物性之間的函數(shù)關(guān)系(例如：在多孔介質(zhì)中電阻率與地震速度間巖石物性特征(Tillmann and St?cker，2000))所以在聯(lián)合反演解釋中起到明顯的作用(Moorkamp et al．，2011)，但是在很多情況下，這種參數(shù)直接關(guān)系是不確定的，而且其穩(wěn)定性也是未知的．結(jié)構(gòu)相似性關(guān)系反映的是不同物性的空間變化關(guān)系，其基礎(chǔ)是多種地球物理方法都能反映相同的地下構(gòu)造及地質(zhì)體(Gallardo and Meju，2011)．結(jié)構(gòu)耦合聯(lián)合反演最初是通過模型曲率來約束結(jié)構(gòu)相似性(Haber and Oldenburg，1997；Zhang and Morgan，1997)．基于模型曲率的結(jié)構(gòu)約束適合于邊緣檢測問題，但這種方法并沒有考慮模型梯度的方向．Gallardo和Meju(2003)提出的另一種方法是用不同物性梯度的叉乘來識別結(jié)構(gòu)邊界的相似性．這種交叉梯度聯(lián)合反演方法已經(jīng)被用到地球物理綜合解釋中(Gallardo and Meju，2004；Gallardo，2007；Gallardo et al．，2012；Abubakar et al．，2012；Fregoso and Gallardo，2009；Moorkamp et al．，2013；Linde et al．，2008；Um et al．，2014；彭淼等，2013；Hamdan et al．，2012)及水文問題(Lochbühler et al．，2013)解釋中．

交叉梯度聯(lián)合反演在結(jié)構(gòu)相似約束方式上有兩種思路，一種思路為尋找不同模型之間嚴密的結(jié)構(gòu)相似性，另一種思路為通過最小化結(jié)構(gòu)相似項的模來實現(xiàn)結(jié)構(gòu)約束(Gallardo and Meju，2011)．目前提出的嚴密結(jié)構(gòu)相似方法都只在相同模型范圍內(nèi)進行聯(lián)合反演．這意味著有以下兩個問題：

第一，對于單獨反演，地下模型的網(wǎng)格剖分規(guī)律取決于解釋什么樣的物理場．例如對直流電阻率法和大地電磁測深法反演，為了考慮邊界條件的影響，二維問題時需要沿著4個模型邊緣進行模型擴展來減少邊界對主要解釋區(qū)域的影響，然而初至波成像或重磁法反演不需要過于擴展模型．但在進行聯(lián)合反演時，密度、磁化率或者縱波速度模型必須跟著電阻率模型進行擴展．

第二，當一種探測方法的測線可能是其他探測方法測線的一部分時，此時只能用重疊區(qū)域的觀測數(shù)據(jù)進行聯(lián)合反演．例如在Hamdan等(2012)的野外實測中，初至波測線是ERT測線的一部分，所以在進行聯(lián)合反演時只用到ERT觀測數(shù)據(jù)的一部分和全部的初至波走時數(shù)據(jù)．

本文首先介紹了基于嚴密結(jié)構(gòu)相似的交叉梯度聯(lián)合反演方法，然后提出了部分區(qū)域約束下的交叉梯度多重地球物理數(shù)據(jù)聯(lián)合反演方法．最后通過模型試驗和野外數(shù)據(jù)應(yīng)用證明了該方法的合理性和實用性．

2　交叉梯度聯(lián)合反演方法

(1)

約束τjl=0:?j≠l，其中，m0j為mj的先驗?zāi)Ｐ蛥?shù)，Cdi和Cmj為觀測數(shù)據(jù)和模型向量的協(xié)方差，di為第i種探測方法的觀測數(shù)據(jù)，αj為阻尼系數(shù)．我們忽略了不同物性參數(shù)之間的巖石物性關(guān)系，只考慮結(jié)構(gòu)上相似程度．在二維的情況下交叉梯度函數(shù)

(2)

只在走向方向含有非零成分．聯(lián)合反演是對不同單位和不同量級的物性參數(shù)之間的耦合，因此，為了防止數(shù)值計算的不穩(wěn)定性，在式(2)中加入系數(shù)κ來標準化模型向量．一般情況下設(shè)定κj為向量mj的預(yù)期變化幅度．聯(lián)合反演之前，分別進行單獨反演，然后單獨反演結(jié)果的每一種模型參數(shù)的最大值和最小值的差即為對應(yīng)的κ值．為了簡單的描述公式推導(dǎo)，將式(1)以矩陣形式表示為

(3)

約束τ(m)=0，

其中

Cd=diag[Cd1,Cd2,…],

將式(3)的目標函數(shù)及約束項在當前模型附近進行泰勒級數(shù)展開，聯(lián)合反演可簡化為線性約束條件下的最小二乘問題．再用拉格朗日乘子方法(Tarantola，2004)，第二次模型更新式為

(4)

其中

ΔmLS為無交叉梯度約束的正則化最小二乘模型更新；B為交叉梯度函數(shù)的雅克比矩陣；拉格朗日乘子Λ為

(5)

ΔmLS和其協(xié)方差矩陣CLS為

(6)

(7)

其中

3　部分區(qū)域交叉梯度聯(lián)合反演

3.1部分區(qū)域交叉梯度聯(lián)合反演原理

部分區(qū)域聯(lián)合反演是指在多個模型的重疊區(qū)域里交叉梯度為零的約束條件下進行的聯(lián)合反演．因此目標函數(shù)跟上面的方法一致，只是式(1)的約束條件跟上面方法有區(qū)別．部分區(qū)域聯(lián)合反演的交叉梯度約束條件為

(8)

圖1　重疊模型區(qū)域內(nèi)的交叉梯度約束Fig.1　Crossgradient constraint in the overlapped region of models

(9)

將式(9)在當前模型附近進行泰勒級數(shù)展開得到

其中

第j種和第l種模型間的交叉梯度對第j種模型參數(shù)的雅克比矩陣為

(10)

(11)

整個交叉梯度雅克比矩陣B為

(12)

圖2　在包含電磁法數(shù)據(jù)的聯(lián)合反演中電阻率模型與其他物性模型的關(guān)系Fig.2　Relations between resistivity model and other geophysical models for joint inversions including electromagnetic data

圖3　部分區(qū)域交叉梯度聯(lián)合反演算法Fig.3　Joint inversion algorithm with sub-region crossgradient constraints

我們可以用式(12)及式(4)—(7)計算出第二次迭代的模型近似m′．式(11)中的B和式(4)中的BTΛ只在對應(yīng)重疊區(qū)域的模型單元內(nèi)存在非零元素項，但在式(4)中，由于受到平滑模型約束的CLS影響，CLSBTΛ對不在重疊區(qū)域的模型單元內(nèi)也存在非零數(shù)值，所以部分區(qū)域交叉梯度約束不僅對重疊區(qū)域元素而且對其他區(qū)域的元素都產(chǎn)生影響．換句話說聯(lián)合反演在其他區(qū)域能得到觀測數(shù)據(jù)擬合差和平滑模型約束的模型，而在重疊區(qū)域不但能得到以上兩種約束條件下的模型還能得出結(jié)構(gòu)相似的交叉梯度約束的模型．因此部分區(qū)域交叉梯度聯(lián)合反演方法不但利用所有觀測數(shù)據(jù)而且使模型在重疊區(qū)域滿足結(jié)構(gòu)相似性約束條件，還能在電磁法和重磁法及初至波走時的聯(lián)合反演時，不需要密度、磁化率及速度模型跟著電阻率模型擴展(圖2)．

3.2部分區(qū)域交叉梯度聯(lián)合反演算法

部分區(qū)域交叉梯度聯(lián)合反演的算法如圖3所示由二重迭代組成．外部迭代以一樣的比率逐漸減少阻尼系數(shù)，內(nèi)部迭代對固定的阻尼系數(shù)尋求滿足交叉梯度約束條件的解．在內(nèi)部迭代為了提高計算效率不再計算模型響應(yīng)fi(mj)和雅克比矩陣Aij(Gallardo et al．，2012)．

如果觀測數(shù)據(jù)類別和模型參數(shù)類別個數(shù)相同并且每一種觀測數(shù)據(jù)只跟一種物性有關(guān)，如圖2所示，電阻率和大地電磁數(shù)據(jù)、密度和布格異常、磁化率和磁異常及縱波速度和初至波走時的聯(lián)合反演時，因為Aij=0(?i≠j)所以矩陣A變成分塊對角矩陣．通過式(7)可以得到矩陣CLS也變成分塊對角矩陣

(13)

為此式(6)可改寫為

(14)

其中ΔmLSi為被第i種觀測數(shù)據(jù)和正則化約束決定的最小二乘模型．用式(5)計算拉格朗日乘子時，系數(shù)矩陣BCLSBT維數(shù)較大，然而矩陣CLS為分塊對角矩陣時可以用迭代解法(比如Minimal Residual Method(Barrett et al．，1994))避免構(gòu)成大型逆矩陣．

做聯(lián)合反演之前，對每一種探測方法做單獨反演，然后從單獨反演中得出阻尼系數(shù)的初值．阻尼系數(shù)的減少比率是通過阻尼系數(shù)在外部循環(huán)的最后迭代達到最小值得到的．根據(jù)我們的經(jīng)驗，αi最大值與最小值之間比率通常為100～1000，外部循環(huán)和內(nèi)部循環(huán)迭代次數(shù)是5～20、3～5次合適．

4　模型試驗與應(yīng)用

4.1理論模型(音頻大地電磁法和初至波走時數(shù)據(jù)的聯(lián)合反演)

為了展示本文方法的優(yōu)勢，構(gòu)建了與文獻Hamdan等(2012)的野外數(shù)據(jù)例子相似的模型(圖4)．電阻率模型為長度 8 km(水平距離1～9 km)、深度2 km，速度模型在位于水平距離約5～8 km的地方，深度為1.5 km．背景介質(zhì)設(shè)定為電阻率300 Ωm，速度隨著深度的變化線性增加的模型(圖5)．兩個矩形異常體的電阻率為30 Ωm及10000 Ωm，頂面埋深為100 m，底部埋深為300 m，其中第一個異常體不在速度模型區(qū)域里(圖4)．在速度模型區(qū)域里的第二個異常體的速度設(shè)定為6000 m·s-1．大地電磁測量數(shù)據(jù)為TE、TM兩種極化模式下的視電阻率和阻抗相位(9個測點(圖4中黑色菱形)、15個頻點(10～5000 Hz的對數(shù)間隔)，圖6a所示的其中TE模式視電阻率的擬斷面)，地震數(shù)據(jù)是由10個地震發(fā)射器(圖5中倒三角)和30個接收器(圖5中正三角)組成的觀測系統(tǒng)記錄的初至波走時(圖6b)，并在所有數(shù)據(jù)中加入了2%的隨機噪聲．反演中初始電阻率模型為300 Ωm的均勻半空間，對于初始速度模型，為了保證折射波傳播，設(shè)定為速度沿著深度方向遞增的介質(zhì)．用于大地電磁及初至波走時正演與雅克比矩陣計算的程序為Occam2DMT(deGroot-Hedlin and Constable，1990)和FAST(Zelt and Barton，1998)．

圖4　理論模型電阻率網(wǎng)格剖分為灰色網(wǎng)格部分，虛線外部為擴展的電阻率單元，黑色實線為速度模型邊緣，黑色為異常體的邊緣．Fig.4　Synthetic modelThe mesh of resistivity model is plotted with gray lines.Outside of dashed lines are padded resistivity cells.Black solid lines are the edge of velocity model.The designed anomalies are outlined with black rectangles．

圖5　速度模型射線分布Fig.5　Ray distribution in the velocity model

反演時，為了比較本文方法和常規(guī)交叉梯度聯(lián)合反演方法，將速度模型的范圍擴充至電阻率模型的范圍，構(gòu)建了范圍一致的電阻率模型和速度模型．再考慮FAST方法對網(wǎng)格剖分的均勻大小要求和電磁場邊界條件的影響，需要將電阻率模型和速度模型同時沿著水平方向和深度方向進行適當?shù)哪Ｐ蛿U展(分別10個模型單元)．對同一區(qū)域模型及不同區(qū)域模型分別進行了單獨反演和聯(lián)合反演．

無論是同一區(qū)域模型還是不同區(qū)域模型，電阻率單獨反演都很明顯的恢復(fù)了兩個異常體(圖7a和圖8a)，但由于大地電磁測點之間的距離較大，第二個高阻異常體的寬度小于實際模型寬度．速度模型的射線分布主要集中在水平距離5～8 km和深度300 m以上(圖5)，其他范圍內(nèi)沒有射線分布，導(dǎo)致反演會在射線分布以外的區(qū)域敏感度降低．在速度單獨反演模型中，由于速度異常體位于射線分布區(qū)域內(nèi)，所以在地震單獨反演剖面上可以很明顯地看到高速異常體(圖7b和圖8b)．

在聯(lián)合反演的電阻率剖面上，無論同一區(qū)域模型(圖7c)還是不同區(qū)域模型(圖8c)高阻異常體的寬度都比單獨反演更接近模型真實值，這說明電阻率模型受到了速度模型結(jié)構(gòu)的約束．對低阻異常體而言，因為速度模型上沒有異常體從而未受到結(jié)構(gòu)約束，所以聯(lián)合反演模型比單獨反演的沒有改善．聯(lián)合反演的速度模型(圖7d和圖8d)跟單獨反演進行對比，速度異常體與背景介質(zhì)的分離更明顯，這說明了速度模型也受到了電阻率模型的高阻異常體結(jié)構(gòu)的影響．在同一區(qū)域模型和不同區(qū)域模型的聯(lián)合反演中，速度模型在深度1 km以下部分的速度分布差異(圖7d和圖8d)可以推測為不同模型范圍的電阻率和速度靈敏度和模型平滑約束共同引起的．

圖6　理論模型響應(yīng)(a)TE模式視電阻率;(b)初至波走時.Fig.6　Synthetic model responses(a)TE mode apparent resistivity;(b)First-arrival traveltime.

圖7　擴充模型的單獨反演(a和b)與聯(lián)合反演(c和d)結(jié)果．(a和c)電阻率反演模型，(b和d)速度反演模型Fig.7　Separate (a and b)and joint (c and d)inversion results of extended models.(a and c)resistivity models,(b and d)velocity models

圖8　圖4模型的單獨反演(a和b)與聯(lián)合反演(c和d)結(jié)果．(a和c)電阻率反演模型，(b和d)速度反演模型Fig.8　Separate (a and b)and joint (c and d)inversion results of Fig.4 models.(a and c)resistivity models,(b and d)velocity models

圖9　本溪—吉安10號測線觀測數(shù)據(jù)及反演模型響應(yīng)(a)觀測TM視電阻率；(b)圖9a中電阻率模型的TM視電阻率響應(yīng)；(c)和(d)分別為觀測布格異常和磁異常及圖9a中密度、磁化率模型的正演響應(yīng).Fig.9　Observations and model responses of the survey line No.10 in Benxi-Ji′an area(a)Observed TM mode apparent resistivity；(b)Response from the resistivity model of Fig.9a；(c)and(d)Observed Bouguer and TMI (Total Magnetic Intensity)anomaly,and responses from the density and susceptibility models of Fig.9a．

總體上，常規(guī)交叉梯度聯(lián)合反演方法和本文方法都比單獨反演更好地恢復(fù)異常體．但兩個方法中存在差別．常規(guī)方法中，電阻率模型為了符合速度模型網(wǎng)格剖分規(guī)律用均勻大小的單元擴展了模型，而本文方法中，電阻率模型隨著水平距離和深度的增加用逐漸變大的單元擴展模型，這是更符合大地電磁正反演精度的要求．在第11次迭代時，兩種反演方法的擬合差都達到期望值1以下，反演時間分別為15分52秒，12分14秒(Intel Core2Duo E7200 2.53 GHz)，時間的增大是由速度模型擴展引起的.

4.2野外數(shù)據(jù)應(yīng)用

本溪—集安地區(qū)橫跨兩個不同地質(zhì)構(gòu)造單元，一是太古代古陸核，由太古代變質(zhì)巖組成，其控礦作用明顯；二是元古代遼吉古裂谷帶，其內(nèi)礦產(chǎn)資源十分豐富．研究區(qū)太古代的含礦建造除了在研究區(qū)的北部出露地表外，大部都處于隱伏狀態(tài)，元古代遼吉古裂谷帶沉積巨厚并且被古生代、中新生代沉積蓋層覆蓋，且無論是太古代含礦建造還是遼吉古裂谷帶都經(jīng)歷了多期構(gòu)造與巖漿侵入活動．查明該區(qū)深部地質(zhì)結(jié)構(gòu)與構(gòu)造，將會對研究區(qū)找礦具有非常重要的意義(張宏嘉等，2013)．為此在研究區(qū)共完成了十幾條重、磁、大地電磁綜合剖面．使用的儀器是加拿大鳳凰公司V5-2000系列電磁儀、美國BURRIS型高精度金屬彈簧重力儀和加拿大GEM公司GSM-19T微機質(zhì)子磁力儀．此外對研究區(qū)涉及的主要地質(zhì)體的物性也進行了測量．

用本文方法處理了其中10號測線數(shù)據(jù)(圖9)．測線方向大致從北西到南東，測線長度為97 km．在測線上共有大地電磁法48個測點，測點之間平均距離為2 km，測量數(shù)據(jù)中采取了15個頻點(360.0～0.02 Hz的對數(shù)間隔)的TM極化模式視電阻率(圖9a)．從重磁觀測數(shù)據(jù)得到布格異常和磁異常，其中采取189個測點數(shù)據(jù)(圖9c和圖9d中的紅色圓形)．

密度與磁化率模型的深度設(shè)定為20 km．對于電阻率模型，考慮電磁場的衰減，用對數(shù)間隔單元沿著水平方向和深度方向擴展了模型(圖10)．主要解釋區(qū)域(圖10中虛線表示的區(qū)域)網(wǎng)格單元大小設(shè)定為寬1100 m，高200～4000 m．反演初始模型為300 Ωm、0 g·cm-3及0 SI的均勻介質(zhì)．重磁正演及雅克比矩陣計算根據(jù)多邊形異常體公式(Singh，2002)．

由于此次研究主要關(guān)心深度10 km以內(nèi)的地質(zhì)結(jié)構(gòu)，因此圖11a中只畫出了深度到10 km的剖面．聯(lián)合反演過程在第6次迭代時以1.52(電)、0.81(重)、1.38(磁)的數(shù)據(jù)擬合差結(jié)束．最終模型響應(yīng)分別表示為圖9b的擬斷面(TM視電阻率)和圖9c(布格異常)、圖9d(磁異常)的藍色實線．

聯(lián)合反演獲得了結(jié)構(gòu)上相似的電阻率、密度及磁化率模型(圖11a)．根據(jù)物性范圍，可清晰地劃分出不同的地質(zhì)體．結(jié)合地表地質(zhì)情況(圖11b)和物性測量數(shù)據(jù)作進一步解釋，得到每個地質(zhì)體的位置和物性范圍，推斷了各地質(zhì)體的屬性以及斷層構(gòu)造．綜合解釋剖面圖(圖11c)與10號測線附近地表地質(zhì)(圖11b)和聯(lián)合反演剖面(圖11a)基本吻合，其中水平位置72～76 km、深度1～6 km的地質(zhì)體沒有出現(xiàn)地表地質(zhì)圖上，但是根據(jù)其物性范圍和在地質(zhì)圖外部浪子山組鄰接的地層及侵入巖屬性，推測該地質(zhì)體為古元古代蓋縣組．10號測線剖面可以分出兩個地質(zhì)單元，第一段為0～50 km區(qū)間的桓仁地塊，主要巖性由早白堊紀梨樹溝組、小嶺組火山沉積層、古元古代大石橋組組成；第二段為50～95 km區(qū)間的集安地塊，主要巖性為早白堊紀小嶺組火山沉積層、古元古代蓋縣組、浪子山組、早白堊紀侵入巖．測線下伏巖體為古元古代時期生成的遼吉花崗巖．根據(jù)綜合解釋剖面和聯(lián)合反演結(jié)果推測的每個地質(zhì)體的物性范圍如表1所示．

圖10　聯(lián)合反演模型網(wǎng)格剖分虛線為密度及磁化率模型的邊緣，灰色網(wǎng)格部分為電阻率模型剖分，交叉梯度約束只適用于虛線區(qū)域內(nèi)．Fig.10　Mesh generation of joint inversion modelDashedline is the edge of density and susceptibility model.Gray line is the mesh of resistivity model.Crossgradient constraint is only applied to the dashed line rectangle.

圖11　10號測線聯(lián)合反演結(jié)果(a)、地質(zhì)圖(b)及聯(lián)合反演綜合解釋剖面圖(c)Fig.11　(a)Joint inversion results of the survey line No.10，(b)geological map and (c)comprehensively interpreted section of the joint inversion result

地質(zhì)體電阻率log10Ωm密度異常g·cm-3磁化率異常SI地質(zhì)屬性K1l2.0～2.3-0.03～0.01-0.017～0.000早白堊紀梨樹溝組Jkx1.7～2.2-0.12～0.000.000～0.030晚侏羅紀—早白堊紀小嶺組Pt1gx3.7～3.9-0.02～-0.01-0.001～0.000古元古代蓋縣組Pt1d1.8～2.50.02～0.14-0.029～-0.004古元古代大石橋組Pt1l3.6～3.90.00～0.050.000～0.007古元古代浪子山組K1Fh3.5～3.6-0.02～-0.010.003～0.025早白堊紀花崗巖K1Lj3.5～3.8-0.01～0.030.000～0.015早白堊紀連江輝石閃長巖Pt2κγ3.3～3.5-0.06～-0.01-0.021～0.001中元古代堿長花崗巖Pt1Hrγ3.3～3.8-0.03～-0.02-0.035～-0.013古元古代桓仁花崗雜巖

5　結(jié)論

針對聯(lián)合反演中不同地球物理模型范圍不一致的問題，本文提出并實現(xiàn)了部分區(qū)域約束下的交叉梯度多重聯(lián)合反演，并進行了模型試算與應(yīng)用試驗，得到了如下結(jié)論：

(1)本文提出的算法僅在重疊區(qū)域內(nèi)應(yīng)用交叉梯度約束，克服了所有模型區(qū)域需要一致的限制，能對不相同模型范圍進行結(jié)構(gòu)耦合聯(lián)合反演．

(2)理論模型結(jié)果表明：當觀測數(shù)據(jù)覆蓋的區(qū)域不完全一致時，本文算法不但可以用到所有觀測數(shù)據(jù)，而且能在重疊區(qū)域內(nèi)很好地恢復(fù)結(jié)構(gòu)耦合的結(jié)果．除此之外，由于平滑模型約束的影響，在非重疊區(qū)域與重疊區(qū)域的邊界處仍然可以得到平滑變化的模型．

(3)應(yīng)用試驗得到的電阻率、密度及磁化率模型結(jié)構(gòu)上很相似，較好地反映了研究地區(qū)的深部構(gòu)造，為確定深部地質(zhì)體性質(zhì)提供了有力證據(jù)．

致謝感謝匿名評審專家為本文修改提出的細致而寶貴的意見．

References

Abubakar A，Gao G，Habashy T M，et al．2012．Joint inversion approaches for geophysicalelectromagnetic and elastic full-waveform data．Inverse Problems，28(5)：055016，doi：10.1088/0266-5611/28/5/055016．

Barrett R，Berry M，Chan T F，et al．1994．Templates for the Solution of Linear Systems:Building Blocks for Iterative Methods．Philadelphia：SIAM．

Constable S C，Parker R L，Constable C G．1987．Occam′s Inversion：A practical algorithm for generating smooth models from electromagnetic sounding data．Geophysics，52(3)：289-300，doi:10.1190/1.1442303．

deGroot-Hedlin C D，Constable S C．1990．Occam′s Inversion to Generate Smooth Two Dimensional Models from Magnetotelluric Data．Geophysics，55(12)：1613-1624，doi：10.1190/1.1442813．

Fregoso E，Gallardo L A．2009．Cross-gradients joint 3D inversion with applicationsto gravity and magnetic data．Geophysics，74(4)：L31-L42，doi：10.1190/1.3119263．Gallardo L A，Meju M A．2003．Characterization of heterogeneousnear-surface materials by joint 2D inversion of DC resistivity and seismicdata．Geophys．Res．Lett．，30(13)，doi：10.1029/2003GL017370．Gallardo L A，Meju M A．2004．Joint two-dimensional DC resistivityand seismic travel time inversion with cross-gradients constraints．J．Geophys．Res．，109(B3)：B03311，doi：10.1029/2003JB002716．

GallardoLA．2007．Multiple cross-gradient joint inversion for geospectralimaging．Geophys．Res．Lett．，34(19)：L19301-5，doi：10.1029/2007GL030409．Gallardo LA，MejuMA．2011．Structure-coupled multiphysics imagingin geophysical sciences．Rev．Geophys．，49(1)：RG1003，doi：10.1029/2010RG000330．Gallardo LA，F(xiàn)ontesSL，MejuMA，et al．2012．Robust geophysical integration through structure-coupled joint inversionand multispectral fusion of seismic reflection，magnetotelluric，magnetic，and gravity images:Example from Santos Basin，offshore Brazil．Geophysics，77(5)：B237-B251，doi：10.1190/GEO2011-0394.1．Haber E，Oldenburg D．1997．Joint inversion:A structuralapproach．Inverse Problems，13(1)：63-77，doi：10.1088/0266-5611/13/1/006．

Hamdan H，Economou N，Kritikakis G，et al．2012．2D and 3D imaging of the metamorphic carbonates at Omalos plateau/polje，Crete，Greece by employing independent and joint inversion on resistivity and seismic data．International Journal of Speleology，41(2)：199-209，doi：10.5038/1827-806X.41.2.7．

Linde N，Tryggvason A，Peterson J E，et al．2008．Joint inversion of crosshole radar and seismic traveltimesacquired at the South Oyster Bacterial Transport Site．Geophysics，73(4)：G29-G37，doi：10.1190/1.2937467．

Lochbühler T，Doetsch J，Brauchler R，et al．2013．Structure-coupled joint inversion of geophysical and hydrological data．Geophysics，78(3)：ID1-ID14，doi：10.1190/GEO2012-0460.1．

Moorkamp M，Heincke B，Jegen M，et al．2011．A framework for 3D joint inversion of MT，gravity and seismicrefraction data．Geophys．J．Int．，184(1)：477-493，doi：10.1111/j.1365-246X.2010.04856.x．

Moorkamp M，Roberts A W，Jegen M，et al．2013．Verification of velocity-resistivity relationships derived fromstructural joint inversion with borehole data．Geophys．Res．Lett．，40(14)：3596-3601，doi：10.1002/grl.50696．

Peng M，Tan H D，Jiang M，et al．2013．Three-dimensional joint inversion of magnetotelluric and seismic traveltime data with cross-gradient constraints．Chinese J．Geophys．(in Chinese)，56(8)：2728-2738，doi：10.6038/cjg20130821．

Sasaki Y．1989．Two-dimensional joint inversion of magnetotelluric and dipole-dipole resistivity data．Geophysics，54(2)：254-262，doi：10.1190/1.1442649．

Singh B．2002．Simultaneous computation of gravity and magnetic anomalies resulting from a 2-D object．Geophysics，67(3)：801-806，doi：10.1190/1.1484524．

TarantolaA．2004．Inverse Problem Theory and Methods for Model Parameter Estimation．Philadelphia：SIAM．

Tillmann A，St?cker T．2000．A new approach for the joint inversion of seismic and geoelectric data．//Extended Abstracts，63rdEAGE Conference and Technical Exhibition．Amsterdam：EAGE．

Um E S，Commer M，Newman G A．2014．A strategy for coupled 3D imaging of large-scale seismic andelectromagnetic data sets:Application to subsalt imaging．Geophysics，79(3)：ID1-ID13，doi：10.1190/GEO2013-0053.1．

Vozoff K，Jupp DLB．1975．Joint inversion of geophysical data．Geophys．J．Int.，42(3)：977-991，doi：10.1111/j.1365-246X.1975.tb06462.x．

Wang J,Meng X H,Chen Z X.2013.The theory of cross-gradient and its application in geophsical joint inversion.Chinese J.Geophys.(inChinese),28(4):2094-2103,doi:10.6038/pg20130454.Wu X P，Xu G M．2000．Study on 3-D resistivity inversion using conjugate gradient method．Chinese J．Geophys．(in Chinese)，43(3)：420-427．

Zelt C A，Barton P J．1998．Three-dimensional seismic refraction tomography：A comparison of two methods applied to datafrom the Faeroe Basin．J．Geophys．Res．，103(B4)：7187-7210，doi：10.1029/97JB03536．

Zhang H J，Cao L，Li F W，et al．2013．GME comprehensive explanation of 3D geological structure of Benxi-Ji′an region．Jilin Geology (in Chinese)，32(1)：88-94．

Zhang J，Morgan FD．1997．Joint seismic and electrical tomography.//Symposium on the Applications of Geophysicsto Engineering and Environmental Problems．Keystone，Colo：EEGS，391-396．

附中文參考文獻

彭淼，譚捍東，姜枚等．2013．基于交叉梯度耦合的大地電磁與地震走時資料三維聯(lián)合反演．地球物理學(xué)報，56(8)：2728-2738，doi：10.6038/cjg20130821．

王俊,孟小紅,陳召曦等.2013.交叉梯度理論及其在地球物理聯(lián)合反演中的應(yīng)用.地球物理學(xué)進展,28(4):2094-2103,doi:10.6038/pg20130454.

吳小平，徐果明．2000．利用共軛梯度法的電阻率三維反演研究．地球物理學(xué)報，43(3)：420-427．

張宏嘉，曹亮，李福文等．2013．本溪—集安地區(qū)三維地質(zhì)結(jié)構(gòu)重磁電綜合解釋．吉林地質(zhì)，32(1)：88-94．

(本文編輯劉少華)

Multiple joint inversion of geophysical data with sub-region crossgradient constraints

LI Tong-Lin1，ZHANG Rong-Zhe1，PAK Yong-Chol1,2*，KIM Gang-Sop2，KIM Jang-Son2，CHOE Byong-Min2

1 College of Geo-Exploration Sciences and Technology,Jilin University,Changchun 130026,China 2 Resource Exploration and Engineering Department,Kimchaek University of Technology,Pyongyang 999093, Democratic People′s Republic of Korea

Crossgradient joint inversions,achieving structure coupling of multiple geophysical models,improve the reliability of inversion results,reduce the uncertainty of inversions,and reduce conflicts between results of different methods,in the case of unknown petrophysical relationships.When coverage of different observations is inconsistent,the crossgradient joint inversion usually need extract data corresponding to the overlapped region,and also expands the scope of some models.This paper proposed and implemented a multiple joint inversion of geophysical data with sub-region crossgradient constraints and followed by a synthetic model example.We have applied this method to comprehensive geological and geophysical interpretation of MT,gravity and magnetic field data sets for deep geological survey of Benxi-Ji′an.The results show that this algorithm not only can well restore similar structure model in the overlapped region,but the model parameters are varied smoothly on the boundary between overlapped and non-overlapped regions.The resistivity,density and susceptibility models,obtained from the joint inversion of Benxi-Ji′an survey line No.10,are very similar in their structures,well reflect the deep geological structure of this area,and provide strong evidences to determine the nature of deep geological bodies.

Crossgradient；Joint inversion；Magnetotellurics；First-arrival traveltime；Gravity and magnetic method

李桐林，張镕哲，樸英哲等.2016.部分區(qū)域約束下的交叉梯度多重地球物理數(shù)據(jù)聯(lián)合反演.地球物理學(xué)報，59(8):2979-2988,

10.6038/cjg20160821.

Li T L,Zhang R Z,Pak Y C,et al.2016.Multiple joint inversion of geophysical data with sub-region crossgradient constraints.Chinese J.Geophys.(in Chinese),59(8):2979-2988,doi:10.6038/cjg20160821.

中國地質(zhì)調(diào)查局項目(12120113098400)，國家重大科學(xué)儀器設(shè)備開發(fā)專項(2011YQ05006009).

李桐林，男，1962年生，教授，主要從事電磁法理論與應(yīng)用教學(xué)和研究.E-mail:lilaoshizh@163.com

樸英哲，男，朝鮮人，1980年生，博士研究生，主要從事電磁正反演理論與應(yīng)用研究，E-mail:piaoyz525@163.com

10.6038/cjg20160821

P318，P319

2015-12-12，2016-05-15收修定稿