李兆祥， 譚捍東， 付少帥， 唐飛， 朱德朋， 高敬語

1 中國地質(zhì)大學(xué)（北京）地球物理與信息技術(shù)學(xué)院， 北京 1000832 江蘇省地質(zhì)勘查技術(shù)院， 南京 210000

?

基于交叉梯度約束的時間域激發(fā)極化法二維同步反演

李兆祥1,2， 譚捍東1*， 付少帥1， 唐飛1， 朱德朋1， 高敬語1

1 中國地質(zhì)大學(xué)（北京）地球物理與信息技術(shù)學(xué)院， 北京 1000832 江蘇省地質(zhì)勘查技術(shù)院， 南京 210000

時間域激發(fā)極化法（Time-domain induced polarization method，簡稱為TDIP）已有的反演算法采用的是分步反演的思路，即先由視電阻率資料反演電阻率，固定電阻率再由視極化率資料反演極化率，這樣就存在極化率結(jié)果嚴重依賴于電阻率反演結(jié)果的問題.為了有效解決這一問題，本文實現(xiàn)了TDIP二維數(shù)據(jù)空間分步反演算法，提出了基于交叉梯度約束的TDIP二維同步反演策略，實現(xiàn)了交叉梯度約束的電阻率和極化率二維同步反演算法.分別用電阻率和極化率結(jié)構(gòu)一致和不一致的二維模型合成數(shù)據(jù)進行了分步和同步反演試算，對不同模型試算結(jié)果進行了對比分析.結(jié)果表明：對于電阻率和極化率結(jié)構(gòu)一致和不一致模型，同步反演結(jié)果比分步反演結(jié)果能更好地確定異常體的空間分布范圍，反演得到的電阻率和極化率值更接近真值.理論模型算例表明本文提出的同步反演算法有效解決了分步反演的問題，優(yōu)于分步反演算法，具有更好的實用性.

時間域激發(fā)極化法； 數(shù)據(jù)空間反演； 交叉梯度； 同步反演

1 引言

近年來，隨著地球物理各種方法的正反演技術(shù)日趨成熟，利用多種地球物理參數(shù)對同一地質(zhì)體進行綜合解釋已經(jīng)得到了廣泛的應(yīng)用（Aric et al．，1997；姜枚等，2012a,2012b；高德章等，2004；胡平等，2010）.但由于受制于單方法本身的多解性，利用單方法單參數(shù)反演結(jié)果進行綜合解釋往往很難獲得較為統(tǒng)一的地質(zhì)-地球物理模型.為了解決這一問題，國內(nèi)外學(xué)者對多方法多參數(shù)聯(lián)合反演方法進行了深入研究（Moorkamp et al．，2007；彭淼等，2012；Colombo and De Stefano，2007；于鵬等，2009）.其中，采用模型結(jié)構(gòu)耦合的聯(lián)合反演已經(jīng)成為國際上的趨勢（Haber and Oldenburg，1997），由Gallardo和Meju（2003）提出的交叉梯度（cross-gradient）耦合方法是比較有代表性的成果.Gallardo和Meju（2007）進行了大地電磁和折射地震二維聯(lián)合反演算法的研究，通過在反演的目標函數(shù)中引入交叉梯度函數(shù)來約束不同的模型參數(shù)，比傳統(tǒng)的直接關(guān)系耦合更具有普遍適用性.國內(nèi)學(xué)者也有相關(guān)的研究成果，如彭淼等（2013）實現(xiàn)了大地電磁和天然地震的三維聯(lián)合反演算法.

地球物理反演方法很多，主要包括最小二乘法、共軛梯度法、OCCAM法以及基于現(xiàn)代優(yōu)化算法技術(shù)的組合最優(yōu)化算法等.其中OCCAM法反演得到的模型綜合了數(shù)據(jù)擬合差和光滑度兩項要素（Constable et al．，1987），數(shù)據(jù)空間法（data-space inversion）是在OCCAM法基礎(chǔ)上發(fā)展起來的一種反演方法（Siripunvaraporn and Egbert，2000; Siripunvaraporn et al.，2005）.

TDIP是以巖、礦石的激電效應(yīng)差異為物質(zhì)基礎(chǔ)，通過觀測和研究大地激電效應(yīng)，來探查地下地質(zhì)情況的一種電法分支方法.該方法在尋找金屬礦、能源和工程地質(zhì)等不同領(lǐng)域均有著廣泛的應(yīng)用，尤其在尋找浸染狀礦產(chǎn)資源方面效果顯著（傅良魁，1982）.國內(nèi)外學(xué)者對TDIP二維反演算法進行了大量研究，但是都是采用分步反演，即先反演出電阻率之后，在所反演出的電阻率模型基礎(chǔ)之上再反演極化率參數(shù)（Douglas and Li，1994，1999），這樣就存在著極化率反演結(jié)果嚴重依賴于電阻率反演結(jié)果等問題.

為了解決上述問題，本文在反演算法中引入了交叉梯度約束，提出了TDIP電阻率和極化率二維同步反演策略，通過設(shè)計電阻率和極化率物性結(jié)構(gòu)一致模型和不一致的模型進行測試，討論該同步反演算法的可靠性和適用性.

2 TDIP二維正反演算法

2.1 TDIP二維正演算法

本文采用等效電阻率法（Douglas and Li，1994）進行TDIP正演數(shù)值模擬.等效電阻率法以電阻率法的正演計算為基礎(chǔ)，該方法大體可以分成三步完成：首先根據(jù)各個模型節(jié)點上的電阻率和極化率數(shù)值，計算出節(jié)點上的等效電阻率；然后進行等效電阻率的正演模擬（把電阻率法正演過程略微修改即可），計算出等效視電阻率；最后根據(jù)電阻率正演和等效電阻率正演得到的視電阻率和等效視電阻率，即可計算出視極化率，這樣就完成了極化率數(shù)據(jù)數(shù)值模擬的整個過程.本文采用的電阻率法正演算法為有限單元法.

2.2 TDIP二維分步反演算法

Siripunvaraporn等人在OCCAM的基礎(chǔ)上進行了改進，通過一系列數(shù)學(xué)變換，將M（模型參數(shù)個數(shù)）行M列的模型空間的叉積矩陣轉(zhuǎn)換為N（觀測數(shù)據(jù)個數(shù)）行N列的數(shù)據(jù)空間叉積矩陣，演變成了數(shù)據(jù)空間反演法.對于二維和三維反演算法來說，N往往比M小得多，并且該方法不需要求取模型協(xié)方差矩陣的逆矩陣，因而大大減少了所需的計算時間和計算機內(nèi)存，提高了反演效率，是一種實用的反演算法.因此，本文采用的電阻率法反演方法即為數(shù)據(jù)空間反演（SiripunvarapornandEgbert，2000;Siripunvarapornetal.，2005）.

本文所采用的TDIP反演方法也是數(shù)據(jù)空間反演法，其基本過程與電阻率法反演一致，只需要完成TDIP反演的雅可比矩陣求取即可.在電阻率法反演的雅可比矩陣基礎(chǔ)上，可以計算得到TDIP反演的雅可比矩陣（DouglasandLi，1994，1999）.

3 引入交叉梯度約束的TDIP二維同步反演

3.1 同步反演目標函數(shù)

本文的同步反演核心思想是，每次迭代中同步更新電阻率和極化率參數(shù)，因此只定義一個目標函數(shù)，同時反演電阻率和極化率，并加入交叉梯度約束項.目標函數(shù)具體形式如下所示：

+λ2Φm（mη）+βΦcg（mρ,mη）,

（1）

式中：

（2）

（3）

Φcg=t（mρ,mη）Tt（mρ,mη）,

（4）

其中，Φd是數(shù)據(jù)擬合項，Φm是模型光滑項，Φcg是交叉梯度項；λ1和λ2是正則化參數(shù)，β是交叉梯度的加權(quán)因子，本文選取的加權(quán)因子和正則化參數(shù)在分步反演和同步反演的所有算例中是一致的；mρ和mη分別為電阻率和極化率模型向量，m0為初始模型，Cm為模型協(xié)方差矩陣，Cd為數(shù)據(jù)協(xié)方差矩陣，起到利用數(shù)據(jù)誤差給觀測數(shù)據(jù)加權(quán)的作用；dobs為實測數(shù)據(jù)向量，f（·）為正演算子；t（·）為交叉梯度算子（GallardoandMeju，2003，2007），在同步反演中引入交叉梯度函數(shù)可以實現(xiàn)電阻率參數(shù)和極化率參數(shù)的相互約束.

電阻率和極化率交叉梯度定義如下：

（5）

（6）

采用五點中心差分格式，離散式（6），得到：

其中，i和j分別表示網(wǎng)格z方向下標和x方向下標.

對目標函數(shù)式（1）求其關(guān)于模型參數(shù)的偏導(dǎo)數(shù)令其等于0，并且考慮到公式（GallardoandMeju，2007）

（8）

其中B為交叉梯度函數(shù)相對模型參數(shù)的偏導(dǎo)數(shù)矩陣，mρ0和mη0分別是電阻率參數(shù)和極化率參數(shù)的初始模型矩陣，參照數(shù)據(jù)空間反演的數(shù)學(xué)變換式，得到模型迭代更新的公式：

mk+1-m0=

（9）

式中Jk是第k次反演迭代過程中的雅可比矩陣，k=dobs-f[mk]+Jkmk，mk是第k次反演得到的模型向量.

式（9）即為TDIP同步反演第k+1次迭代模型公式.

由此可見，與分步反演不同，同步反演目標函數(shù)中包含了電阻率項和極化率項，并且加入了交叉梯度函數(shù)，這樣每迭代一次都同步更新電阻率模型和極化率模型參數(shù)，并且計算交叉梯度耦合.

3.2 同步反演流程

本文所實現(xiàn)的TDIP同步反演算法的流程大體如下：

（1）設(shè)置初始模型，輸入初始電阻率模型和初始極化率模型；

（2）計算交叉梯度函數(shù)，并同時實現(xiàn)電阻率和極化率的模型更新；

（3）計算視電阻率和視極化率數(shù)據(jù)擬合差之和；

（4）判斷擬合差之和是否滿足目標函數(shù)設(shè)置的收斂條件（收斂條件為兩次的擬合差之差小于10-6），不設(shè)置最大迭代次數(shù)，本文所列出的所有分步反演和同步反演算例的收斂條件都是一致的.擬合差的計算公式為

其中，dobs為含有高斯誤差的觀測數(shù)據(jù)，df為正演計算數(shù)據(jù)，e為數(shù)據(jù)誤差項.當模型能夠完全擬合不含誤差的觀測數(shù)據(jù)時擬合差為1,如果不滿足，迭代繼續(xù)，返回步驟（2）重新計算；如果滿足收斂條件，則完成同步反演，跳出迭代，輸出反演結(jié)果.

4 算例分析

基于以上理論部分，本文實現(xiàn)了TDIP二維同步反演算法.為了檢驗算法的正確性和可靠性，同時對比同步反演結(jié)果與分步反演結(jié)果的差別，本文進行了理論模型合成數(shù)據(jù)的反演，分別設(shè)計了電阻率和極化率結(jié)構(gòu)一致和不一致的二維模型，下面依次作分析討論.

4.1 電阻率和極化率結(jié)構(gòu)一致模型

本文設(shè)計了雙棱柱高阻極化體的物性結(jié)構(gòu)一致模型.左側(cè)模型異常體電阻率為1000 Ωm，極化率為15%，右側(cè)模型異常體電阻率為500 Ωm，極化率為10%，二者尺寸均為40 m×30 m，頂面埋深均為40 m，二者間距為80 m；異常體埋藏于電阻率為100 Ωm，極化率為2%的均勻半空間中，具體可見圖1a和圖1d所示.觀測系統(tǒng)與前文一致，在正演數(shù)據(jù)中加入2.0%的高斯隨機誤差.對模型合成數(shù)據(jù)分別進行分步反演和同步反演，初始模型采用電阻率為100 Ωm，極化率為1%的均勻半空間，分步反演的結(jié)果見圖1b和圖1e，同步反演的結(jié)果見圖1c和圖1f.

對反演結(jié)果進行分析，結(jié)果顯示：電阻率的分步反演結(jié)果對棱柱體的空間形態(tài)有一定的恢復(fù)，但是左側(cè)棱柱體電阻率恢復(fù)到150 Ωm左右，與電阻率真值（1000 Ωm）有著明顯的差距，右側(cè)棱柱體電阻率恢復(fù)到130 Ωm左右，與電阻率真值（500 Ωm）差距明顯；極化率的分步反演結(jié)果能恢復(fù)出異常體的大致空間形態(tài)，但是對左棱柱體極化率的真實數(shù)值（15%）最高只能恢復(fù)到3%左右，對右棱柱體極化率的真實數(shù)值（10%）最高只能恢復(fù)到3%左右.相比于分步反演，二者的埋深均比實際埋深偏淺.圖1c和圖1f中的同步反演結(jié)果對棱柱體的空間形態(tài)有更好的恢復(fù)，左棱柱體的反演電阻率可達到200 Ωm左右，右棱柱體的反演電阻率可達到180 Ωm左右，均比分步反演的結(jié)果更加接近于真實電阻率數(shù)值；而極化率的同步反演結(jié)果可以恢復(fù)出兩個異常體的埋深，并且兩個異常體的極化率均恢復(fù)到6.5%左右，相比于分步反演的結(jié)果改善顯著.

電阻率法分步反演迭代8次，擬合差從8.95下降到1.45，TDIP分步反演迭代8次，擬合差從26.96下降到1.81，同步反演迭代了10次，最后收斂時視電阻率的數(shù)據(jù)擬合差為1.37，視極化率的數(shù)據(jù)擬合差為1.54，二者均小于分步反演結(jié)果的數(shù)據(jù)擬合差.

得到分步反演和同步反演的結(jié)果之后，計算電阻率模型和極化率模型交叉梯度數(shù)值，對于該組模型而言，分步反演結(jié)果的交叉梯度數(shù)值為9.87×10-5，同步反演結(jié)果的交叉梯度數(shù)值為2.12×10-7，可見交叉梯度項的引入，使得同步反演結(jié)果的交叉梯度項數(shù)值更小.

對于電阻率和極化率結(jié)構(gòu)一致的模型來說，電阻率和極化率的同步反演結(jié)果無論從模型逼近程度還是從擬合差大小角度而言都優(yōu)于分步反演結(jié)果.

4.2 電阻率和極化率結(jié)構(gòu)不一致模型算例1

為了更好地討論同步反演相對于分步反演的優(yōu)越性，建立電阻率和極化率結(jié)構(gòu)不一致的模型，即電阻率異常體和極化體異常體的邊界并不完全一致.首先設(shè)計大尺寸低阻和小尺寸高極化的棱柱異常體模型，其中異常體電阻率為10 Ωm，尺寸為80 m×30 m，異常體極化率為12%，尺寸為40 m×30 m，二者頂面埋深均為40 m，埋藏于電阻率為100 Ωm，極化率為2%的均勻半空間中，具體可見圖2a和圖2d所示.觀測系統(tǒng)與前文一致，在正演數(shù)據(jù)中加入2%的高斯隨機誤差.對模型合成數(shù)據(jù)分別進行分步反演和同步反演，初始模型采用電阻率為100 Ωm，極化率為2%的均勻半空間，圖2b和圖2e是分步反演的結(jié)果，圖2c和圖2f是同步反演的結(jié)果.

對反演結(jié)果進行分析，結(jié)果顯示：電阻率的分步反演結(jié)果對棱柱體的空間形態(tài)有較好的恢復(fù)，異常體的埋深和幾何位置都較準確，但是棱柱體電阻率只能恢復(fù)到30 Ωm左右，與其電阻率真值（10 Ωm）有著一定的差距；極化率的分步反演結(jié)果能恢復(fù)出異常體的大致空間形態(tài)，異常體埋深較為準確，然而位置比真實情況偏左側(cè)，對棱柱體極化率真值（12%）最高只能恢復(fù)到9%，異常體兩側(cè)也存在著高極化假異常.相比分步反演，圖2c和圖2f的同步反演結(jié)果對棱柱體的空間形態(tài)有更好的恢復(fù)，棱柱體的電阻率反演結(jié)果更接近于電阻率真值，可達到18 Ωm，有一定程度的改善，并且電阻率異常體的尺寸也更接近于真實情況；而極化率的同步反演結(jié)果可達到11.5%左右，十分接近極化率真實數(shù)值，兩側(cè)高極化假異常得到一定程度的壓制，異常體的空間分布位置也更加準確，改善顯著.電阻率法分步反演迭代9次，數(shù)據(jù)擬合差從24.91下降到1.66，TDIP分步反演迭代10次，擬合差從17.11下降到2.13；同步反演迭代10次，視電阻率數(shù)據(jù)擬合差最終為1.55，視極化率數(shù)據(jù)擬合差為1.79，均略小于分步反演結(jié)果的數(shù)據(jù)擬合差.

圖1 物性結(jié)構(gòu)一致模型分步反演和同步反演結(jié)果對比圖

圖2 物性結(jié)構(gòu)不一致模型分步反演和同步反演結(jié)果對比圖

得到分步反演和同步反演的結(jié)果之后，分別計算二者電阻率模型和極化率模型交叉梯度數(shù)值，對于該組模型而言，分步反演結(jié)果的交叉梯度數(shù)值為2.04×10-4，同步反演結(jié)果的交叉梯度數(shù)值為2.47×10-5，可見交叉梯度項的引入，使得同步反演結(jié)果的交叉梯度項數(shù)值更小.

4.3 電阻率和極化率結(jié)構(gòu)不一致模型算例2

本文進一步設(shè)計了小尺寸高阻體和小尺寸高極化體的物性結(jié)構(gòu)不一致模型.模型異常體電阻率為1000 Ωm，分布在左側(cè)，極化率異常體的極化率為20%，分布在右側(cè)，二者尺寸均為40 m×30 m，頂面埋深均為40 m；異常體埋藏于電阻率為100 Ωm，極化率為2%的均勻半空間中，具體可見圖3a和圖3d所示.觀測系統(tǒng)與前文一致，在正演數(shù)據(jù)中加入2%的高斯隨機誤差.對模型合成數(shù)據(jù)分別進行分步反演和同步反演，初始模型采用電阻率為100 Ωm，極化率為2%的均勻半空間，分步反演的結(jié)果見圖3b和圖3e，同步反演的結(jié)果見圖3c和圖3f.

對反演結(jié)果進行分析，結(jié)果顯示：電阻率的分步反演結(jié)果對棱柱體的空間形態(tài)有較好的恢復(fù)，但是左側(cè)棱柱體電阻率最高恢復(fù)到160 Ωm左右，與電阻率真值（1000 Ωm）有著明顯差距，并且異常體的分布在左右兩側(cè)超出了實際范圍，尤其是右側(cè)明顯地超出實際邊界；極化率的分步反演結(jié)果能恢復(fù)出異常體的大致空間形態(tài)，但是反演結(jié)果埋深偏淺，對棱柱體極化率的真實數(shù)值（20%）最高只能恢復(fù)到10%左右.相比于分步反演，圖3c和圖3f中的同步反演結(jié)果對棱柱體的空間形態(tài)有更好的恢復(fù)，左側(cè)高阻棱柱體的反演電阻率可達到220 Ωm左右，數(shù)值上有較顯著的改善，并且高阻異常體反演結(jié)果的分布范圍更加準確，左右兩側(cè)偏出異常體實際分布范圍的情況得到了明顯改善，尤其是右側(cè)邊界基本與實際情況吻合，更接近于真實情況；而極化率的同步反演結(jié)果可以恢復(fù)到14.3%左右，相比于分步反演的結(jié)果改善顯著，并且反演結(jié)果埋深也更加準確，分步反演結(jié)果埋深偏淺的情況得到改善.電阻率法分步反演迭代8次，擬合差從6.47下降到1.74，TDIP分步反演迭代7次，擬合差從25.66下降到2.06，同步反演迭代了9次，最后收斂時視電阻率擬合差為1.43，視極化率數(shù)據(jù)擬合差為1.55，均小于分步反演結(jié)果的數(shù)據(jù)擬合差.

得到分步反演和同步反演的結(jié)果之后，計算電阻率模型和極化率模型交叉梯度數(shù)值，對于該組模型而言，分步反演結(jié)果的交叉梯度數(shù)值為1.27×10-4，同步反演結(jié)果的交叉梯度數(shù)值為1.19×10-6，可見交叉梯度項的引入，使得同步反演結(jié)果的交叉梯度項數(shù)值更小.

4.4 加入交叉梯度項與未加入交叉梯度項對比模型

為了更好地突出目標函數(shù)中的交叉梯度項在反演中所起到的作用，本文設(shè)計了一組加入交叉梯度項和未加入交叉梯度項的對比模型.模型異常體電阻率為1000 Ωm，極化率異常體的極化率為20%，異常體尺寸為80 m×30 m，頂面埋深為40 m；異常體埋藏于電阻率為100 Ωm，極化率為2%的均勻半空間中.觀測系統(tǒng)與前文一致，在正演數(shù)據(jù)中加入1.5%的高斯隨機誤差.對模型合成數(shù)據(jù)分別進行分步反演、未加入交叉梯度項的同步反演和含交叉梯度項的同步反演，初始模型均采用電阻率為100 Ωm，極化率為2%的均勻半空間，分步反演的結(jié)果見圖4a和圖4d，不含交叉梯度項同步反演的結(jié)果見圖4b和圖4e，含交叉梯度項同步反演結(jié)果如圖4c和圖4f所示.

對反演結(jié)果進行分析，結(jié)果顯示：電阻率的分步反演結(jié)果對棱柱體的空間形態(tài)有一定程度的恢復(fù)，電阻率最高恢復(fù)到180 Ωm左右，與棱柱體電阻率真值（1000 Ωm）有著明顯差距，并且異常體左右兩側(cè)存在一定程度的假異常；極化率的分步反演結(jié)果能恢復(fù)出異常體的大致空間形態(tài)，但是反演結(jié)果對棱柱體極化率的真實數(shù)值（20%）最高只能恢復(fù)到8%左右.相比于分步反演，圖4b和圖4e中的未加入交叉梯度項同步反演結(jié)果對棱柱體的空間形態(tài)有更好的恢復(fù)，高阻棱柱體的電阻率可達到240 Ωm左右，數(shù)值上有較顯著的改善，但是異常體還是存在一定程度的假異常；而極化率未加入交叉梯度項同步反演結(jié)果可以恢復(fù)到12%左右，相比于分步反演的結(jié)果改善顯著，但是反演結(jié)果分布與實際有一定的差距，左右兩側(cè)存在較大范圍的高極化假異常.圖4c和圖4f是含交叉梯度項同步反演結(jié)果，其中反演得到的電阻率模型更加接近于真實情況，可將電阻率數(shù)值恢復(fù)到280 Ωm左右，并且異常體兩側(cè)的假異常明顯減少；而極化率結(jié)果也得到了改善，數(shù)值可以達到14.5%，兩側(cè)的假異常也明顯減少，因為加入了交叉梯度項，極化率和電阻率模型的相似性得以提高，因此異常體的極化率分布更加準確.電阻率法分步反演迭代7次，擬合差從7.42下降到1.84，TDIP分步反演迭代7次，擬合差從21.81下降到1.96，未加入交叉梯度項的同步反演迭代了9次，最后收斂時視電阻率擬合差為1.63，視極化率數(shù)據(jù)擬合差為1.75，加入了交叉梯度項的同步反演迭代了11次，最后收斂時視電阻率擬合差為1.55，視極化率數(shù)據(jù)擬合差為1.45，是三種反演方法中擬合差最小的.

圖3 物性結(jié)構(gòu)不一致模型分步反演和同步反演結(jié)果對比圖

圖4 無交叉梯度項和含有交叉梯度項同步反演結(jié)果對比圖

得到不含交叉梯度項反演和含交叉梯度項反演的結(jié)果之后，計算電阻率模型和極化率模型交叉梯度數(shù)值，對于該組模型而言，前者的交叉梯度數(shù)值為2.37×10-5，后者的交叉梯度數(shù)值為1.84×10-7，可見交叉梯度項的引入，使得同步反演結(jié)果的交叉梯度項數(shù)值更小.

總之，加入交叉梯度項的電阻率和極化率的同步反演結(jié)果在模型逼近程度和數(shù)據(jù)擬合差大小方面均優(yōu)于分步反演和不含交叉梯度項的同步反演結(jié)果.

5 結(jié)論

本文實現(xiàn)了基于交叉梯度耦合的TDIP二維同步反演算法，并設(shè)計了電阻率和極化率結(jié)構(gòu)一致和不一致的模型對反演算法進行了理論模型測試.通過對比合成數(shù)據(jù)的同步反演和分步反演結(jié)果，可以得出以下結(jié)論：

（1） 無論是電阻率和極化率結(jié)構(gòu)一致還是不一致的模型，電阻率的同步反演結(jié)果比分步反演對異常體的空間形態(tài)都有更好的恢復(fù)，尤其是高阻異常體改善效果更加明顯.對于低阻異常體而言，電阻率的同步反演結(jié)果更接近于電阻率真值，并且異常體分布范圍更加準確，而高阻體電阻率的同步反演結(jié)果不僅更加接近于電阻率的實際數(shù)值，其埋深和空間分布也更加準確.

（2） 對于極化率模型而言，同步反演能改善極化率反演結(jié)果，使其更接近甚至完全達到極化率真值，并且極化率分步反演結(jié)果埋深和空間分布不準確的問題也得到改善.

（3） 兩組電阻率和極化率結(jié)構(gòu)不一致的模型算例表明，同步反演在電阻率和極化率模型邊界不完全一致的時候也能改善反演結(jié)果，尤其是電阻率和極化率模型的空間分布準確性改善顯著.

（4） 對于不同模型而言，同步反演結(jié)果不但在模型逼近程度上優(yōu)于分步反演的結(jié)果，同步反演結(jié)果的數(shù)據(jù)擬合差也要小于分步反演的數(shù)據(jù)擬合差.

（5） 加入了交叉梯度約束，使得同步反演算法結(jié)果交叉梯度項數(shù)值更小，更加符合真實情況，因而該算法具有更普遍的適用性.

致謝 兩位匿名評審專家為本文修改提出了細致而寶貴的意見，使得文章增色不菲，在此表示感謝．

Aric K, Adam A, Smythe D K. 1997. Combined seismic and magnetotelluric imaging of upper crystalline crust in the Southern Bohemian Massif.FirstBreak, 15（8）: 265-271. Colombo D, De Stefano M. 2007. Geophysical modeling via simultaneous joint inversion of seismic, gravity, and electromagnetic data: Application to prestack depth imaging.TheLeadingEdge, 26（3）: 326-331. Constable S C, Parker R L, Constable C C. 1987. Occam′s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data.Geophysics, 52（3）: 289-300. Douglas W O, Li Y G. 1994. Inversion of induced polarization data.Geophysics, 59（9）: 1327-1341.

Douglas W O, Li Y G. 1999. Estimating depth of investigation in DC resistivity and IP surveys.Geophysics, 64（2）: 403-416.

Fu L K. 1982. Induced Polarization Method （in Chinese）. Beijing: Geological Publishing House.

Gallardo L A, Meju M A. 2003. Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismic data.Geophys.Res.Lett., 30（13）: 1658.

Gallardo L A, Meju M A. 2007. Joint two-dimensional cross-gradient imaging of magnetotelluric and seismic traveltime data for structural and lithological classification.Geophys.J.Int., 169（3）: 1261-1272.

Gao D Z, Zhao J H, Bo Y L, et al. 2004. A profile study of gravitative-magnetic and seismic comprehensive survey in the East China Sea.ChineseJ.Geophys. （in Chinese）, 47（5）: 853-861, doi: 10.3321/j.issn:0001-5733.2004.05.017.Haber E, Oldenburg D. 1997. Joint inversion: A structural approach.InverseProblems, 13（1）: 63-77.

Hu P, Liu B J, Bai L X, et al. 2010. Synthetic exploration of the buried faults in Olympic park area.ChineseJ.Geophys. （in Chinese）, 53（6）: 1486-1494, doi: 10.3969/j.issn.0001-5733.2010.06.026.

Jiang M, Peng M, Wang Y X, et al. 2012a. Geophysical evidence for deep subduction of Indian lithospheric plate beneath Eastern Himalayan Syntaxis.ActaPetrologicaSinica（in Chinese）, 28（6）: 1755-1764.

Jiang M, Tan H D, Qian H, et al. 2012b. Geophysical deep structure and genetic model of Jinchuan copper-nickel deposit.MineralDeposits（in Chinese）, 31（2）: 207-215.

Moorkamp M, Jones A G, Eaton D W. 2007. Joint inversion of teleseismic receiver functions and magnetotelluric data using a genetic algorithm: Are seismic velocities and electrical conductivities compatible?Geophys.Res.Lett., 34（16）: L16311. Peng M, Tan H D, Jiang M, et al. 2012. Joint inversion of receiver functions and magnetotelluric data: Application to crustal and mantle structure beneath central Namche Barwa, eastern Himalayan syntaxis.ChineseJ.Geophys. （in Chinese）, 55（7）: 2281-2291, doi: 10.6038/j.issn.0001-5733.2012.07.014. Peng M, Tan H D, Jiang M, et al. 2013. Three-dimensional joint inversion of magnetotelluric and seismic travel time data with cross-gradient constraints.ChineseJ.Geophys. （in Chinese）, 56（8）: 2728-2738, doi: 10.6038/cjg20130821.

Siripunvaraporn W, Egbert G. 2000. An efficient data subspace inversion method for 2-D magnetotelluric data.Geophysics, 65（3）: 791-803.

Siripunvaraporn W, Egbert G, Lenbury Y, et al. 2005. Three-dimensional magnetotelluric inversion: data-space method.Phys.EarthPlanet.Interiors, 150（1-3）: 3-14.

Yu P, Dai M G, Wang J L, et al. 2009. Joint inversion of magnetotelluric and seismic data based on random resistivity and . （in Chinese）, 52（4）: 1089-1097, doi: 10.3969/j.issn.0001-5733.2009.04.026.

附中文參考文獻

傅良魁. 1982. 激發(fā)極化法. 北京: 地質(zhì)出版社. 高德章, 趙金海, 薄玉玲等. 2004. 東海重磁地震綜合探測剖面研究. 地球物理學(xué)報, 47（5）: 853-861, doi: 10.3321/j.issn:0001-5733.2004.05.017. 胡平, 劉保金, 白立新等. 2010. 奧林匹克公園地區(qū)隱伏斷裂綜合

探測. 地球物理學(xué)報, 53（6）: 1486-1494, doi: 10.3969/j.issn.0001-5733.2010.06.026.

姜枚, 彭淼, 王有學(xué)等. 2012a. 喜馬拉雅東構(gòu)造結(jié)巖石圈板片深俯沖的地球物理證據(jù). 巖石學(xué)報, 28（6）: 1755-1764.

姜枚, 譚捍東, 錢輝等. 2012b. 金川銅鎳礦床的地球物理深部結(jié)構(gòu)與成因模式. 礦床地質(zhì), 31（2）: 207-215.

彭淼, 譚捍東, 姜枚等. 2012. 利用接收函數(shù)和大地電磁數(shù)據(jù)聯(lián)合反演南迦巴瓦構(gòu)造結(jié)中部地區(qū)殼幔結(jié)構(gòu). 地球物理學(xué)報, 55（7）: 2281-2291, doi: 10.6038/j.issn.0001-5733.2012.07.014. 彭淼, 譚捍東, 姜枚等. 2013. 基于交叉梯度耦合的大地電磁與地震走時資料三維聯(lián)合反演. 地球物理學(xué)報, 56（8）: 2728-2738, doi: 10.6038/cjg20130821.

于鵬, 戴明剛, 王家林等. 2009. 電阻率和速度隨機分布的MT與地震聯(lián)合反演. 地球物理學(xué)報, 52（4）: 1089-1097, doi: 10.3969/j.issn.0001-5733.2009.04.026.

（本文編輯 何燕）

Two-dimensional synchronous inversion of TDIP with cross-gradient constraint

LI Zhao-Xiang1,2, TAN Han-Dong1*, FU Shao-Shuai1, TANG Fei1, ZHU De-Peng1, GAO Jing-Yu1

1SchoolofGeophysicsandInformationTechnology,ChinaUniversityofGeosciences,Beijing100083,China2GeologicalExplorationTechnologyInstituteofJiangsuProvince,Nanjing210000,China

The existing inversion algorithm of time-domain induced polarization（TDIP）data is a step-wise one, in which the polarizability value is gained after the calculation of the resistivity value is completed with the resistivity value fixed, which results in that the polarizability value depends on the resistivity value．To solve this problem, we study a two-dimensional synchronous inversion algorithm of TDIP data with cross-gradient constraint．Based on well-developed two-dimensional forward and inversion algorithms for DC resistivity，we propose two-dimensional forward and inversion algorithms for induced polarization and synchronous inversion algorithm with common inversion grids and cross-gradient as the structural constraint. Both resistivity and polarizability models are able to be synchronously updated．Trial computation of synchronous inversion has been carried out on the synthetic data from a model with the same physical structure and different physical structures. Results indicate that compared with step-wise inversion，synchronous inversion is better in recovering both the various models′ numerical values and spatial morphology．The resistivity and polarizability values of the conductive-high polarizability models are much closer to the true model in synchronous inversion，while by the resistive-high polarizability model，synchronous inversion not only can improve the values of TDIP data to some extent，but also has large improvement in increasing the accuracy of spatial distribution．Moreover，synchronous inversion can improve the result of the model with different physical structures and thus has wide applicability．

Time-domain induced polarization method（TDIP）； Data-space inversion； Cross-gradient； Synchronous inversion

國家自然科學(xué)基金項目（41374078）和國土資源部地質(zhì)調(diào)查項目（12120113086100，12120113101300）聯(lián)合資助.

李兆祥，男，1989年生，碩士研究生，主要從事地球物理地電學(xué)算法研究. E-mail：lizhaoxiang1111@163.com

*通訊作者 譚捍東，男，1966年生，教授，主要從事電法勘探理論及應(yīng)用研究. E-mail：thd@cugb.edu.cn

10.6038/cjg20151232.

10.6038/cjg20151232

P631

2014-12-22，2015-11-09收修定稿

李兆祥， 譚捍東， 付少帥等. 2015. 基于交叉梯度約束的時間域激發(fā)極化法二維同步反演.地球物理學(xué)報，58（12）:4718-4726,

Li Z X, Tan H D, Fu S S, et al. 2015. Two-dimensional synchronous inversion of TDIP with cross-gradient constraint.ChineseJ.Geophys. （in Chinese）,58（12）:4718-4726,doi:10.6038/cjg20151232.