葉益信，李澤林，付 宸，丁尚見

（1.東華理工大學(xué) 放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室，南昌 330013；2.中國地質(zhì)大學(xué) 地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室，武漢 430074）

利用阻尼型高斯牛頓法的激發(fā)極化數(shù)據(jù)聚焦反演

葉益信1，2，李澤林1，付 宸1，丁尚見1

（1.東華理工大學(xué) 放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室，南昌 330013；2.中國地質(zhì)大學(xué) 地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室，武漢 430074）

利用阻尼型高斯牛頓法求解反演目標(biāo)函數(shù)，并將粗糙度矩陣約束和最小支撐泛函約束以及參考模型和模型參數(shù)界限約束引入到目標(biāo)函數(shù)中，同時(shí)結(jié)合預(yù)條件共軛梯度法進(jìn)行反演迭代，實(shí)現(xiàn)了激發(fā)極化數(shù)據(jù)三維聚焦反演快速計(jì)算。通過對(duì)幾個(gè)典型模型的試算，并與傳統(tǒng)光滑模型約束的反演結(jié)果進(jìn)行比較，表明該反演算法的可靠性和穩(wěn)定性較好，并且算法速度快，占用內(nèi)存低且易于并行化。

激發(fā)極化；高斯牛頓法；聚焦反演；最小支撐泛函

0 引言

激發(fā)極化（IP）法的應(yīng)用范圍非常廣泛，無論是在金屬礦和非金屬礦產(chǎn)勘查［1－2］，還是在尋找地下水資源和地?zé)崽锓矫妫?－4］，都獲得了成功地應(yīng)用，近年來又拓展到環(huán)境監(jiān)測(cè)領(lǐng)域［5－6］。國內(nèi)、外有關(guān)學(xué)者對(duì)激發(fā)極化法正反演也有較多的研究，Pelton等［7］利用參考數(shù)據(jù)庫中數(shù)據(jù)插值計(jì)算偏導(dǎo)數(shù)，實(shí)現(xiàn)了激發(fā)極化二維最小二乘反演方法；日本的Sasaki博士［8］在前人的基礎(chǔ)上發(fā)展了有限元電阻率／極化率數(shù)據(jù)的二維最小二乘反演方法；Ellis等［9］也將解的先驗(yàn)信息加入目標(biāo)函數(shù)中，使反演結(jié)果更符合實(shí)際；Oldenburg等［10］在電阻率反演的基礎(chǔ)上，討論了幾種激發(fā)極化反演方法；Li等［11］研究了激發(fā)極化數(shù)據(jù)的三維反演；阮百堯等［12］在前人研究的基礎(chǔ)上，實(shí)現(xiàn)了電阻率和極化率分塊連續(xù)變化的激發(fā)極化數(shù)據(jù)二維最小二乘反演算法，并編制了相應(yīng)的程序，在各種金屬礦和非金屬礦勘探的數(shù)據(jù)處理中起了非常大的作用；吳小平［13］利用共軛梯度法實(shí)現(xiàn)了激發(fā)極化數(shù)據(jù)的三維快速反演，大大提高了反演速度；黃俊革［14］利用異常電位法實(shí)現(xiàn)了電阻率和激發(fā)極化法三維有限元正反演，提高了正演和反演精度；陳進(jìn)超等［15］應(yīng)用非結(jié)構(gòu)化三角網(wǎng)格剖分實(shí)現(xiàn)了激發(fā)極化法2.5維有限元正演，模擬具有較高的計(jì)算精度和較快的計(jì)算速度；蔣首進(jìn)［16］研究了激發(fā)極化法二維最小二乘反演方法并將其應(yīng)用到礦產(chǎn)勘查中，取得了較好的應(yīng)用效果。

盡管對(duì)激發(fā)極化數(shù)據(jù)的正反演研究較多，但還需從不同側(cè)面進(jìn)行研究，以期提高激發(fā)極化法的應(yīng)用效果。這里對(duì)如何提高激發(fā)極化數(shù)據(jù)反演分辨率問題，給出了具體解決方案。粗糙度約束有效排除了不真實(shí)的多余結(jié)構(gòu)，引入最小支撐泛函約束［17］即聚焦突出異常體邊界，異常體范圍小，反演結(jié)果更趨于真值。同時(shí)利用阻尼型高斯牛頓法求解反演目標(biāo)函數(shù)，并且采用共軛梯度法進(jìn)行反演迭代，實(shí)現(xiàn)了激發(fā)極化數(shù)據(jù)三維快速反演計(jì)算，解決了內(nèi)存占用量巨大的問題，且易于并行化。通過對(duì)典型模型的反演試算，表明了該反演算法的可靠性和穩(wěn)定性較好。

1 反演理論

首先給出一個(gè)基本的正則化反演目標(biāo)函數(shù)：

其中：Wd、Wm分別表示數(shù)據(jù)和模型的權(quán)系數(shù)矩陣；dobs為觀測(cè)數(shù)據(jù)；d（m）為正演理論值；m為模型參數(shù)；mref為參考模型；α是正則化因子。

在本算法中，Wd采用式（2）的對(duì)角矩陣：

式中：S（dobs）為觀測(cè)數(shù)據(jù)的標(biāo)準(zhǔn)偏差；ε為權(quán)重項(xiàng)，避免在很小的數(shù)據(jù)上權(quán)重過大。

為了減小多解性，使反演結(jié)果更加真實(shí)，光滑反演通常采用一階正則化或二階正則化粗糙度矩陣作為模型的權(quán)系數(shù)矩陣，一階正則化粗糙度矩陣為式（3）：

式中：Rx、Ry和Rz分別為x、y和z方向的梯度。

雖然采用一階正則化粗糙度約束能使反演結(jié)果呈現(xiàn)良好的平緩結(jié)構(gòu)，有效排除不真實(shí)的多余突變結(jié)構(gòu)，但地質(zhì)體之間通常存在明顯的分界面，不是平緩過渡的。因此光滑反演結(jié)果通常過于平緩，故不利于地質(zhì)解釋。根據(jù)聚焦反演原理［17］，引入最小支撐泛函作為模型約束條件（式（4））。

其中β為不等于“0”的小數(shù)，根據(jù)最小支撐泛函的性質(zhì)，在反演過程中，模型參數(shù)中異常體剖面的面積會(huì)盡可能聚焦到最小，從而可以提高模型結(jié)構(gòu)的分辨率。然后將模型參數(shù)轉(zhuǎn)化為加權(quán)模型空間形式，令，由此新的反演目標(biāo)函數(shù)可改寫為式（5）：

有了目標(biāo)函數(shù)（式（5））后，用高斯牛頓法選擇合適的正則化參數(shù)來反演出符合觀測(cè)數(shù)據(jù)的模型，首先定義一個(gè)初始模型mwi，通常采用最接近真實(shí)模型的均勻半空間，然后對(duì)式（5）線性化后得到式（6）。

為了得到式（6）的最小值，將式（6）對(duì)mw求導(dǎo)，并使其等于零，得到高斯牛頓更新公式為式（7）。

求解式（7）可得到模型的修正量δmw，由此可得到一個(gè)新的模型（式（8））。

式中：α是線性搜索參數(shù)，這里采用Armijo步長準(zhǔn)則對(duì)α線性搜索，確保目標(biāo)函數(shù)能夠充分減小。

迭代后，未加權(quán)的模型參數(shù)可由式（9）求得。

重復(fù)上述步驟，直到目標(biāo)函數(shù)減小到一定值或模型修正量小于一定值。

同時(shí)為了縮小解的范圍，減少解的非唯一性，提高反演分辨率，對(duì)模型參數(shù)施加界限約束，假設(shè)在反演中地形介質(zhì)電性變化的上邊界η＋（r）和下邊界η－（r），這樣，在每次迭代過程中可將參數(shù)的變化量控制在界限范圍內(nèi)，而且也不會(huì)出現(xiàn)負(fù)的極化率值，算法可寫為式（10）。

這樣參數(shù)范圍被限定在：［η－（r），η＋（r）］，背景均勻參數(shù)和最大異常參數(shù)可由先驗(yàn)信息獲得。

2 模型反演試驗(yàn)

反演中的觀測(cè)數(shù)據(jù)為E－SCAN方式測(cè)量的雙極－雙極電位值φobs。由于它們的變化范圍很大，一般用對(duì)數(shù)來標(biāo)定觀測(cè)數(shù)據(jù)及模型參數(shù)，即dobs＝lnφobs及m＝lnσ。視極化率數(shù)據(jù)可依據(jù)Seigel［18］理論得到，可表示為式（11）。

寫成矩陣形式為式（12）。

其中J為電阻率反演中的偏導(dǎo)數(shù)矩陣的負(fù)矩陣。由于J已在電阻率反演中得到，因此在電阻率三維反演基礎(chǔ)上，只需少量計(jì)算即可獲得激發(fā)極化數(shù)據(jù)的反演結(jié)果。電極系布置如圖1所示，共4 332個(gè)數(shù)據(jù)。

圖1 地表電極排列布置示意圖Fig.1 The configuration of the electrodes array on the earth surface

2.1 單個(gè)異常體模型

模型一為一個(gè)大小為20m×20m×17m、電阻率為200Ω·m、極化率為1%的均勻半空間含有一個(gè)大小為2.5m×2.5m×4m、電阻率為10Ω·m、極化率為10%的異常體，其頂部埋深為2m，其z＝4m處的水平切片（XOY平面）如圖2（a）所示，y＝0 m處的垂直切片（XOZ平面）如圖2（b）所示。模型網(wǎng)格剖分成40×40×20共32 000個(gè)單元。

反演初始模型為電阻率200Ω·m、極化率1%的均勻半空間，所有計(jì)算均在雙核2.9GHz CPU 和8GRAM個(gè)人計(jì)算機(jī)上完成。共迭代10次，共耗時(shí)13.5min。電阻率和極化率反演結(jié)果分別如圖3和圖4所示，并與光滑模型約束的反演結(jié)果進(jìn)行了比較，黑線框表示模型的邊界位置。對(duì)比兩種模式的反演結(jié)果可看出，光滑模型約束的反演結(jié)果在異常體位置出現(xiàn)了低阻高極化異常值，但低阻異常值約80Ω·m左右，極化率異常值為7%左右，基本能夠表征異常體的存在，但是異常范圍偏大，尤其是電阻率的反演結(jié)果，與異常體真實(shí)電阻率值相差較大，而聚焦反演結(jié)果電阻率極化率異常更明顯，電阻率異常最小值接近10Ω·m，極化率異常最大幅值接近10%，與真實(shí)異常體特征更接近，改善了異常體的聚焦效果，對(duì)異常體的定位更準(zhǔn)確。

圖5為兩種反演方法的歸一化迭代擬合差隨迭代次數(shù)的變化情況，從圖5中可看出，兩種反演算法都能穩(wěn)定收斂，但是聚焦反演的擬合差更小，收斂速度較快。

2.2 兩個(gè)異常體模型

模型二為一個(gè)大小為20m×20m×17m、電阻率為200Ω·m、極化率為1%的均勻半空間含有兩個(gè)大小為2.5m×2.5m×4m、電阻率為10 Ω·m、極化率為10%的異常體，其頂部埋深為2m，其z＝4m處的水平切片（XOY平面）如圖6（a）所示，y＝0m處的垂直切片（XOZ平面）如圖6（b）所示。網(wǎng)格剖分與模型一相同。

圖2 模型一切片圖Fig.2 A depth slice of the model 1

圖3 模型一電阻率光滑反演與聚焦反演結(jié)果比較Fig.3 Comparisons of focusing inversion with smoothing inversion for model 1

圖4 模型一極化率光滑反演與聚焦反演結(jié)果比較Fig.4 Comparisons of focusing inversion with smoothing inversion for model 1

反演初始模型為電阻率200Ω·m、極化率1%的均勻半空間，反演共迭代10次，共耗時(shí)13.6 min。圖7和圖8分別為電阻率和極化率的反演結(jié)果，并與光滑模型約束的反演結(jié)果進(jìn)行了比較，黑線框表示模型的邊界位置。對(duì)比兩種模式的反演結(jié)果可看出，光滑模型約束的反演結(jié)果在異常體位置出現(xiàn)了低阻高極化異常，低阻異常值約90Ω·m左右，極化率異常值為6%左右，基本能夠表征異常體的存在，但是整體聚焦效果較差，圖像較模糊，尤其是電阻率反演結(jié)果，兩個(gè)異常體的電阻率值與真實(shí)值相差較大，而聚焦反演結(jié)果低阻高極化異常更明顯，電阻率異常最小值接近20Ω·m，極化率異常最大值接近10%，與真實(shí)異常體特征更接近，異常體的聚焦效果明顯提高。

圖5 模型一反演歸一化擬合差曲線Fig.5 The normalized misfit curves with iteration inversions for model 1

3 結(jié)論

作者將最小支撐泛函引入到激發(fā)極化數(shù)據(jù)三維反演的目標(biāo)函數(shù)中，然后采用阻尼型高斯牛頓法求解反演目標(biāo)函數(shù)最優(yōu)化問題，同時(shí)結(jié)合預(yù)條件共軛梯度法得到激發(fā)極化數(shù)據(jù)三維聚焦反演結(jié)果。通過兩個(gè)模型的試算，算法反演結(jié)果增強(qiáng)了異常體成像的聚焦效果，也更好地框定了異常體的范圍，反演結(jié)果更接近于真實(shí)模型。相比較而言，電阻率的反演聚焦效果更明顯，極化率的反演無論是光滑反演還是聚焦反演，對(duì)異常體的分辨均較好，但對(duì)模型的邊界位置，尤其是下底邊界擬合較差。同時(shí)該算法具有計(jì)算速度快，占用內(nèi)存低且易于并行化等優(yōu)點(diǎn)。

圖6 模型二切片圖Fig.6 A depth slice of the model 2

圖7 模型二電阻率光滑反演與聚焦反演結(jié)果比較Fig.7 Comparisons of focusing inversion with smoothing inversion for model 2

圖8 模型二極化率光滑反演與聚焦反演結(jié)果比較Fig.8 Comparisons of focusing inversion with smoothing inversion for model 2

［1］ 魏玉山，朱春生，李凱春.激發(fā)極化法在金星金礦勘查的應(yīng)用效果［J］.吉林地質(zhì)，2010，29（2）：89－93.WEI Y S，ZHU CH S，LIU K CH.The application of induce polarization method in Jinxing gold deposit［J］.Jilin Geology，2010，29（2）：89－93.（In Chinese）

［2］ 王平康，冉波，龍鋒，等.激發(fā)極化法在內(nèi)蒙某鉛鋅礦點(diǎn)成礦預(yù)測(cè)中的應(yīng)用［J］.中國礦業(yè)，2010，19（7）：108－110.WANG P K，LAN B，LONG F，et al.Application of induce polarization in a lead－zinc deposit metallogenetic prognostication，Inner Mongolia［J］.China Mining Magazine，2010，19（7）：108－110.（In Chinese）

［3］ 陳青松，孫啟斌，郭成有.激發(fā)極化法在新疆地區(qū)找水的應(yīng)用［J］.西部探礦工程，2002，79（6）：55－57.CHEN Q S，SUN Q B，GUO CH Y.The application of induce polarization in finding water in XinJiang province［J］.West－China Exploration Engineering，2002，79（6）：55－57.（In Chinese）

［4］ 杜炳銳，傅順，楊威.激發(fā)極化法在磨盤山地?zé)峥辈橹械膽?yīng)用［J］.物探化探計(jì)算技術(shù)，2010，32（5）：514 －517.DU B R，F(xiàn)U S，YANG W.Application of induce polarization in Mopanshan geothermal resources exploration［J］.Computing Techniques of Geophysical and Geochemical Exploration，2010，32（5）：514－517.（In Chinese）

［5］ BARKER R D.Investigation of groundwater salinity by geophysical methods，in Ward，S.H.，Ed.，Geotechnical and environmental geophysics［J］.Investigations in Geophysics，1990，2（5）：201－211.

［6］ SLATER L D，LESMES D.IP interpretation in environmental investigations［J］.Geophysics，2002，67 （1）：77－88.

［7］ PELTON W H，RIJO L，SWIFT C M.Inversion of two－dimensional resistivity and induced polarization data［J］.Geophysics，1978，43（4）：788－803.

［8］ SASAKI Y.Automatic interpretation of induced polarization data over two－dimensional structures［J］.Memoirs of the Faculty of Engineering，Kyushu University，1982，42（1）：59－74.

［9］ ELLIS R G，OLDENBURG D W.Applied geophysical inversion［J］.Geophysical Journal International，1994，116：5－10.

［10］OLDENBURG D W，LI Y.Inversion of induced po－larization data［J］.Geophysics，1994，59（9）：1327－1341.

［11］LI Y，OLDENBURG D W.3－D Inversion of induced polarization data［J］.Geophysics，2000，65（6）：1931 －1945.

［12］阮百堯，村上裕，徐世浙.電阻率／激發(fā)極化率數(shù)據(jù)的二維反演程序［J］.物探化探計(jì)算技術(shù)，1999，21（2）：116－125.RUAN B Y，YUTAKA M，XU SH ZH.2Dinversion programs of induced polarization data［J］.Computing Techniques of Geophysical and Geochemical Exploration，1999，21（2）：116－125.（In Chinese）

［13］吳小平.利用共軛梯度方法的激發(fā)極化三維快速反演［J］.煤田地質(zhì)與勘探，2004，32（5）：62－64.WU X P.Rapid 3Dinversion of induced polarization data using conjugate gradient method［J］.Coal Geology ＆Exploration，2004，32（5）：62－64.（In Chinese）

［14］黃俊革.三維電阻率／極化率有限元正演模擬與反演成像［D］.長沙：中南大學(xué)，2003.HUANG J G.3Dresistivity／IP modeling and inversion based on FEM［D］.Changsha：Central South University，2003.（In Chinese）

［15］陳進(jìn)超，王緒本，王麗坤.時(shí)間域激發(fā)極化法非結(jié)構(gòu)化三角網(wǎng)格有限元正演模擬［J］.物探化探計(jì)算技術(shù)，2011，33（4）：411－417.CHEN J C，WANG X B，WANG L K.Unstructured triangular grid finite element method for time domain induce polarization forward modeling［J］.Computing Techniques of Geophysical and Geochemical Exploration，2011，33（4）：411－417.（In Chinese）

［16］蔣首進(jìn).激發(fā)極化二維正反演研究及其在礦產(chǎn)中的應(yīng)用［D］.成都：成都理工大學(xué)，2013.JIANG SH J.To excitation polarization numerical simulation and its application in mineral［D］.Chengdu：Chengdu University of Technology，2013.（In Chinese）

［17］ZHDANOV M S，ELLIS R，MUKHERJEE S.Three －D regularized focusing inversion of gravity gradient tensor data［J］.Geophysics，2004，69：925－937.

［18］SEIGEL H O.Mathematical formulation and type curves for induced polarization［J］.Geophysics，1959，24：547－565.

Regularized focusing inversion of induced polarization data using damped Gauss－Newton method

YE Yi－xin1，2，LI Ze－lin1，F(xiàn)U Chen1，DING Shang－jian1

（1.Fundamental Science on Radioactive Geology and Exploration Technology Laboratory，East China Institute of Technology，Nanchang 330013，China；2.Hubei Subsurface Multi－scale Imaging Lab（SMIL），China University of Geosciences（Wuhan），Wuhan 430074，China）

In this paper，roughness，minimum support functional constraints and reference model，model parameter range constraints are incorporated into the objective function of induced polarization（IP）data inversion.Then we use the damped Gauss－Newton method to find the optimization of the objective function，with the model update being calculated using apreconditioned conjugate gradient algorithm.The method is demonstrated on two synthetic models by comparisons with the results obtained by the traditional inversion method based on smoothness stabilizing functional constraint.The model tests show that the proposed inversion algorithm has good reliability，stability，and high speed.The comparisons show that the program is more convergent and this approach helps to generate more focused and resolved images of blocky geoelectrical structures than traditional inversion approach.

induce polarization（IP）；Gauss－Newton method；focusing inversion；minimum support functional

P 631.3

：A

10.3969／j.issn.1001－1749.2015.06.02

1001－1749（2015）06－0680－07

2014－11－13改回日期：2015－01－07

國家自然科學(xué)基金（41204055）；中國地質(zhì)大學(xué)（武漢）地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室開放基金項(xiàng)目（SMIL－2014－06）

葉益信（1983－），男，博士，主要從事電法、電磁法正反演研究和應(yīng)用研究工作，E－mail：yixinye321＠126.com。