利用阻尼型高斯牛頓法的激發極化數據聚焦反演

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

（1.東華理工大學 放射性地質與勘探技術國防重點學科實驗室，南昌 330013；2.中國地質大學 地球內部多尺度成像湖北省重點實驗室，武漢 430074）

利用阻尼型高斯牛頓法的激發極化數據聚焦反演

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

（1.東華理工大學 放射性地質與勘探技術國防重點學科實驗室，南昌 330013；2.中國地質大學 地球內部多尺度成像湖北省重點實驗室，武漢 430074）

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

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

0 引言

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

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

1 反演理論

首先給出一個基本的正則化反演目標函數：

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

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

式中：S（dobs）為觀測數據的標準偏差；ε為權重項，避免在很小的數據上權重過大。

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

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

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

其中β為不等于“0”的小數，根據最小支撐泛函的性質，在反演過程中，模型參數中異常體剖面的面積會盡可能聚焦到最小，從而可以提高模型結構的分辨率。然后將模型參數轉化為加權模型空間形式，令，由此新的反演目標函數可改寫為式（5）：

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

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

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

式中：α是線性搜索參數，這里采用Armijo步長準則對α線性搜索，確保目標函數能夠充分減小。

迭代后，未加權的模型參數可由式（9）求得。

重復上述步驟，直到目標函數減小到一定值或模型修正量小于一定值。

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

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

2 模型反演試驗

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

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

其中J為電阻率反演中的偏導數矩陣的負矩陣。由于J已在電阻率反演中得到，因此在電阻率三維反演基礎上，只需少量計算即可獲得激發極化數據的反演結果。電極系布置如圖1所示，共4 332個數據。

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

2.1 單個異常體模型

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

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

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

2.2 兩個異常體模型

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

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

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

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

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

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

3 結論

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

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

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

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

［1］ 魏玉山，朱春生，李凱春.激發極化法在金星金礦勘查的應用效果［J］.吉林地質，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］ 王平康，冉波，龍鋒，等.激發極化法在內蒙某鉛鋅礦點成礦預測中的應用［J］.中國礦業，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］ 陳青松，孫啟斌，郭成有.激發極化法在新疆地區找水的應用［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］ 杜炳銳，傅順，楊威.激發極化法在磨盤山地熱勘查中的應用［J］.物探化探計算技術，2010，32（5）：514 －517.DU B R，FU 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］阮百堯，村上裕，徐世浙.電阻率／激發極化率數據的二維反演程序［J］.物探化探計算技術，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］吳小平.利用共軛梯度方法的激發極化三維快速反演［J］.煤田地質與勘探，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］.長沙：中南大學，2003.HUANG J G.3Dresistivity／IP modeling and inversion based on FEM［D］.Changsha：Central South University，2003.（In Chinese）

［15］陳進超，王緒本，王麗坤.時間域激發極化法非結構化三角網格有限元正演模擬［J］.物探化探計算技術，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］蔣首進.激發極化二維正反演研究及其在礦產中的應用［D］.成都：成都理工大學，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，FU 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

國家自然科學基金（41204055）；中國地質大學（武漢）地球內部多尺度成像湖北省重點實驗室開放基金項目（SMIL－2014－06）

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