基于混合正則化的最小二乘三維電阻率反演成像

尹洪東，楊懷章，薛亞茹，劉得軍

（中國石油大學地球物理與信息工程學院，北京102249）

基于混合正則化的最小二乘三維電阻率反演成像

尹洪東，楊懷章，薛亞茹，劉得軍

（中國石油大學地球物理與信息工程學院，北京102249）

基于經(jīng)典Tikhonov正則化（classical Tikhonov regularization）的最小二乘反演是三維電阻率反演的主要方法。對于電阻率分片連續(xù)的地質(zhì)體，由于光滑反演解的光滑性使得目標區(qū)域與背景區(qū)域間邊界模糊，不能很好地體現(xiàn)異常體的形態(tài)信息和位置信息，成像效果不好。將總變分正則化（total variation regularization）與經(jīng)典Tikhonov正則化結合，提出混合正則化反演方法，通過模型分析比較經(jīng)典Tikhonov正則化、TV正則化、混合正則化在反演結果上的不同，證明了引入混合正則化約束的最小二乘反演既保持了經(jīng)典Tikhonov正則化方法解的穩(wěn)定性，又具有TV正則化方法解的保邊緣性，有效地改善了三維電阻率成像效果。最后將混合正則化的反演方法應用到實際工程，驗證了該方法的有效性和實用性。

TV正則化；經(jīng)典Tikhonov正則化；混合正則化；三維電阻率反演

隨著電阻率法勘探的不斷發(fā)展，人們對電阻率反演問題越來越重視。國內(nèi)外學者對電阻率三維反演問題展開了研究并提出了不同的反演方法，其中Li［1］提出了基于Born近似的直流電法三維反演，對于電阻率變化較大的地質(zhì)體，此方法的反演精度很低，很難應用于復雜電性結構的反演。Zohdy［2］反演是針對溫納或施倫貝格測深曲線的解釋提出的，其原理是通過不斷調(diào)整初始模型直至模型測深曲線和實際測深曲線之間的擬合誤差達到最小，最終的模型參數(shù)即為反演結果，實際是一種最小二乘優(yōu)化反演［2-3］。目前最小二乘反演是主要的三維電阻率反演方法，Sasaki將有限單元法和最小二乘法結合，并且將光滑約束加入到最小二乘反演中，得到了較好的反演結果，奠定了使用有限單元法和最小二乘法進行電阻率反演的基礎［4］。黃俊革等［5-8］分別將“體積因子”、電阻率取值范圍、空間形態(tài)等作為先驗信息施加到電阻率三維反演中，極大地改善了三維電阻率反演問題的多解性，提高了深度較大的地質(zhì)體的分辨能力。吳小平等［9-12］分別在雅克比矩陣的求解、反演方程求解方法等方面對最小二乘反演方法進行了發(fā)展和完善，使得三維電阻率反演速度得到了很大的提高、反演效果有了顯著進步。相鵬［13］提出了一種改進的預條件非線性共軛梯度法，通過構建形狀更接近高斯牛頓Hessian矩陣的預條件算子提高反演精度和計算速度。另外隨著計算技術和最優(yōu)化理論的不斷發(fā)展，非線性反演方法也被引入到電阻率反演中，如程勃等［14-18］實現(xiàn)了遺傳算法、統(tǒng)計學建模法、BP神經(jīng)網(wǎng)絡等優(yōu)化方法在電阻維反演中的應用。上述的研究工作為三維電阻率反演奠定了基礎并且取得了很好的發(fā)展和完善，但對于電阻率值分片連續(xù)的三維地質(zhì)體，最小二乘反演解的光滑性使得反演的結果過度光滑，導致目標區(qū)域與背景區(qū)域邊界模糊，不能很好地體現(xiàn)出異常體的形態(tài)信息和位置信息，降低了成像的分辨率。為解決上述問題，筆者在最小二乘光滑反演的基礎上引入TV正則化方法。TV正則化方法首先由Rudin等提出并應用于圖像處理的圖像去噪中［19］。由于TV正則化函數(shù)可以很好地保持圖像的邊緣信息，因而在圖像重建、地震資料解釋等方面有著廣泛的應用［20-22］，并且已經(jīng)有學者將TV正則化方法應用到反演中去［23-24］，并得到了很好的結果。韓波等［25-26］將TV正則化與經(jīng)典的Tikhonov正則化的混合正則化方法應用到電阻率二維反演中，大大提高了異常體的成像分辨率。筆者在相關研究的基礎上，將混合正則化方法引入到三維電阻率反演中，針對電阻率值分片連續(xù)的地質(zhì)體模型進行算例分析。

1　三維電阻率成像基本原理

1.1數(shù)學模型

在直流電法勘探中，將直流電源的兩端通過電極向大地供電形成的地下電場可作為穩(wěn)定的電流場來處理。要計算地下穩(wěn)定電流場的電勢分布則必須從其滿足的微分方程和邊界條件出發(fā)，求解邊值問題。三維半空間電阻率的數(shù)學模型可以歸結為以下方程：

式中，σ為電導率；φ為電勢；δ為狄拉克函數(shù)；（x0，y0，z0）為源的空間位置；n為邊界外法線方向的單位矢量；r為源點到邊界點的徑向矢量；θ為n和r的夾角；Ω為半空間區(qū)域；Γs為半空間區(qū)域的地面邊界；?！貫榘肟臻g區(qū)域的無窮遠邊界。

1.2正演方法

考慮到有限體積法兼具有有限單元法和有限差分法的優(yōu)點，本文中采用有限體積法模擬三維空間內(nèi)的電阻率分布情況。首先采用單元中心方式的結構網(wǎng)格對三維區(qū)域進行離散，為方便起見將方程（1）中第一個方程的右端用s表示：

對式（2）進行體積分后得

其中離散電阻率σ=（σ1，σ2，…，σn），n表示有限體積單元的個數(shù)。將式（3）離散后得

式中，A為離散正演算子；u為離散電勢；q為點源位置矩陣。在實際工作中，只能測量到部分區(qū)域（如地表）的電位信息d，所以有

其中Q為投影矩陣。

1.3基于混合正則化的最小二乘反演

1.3.1最小二乘混合正則化反演的目標函數(shù)

反演問題是地球物理學中的核心問題，反演的目的是通過地面測量的數(shù)據(jù)獲取地球物理模型。對于三維反演問題，將模型剖分成Nx×Ny×Nz的三維網(wǎng)格，反演要求的參數(shù)就是各網(wǎng)格單元內(nèi)的電阻率值。采用最小二乘光滑反演目標函數(shù)可表示為

式中，d（m）為給定模型m通過正演得到的數(shù)據(jù)；dobs為觀測數(shù)據(jù)；α為正則化參數(shù)；Wm為光滑度矩陣；mref為參考模型。

相比于最小二乘法的目標函數(shù)，該方法加入了經(jīng)典Tikhonov正則化約束項此約束項減弱了反演問題近似解的不穩(wěn)定性，使模型參數(shù)的變化趨于光滑，最終得到光滑的反演結果。對于電阻率值分片連續(xù)的地質(zhì)體，只采用經(jīng)典Tikhonov正則化約束的最小二乘光滑反演往往會導致目標區(qū)域與背景區(qū)域邊界模糊，不能保持異常體的邊緣，為此將TV正則化方法引入到最小二乘反演中。TV正則化函數(shù)如下：

其中D為三維梯度算子矩陣，其作用是使三維模型m中的每一個元素都同其相鄰的元素做梯度運算。對于大小為Nx×Ny×Nz的剖分網(wǎng)格，D可以表示為

其中，i、j、k表示三維模型中每個單元編號。

通過模型試驗發(fā)現(xiàn)，由于TV正則化存在解的不穩(wěn)定性，導致成像時會產(chǎn)生虛假邊界。受文獻［25］、［26］啟發(fā)，本文中將經(jīng)典Tikhonov正則化和TV正則化兩種正則化方法結合后應用到最小二乘三維電阻率反演中來，此時目標函數(shù)可寫為

式中，α、β為正則化參數(shù)，兩者的值起到調(diào)節(jié)不同約束項的作用。

1.3.2基于迭代再加權的模型修正量的求解

在最小二乘意義下，反演方程的求解采用了將非線性問題線性化的方法，因此反演必須經(jīng)過多次迭代來修正模型參數(shù)，才能使得估計解^m不斷逼近真實解m*。模型參數(shù)修正量可以通過迭代再加權方法求解。迭代再加權方法［27］是解無約束優(yōu)化問題的重要方法，其核心思想是通過計算一系列權系數(shù)不斷更新的加權最小二乘問題得到目標函數(shù)的最優(yōu)解。求解時首先對目標函數(shù)（8）進行轉(zhuǎn)化，采用加權的方式將L1范數(shù)松弛為L2范數(shù)形式：

設第k次迭代后模型參量為m（k），此時的w=將d（m）、m在m（k）的鄰域用泰勒級數(shù)展開，并略去二次以上高階項可以近似為

將式（10）、（11）帶入到式（9）中，得

根據(jù)多元函數(shù)的極值理論，要求解模型修正量δm須對式（12）進行求導運算：

記J=d′（m（k））為偏導數(shù)矩陣，Wtv=w·D為再加權矩陣，令，得

這樣得出模型修正量為δm=H-1·g。

2　模型試驗及分析

為了驗證基于混合正則化的最小二乘三維電阻率反演方法兼有經(jīng)典Tikhonov正則化和TV正則化的優(yōu)點，利用RESINVM3D共享代碼［28］建立了兩個模型，合成觀測數(shù)據(jù)開展數(shù)值算例研究。兩個模型電極排列均采用四極E-SCAN高密度測量方式，其中電流源采用雙極供電，共41對，電位測量點93個，共有3813個測量數(shù)據(jù)。

2.1模型1及其分析

模型1是電阻率為40 Ω·m、大小為34 m×35 m×16 m的均勻半空間中存在兩個電導率為200 Ω ·m、大小為6 m×6 m×4 m的高阻異常體，兩個異常體并排放置，兩個異常體的頂部埋深都為2 m，如圖1所示（a）為模型1的水平剖面圖，（b）為模型1的x-z垂直剖面圖，（c）為模型1異常體所在的目標層的三維剖面圖。

圖1 模型1剖面圖Fig.1 Profiles of model 1

圖2為模型1反演結果的目標層剖面圖，其中圖2（a）是基于經(jīng)典Tikhonov正則化的光滑反演結果的目標層剖面圖，可以看出異常體各層形狀均呈橢圓形，且隨著深度增加，異常體邊緣逐漸擴散，目標區(qū)域與背景區(qū)域邊界模糊，成像的對比度不高；（b）是基于TV正則化反演結果的目標層剖面圖，可以看出目標層的第1、2層中異常體的邊緣較為清晰，但在第3、4層幾乎不能判斷出有異常體存在，其原因是TV正則化的保邊緣性使得反演過程中能夠得到跳躍性較大的參數(shù)部分，但是其解的不穩(wěn)定性造成了異常體的深度判斷異常；（c）和（d）分別是TV正則化參數(shù)β=0.000 05和β=0.000 01時的混合正則化反演結果的目標層剖面圖，相比經(jīng)典Tikhonov正則化約束和TV正則化約束，反演結果有了很大的改善：一方面目標區(qū)域與背景區(qū)域邊界清晰，能夠很好地分辨出異常體的邊緣，另一方面異常體內(nèi)部電阻率值比較一致，能夠很好地體現(xiàn)出均勻異常體的特性。由此可以看出：混合正則化既保留了經(jīng)典Tikhonov正則化解的穩(wěn)定性和光滑性，又具有TV正則化解的保邊緣性。

圖2 模型1反演結果目標層剖面圖Fig.2 Horizontal planes of target of three inversion results of model 1

圖3和圖4分別為模型1反演結果在4 m深處的水平剖面圖和y=17 m處的垂直剖面圖，其中（a）、（b）、（c）、（d）分別是經(jīng)典Tikhonov光滑反演、TV正則化反演、TV正則化系數(shù)β=0.00005和β= 0.00001時的混合正則化反演的結果。從圖3（a）和圖4（a）中可以清晰地看出：光滑反演使得異常體邊緣擴散嚴重，尤其在縱深方向，幾乎不能判斷異常體的真實邊界。TV正則化的保邊緣特性能夠在圖3（b）中清晰地看出，但由于TV正則化存在對目標區(qū)域定位不準確的缺點［26］，導致反演的結果中異常體出現(xiàn)重心“上漂”現(xiàn)象［12］（圖4（b））。圖3（c）、圖4（c）很好地體現(xiàn)出了混合正則化的優(yōu)越性：①吸收了TV正則化保邊緣特性，相比于經(jīng)典Tikhonov正則化異常體的邊緣擴散現(xiàn)象得到了很好的改善，能夠很好地分辨出異常體的邊緣；②繼承了經(jīng)典Tikhonov正則化解的穩(wěn)定性的特點，相比于TV正則化，目標區(qū)域定位更準確，異常體中心“上漂”現(xiàn)象消失。

圖3　模型1反演結果水平剖面圖（z=4 m）Fig.3 Horizontal plane of model 1（z=4 m）

圖4　模型1反演結果垂直剖面圖（y=17 m）Fig.4 Vertical plane of model 1（y=17 m）

在將兩種正則化項進行組合過程中，正則化系數(shù)的選取是十分重要的環(huán)節(jié)，針對模型1在保持經(jīng)典Tikhonov正則化項的系數(shù)α不變的前提下，調(diào)節(jié)TV正則化項的系數(shù)β，經(jīng)過多次試驗得出β的取值范圍在0.00001至0.0001之間時得到的結果較為理想：β值大于0.0001時，解會表現(xiàn)出TV正則化解的不穩(wěn)定性，且β值越大，解就越不穩(wěn)定；β值小于0.00001時，TV正則化解的保邊緣的特性不明顯，解的光滑性越來越顯著，且β值越小保邊緣的特性越不明顯。比較圖2中（c）和（d）可以看出，β= 0.00005時異常體的邊緣更加清晰，β=0.000 01時異常體的邊緣和圖2（a）更加接近。同樣，從圖3（c）和（d）、圖4（c）和（d）中能夠更好地體現(xiàn)出這一結論。

2.2模型2及其分析

模型2是一個大小為34 m×35 m×16 m、電阻率值為200 Ω·m的均勻半空間中存在兩個大小為6 m×6 m×4 m的異常體，其中左邊藍色異常體是電阻率為40 Ω·m的低阻體，右邊紅色異常體是電阻率為400 Ω·m的高阻體，兩個異常體并排放置，兩個異常體的頂部埋深都為2 m。圖5中（a）為模型2的水平剖面圖，（b）為模型2的x-z垂直剖面圖，（c）為模型2異常體所在的目標層的三維剖面圖。

圖6為模型2不同方法反演結果的目標層剖面圖，其中（a）是基于經(jīng)典Tikhonov正則化的光滑反演結果的目標層剖面圖，（b）是基于TV正則化反演結果的目標層剖面圖，（c）、（d）分別是TV正則化參數(shù)β=0.00005和β=0.00001時的混合正則化反演結果的目標層剖面圖。圖7和圖8分別為模型2反演結果在4 m深處的水平剖面圖和y=17 m處的垂直剖面圖，圖7和圖8中（a）、（b）、（c）、（d）分別是經(jīng)典Tikhonov光滑反演、TV正則化反演、TV正則化系數(shù)β=0.00005和β=0.00001時的混合正則化反演結果。

圖5　模型2剖面圖Fig.5 Profiles of model 2

圖6　模型2反演結果目標層剖面圖Fig.6 Horizontal planes of target of three inversion results of model 2

圖7　模型2反演結果水平剖面圖（z=4 m）Fig.7 Horizontal plane of model 2（z=4 m）

圖8　模型2反演結果垂直剖面圖（y=17 m）Fig.8 Corss-section of model 2（y=17 m）

從圖6（a）中可以看出：經(jīng)典Tikhonov光滑反演結果中各層高阻異常體電阻率都能夠達到350 Ω· m，與模型中高阻異常體電阻率值400 Ω·m十分接近，低阻異常體的電阻率值約為100 Ω·m，與模型中的低阻異常體電阻率值400 Ω·m也較為接近，但是在異常體的邊緣方面，無論是高阻異常體還是低阻異常體，其形狀都呈橢圓形，這與模型中異常體的形狀不符，同樣的結果在圖7（a）、8（a）中更能夠清晰地看出來。造成這種現(xiàn)象的原因是：經(jīng)典Tikhonov正則化解的光滑性使得電阻率值由高到低過度平滑，不能檢測出跳躍性較大的電阻率參數(shù)，成像后表現(xiàn)出異常體邊緣擴散的現(xiàn)象。圖6（b）是TV正則化反演的結果，其中第一層的高阻異常體電阻率約為300 Ω·m，極少數(shù)的值在500 Ω·m，其余三層幾乎看不出有高阻異常體的存在，低阻異常體也只能在第一、二層明顯地看出來，其值約為100 Ω· m，這樣的結果在異常體的電阻率值方面與模型較為接近，但是在異常體的位置方面與模型不符，這可能是由TV正則化解的不穩(wěn)定性造成的，這種不穩(wěn)定使得反演陷入局部最優(yōu)；另外，在保持異常體邊緣方面，TV正則化反演的結果要好很多，在圖6（b）、7（b）中能夠清晰地看出低阻異常體為長方形，體現(xiàn)出了TV正則化良好的保邊緣性。

同樣，在保持系數(shù)α不變的前提下，分別設TV正則化項的系數(shù)β為0.000 05和0.000 01，將兩種正則化項進行結合。從圖6、圖7、圖8的（c）、（d）中可以看出TV正則化系數(shù)β越大，TV正則化項的作用體現(xiàn)得越明顯；TV正則化系數(shù)β越小，TV正則化項的作用越小，經(jīng)典Tikhonov正則化項越起主導作用。

綜合上述兩個模型的反演結果可以得出：

（1）異常體邊緣信息方面。模型1和模型2中異常體的水平剖面和垂直剖面均為長方形，從反演結果可以看出，經(jīng)典Tikhonov正則化反演異常體水平剖面和垂直剖面近似為圓形或者橢圓形，并且邊緣較為發(fā)散；TV正則化反演異常體的水平剖面極為接近長方形，但是垂直剖面形狀不規(guī)則；混合正則化反演的水平剖面和垂直剖面都非常接近長方形，且異常體邊緣較為收斂。由此可見，無論是高阻體還是低阻體，混合正則化約束都能保留異常體邊緣信息，而且在刻畫異常體邊緣信息和減小邊界擴散方面要比經(jīng)典Tikhonov正則化反演有更好的效果。

（2）異常體位置信息方面。從垂直剖面圖中可以看出TV正則化反演出現(xiàn)異常體重心“上漂”現(xiàn)象，混合正則化反演要比經(jīng)典Tikhonov正則化反演得到的異常體的位置信息更為精確。模型1中兩個高阻體都位于2～6 m深處，經(jīng)典Tikhonov正則化反演結果大致位于3～8 m處，混合正則化反演的結果大致位于2～6 m處。模型2中低阻體和高阻體都位于2～6 m處，經(jīng)典Tikhonov正則化反演中低阻體大致位于2～10 m處，高阻體大致位于3～9 m，混合正則化反演中低阻體大致位于2～6 m深處，高阻體也大致位于2～6m深處。由此看出，混合正則化約束的最小二乘反演方法在體現(xiàn)異常體位置信息方面要比基于經(jīng)典Tikhonov正則化的最小二乘光滑反演方法更加精確。

（3）數(shù)值方面。從兩個模型反演結果的剖面圖中可以看出，經(jīng)典Tikhonov正則化反演中異常體的電阻率值比混合正則化約束最小二乘光滑反演更加接近模型的真實值，但是經(jīng)典Tikhonov正則化反演的結果中產(chǎn)生了假異常，混合正則化反演的結果中則不存在假異常。

3　實際工程應用

在某油田的試驗中，以某一試驗井的金屬套管為供電電極，以1 000 m以外的另一口井的金屬套管為回流電極，采用10 KW井地電位發(fā)送機發(fā)射雙極性方波信號，在地下介質(zhì)中形成一個垂直線電流發(fā)射源。以試驗井為圓心，以正北方為0°方向，以20°為角度間隔，布置18條長400 m的呈放射狀分布的徑向測線，測環(huán)徑向間距為50 m，在徑向測線與測環(huán)交點處共布置144個測量電極（圖9），在試驗井附近選擇一點作為電位參考點，利用井地電位接收機采集觀測井地表各電極的電位差數(shù)據(jù)。

圖9 測網(wǎng)布置示意圖Fig.9 Sketch map of arrangement of survey network

圖10為研究區(qū)測量的電位差分布圖，通過定性分析圖10可以得出，100 m測線電位差相對均勻，在150～400 m測線上300°～0°～80°方向有高電位差異常，說明該范圍上存在較大面積的低電阻率分布，相反方向則存在相對高電阻率分布。因此，在遠井帶的300°～0°～80°方向存在大面積低阻高含水區(qū)，其相反方向則為相對高阻區(qū)。圖11為目標層（800～1000 m）反演的電阻率分布圖，其中（a）為經(jīng)典Tikhonov正則化光滑反演的結果，其定位的低阻區(qū)域在270°～0°～140°方向，大大超出了定性分析的結果，（b）為混合正則化反演的結果，其定位的低阻區(qū)域在300°～0°～120°方向，相對于經(jīng)典Tikhonov正則化的反演結果，更接近定性分析的結果。

圖10　研究區(qū)測量的電位差Fig.10 Potential difference of work area

圖11　800～1000 m電阻率分布Fig.11 Resistivity distribution at 800-1000 m

4　結束語

將TV正則化約束方法引入到最小二乘電阻率反演中，通過兩種正則化的線性組合得到混合正則化的方法。通過對比試驗，混合正則化約束的最小二乘反演既保持了經(jīng)典Tikhonov正則化方法解的穩(wěn)定性又具有TV正則化方法解的保邊緣性，有效地改善了反演中異常體邊緣信息與位置信息不精確的問題，從而提高了辨別異常體位置信息與形態(tài)的精度，提高了成像的分辨率。在處理實際數(shù)據(jù)方面，受數(shù)據(jù)采集的精度和分辨率的影響，該方法只能分辨出大致的電阻率分布趨勢，具體的細節(jié)還要結合地質(zhì)資料進行綜合分析。

［1］ LI Y G，OLDENBURG D W.Approximate inverse mappings in DC resistivity problems［J］.Geophysical Journal International，1992，109（2）：343-362.

［2］ ZOHDY A A R.A new method for the automatic interpretation of schlumberger and wenner sounding curves［J］. Geophysics，1989，54（2）：245-253.

［3］ 王興泰，李曉芹.電阻率圖像重建的佐迪（Zohdy）反演及其應用效果［J］.物探與化探，1996，20（3）：228-233. WANG Xingtai，LI Xiaoqin.Resistivity image reconstruction with Zohdy inversion and its application effect［J］.Geophysical and Geochemical Exploration，1996，20（3）：228-233.

［4］ SASAKI Y.3-D resistivity inversion using the finite-element method［J］.Geophysics，1994，59（12）：1839-1848.

［5］ 黃俊革，阮百堯，鮑光淑.有限單元法三維電阻率最小二乘反演中存在問題的研究［J］.地質(zhì)與勘查，2004，40（4）：70-75. HUANG Junge，RUAN Baiyao，BAO Guangshu.Re-search on 3D resistivity least-square inversion using the finite element method［J］.Geology and Prospecting，2004，40（4）：70-75.

［6］ 劉斌，李術才，李樹忱，等.基于不等式約束的最小二乘法三維電阻率反演及其算法優(yōu)化［J］.地球物理學報，2012，55（1）：260-268. LIU Bin，LI Shucai，LI Shuchen，et al.3D electrical resistivity inversion with least-squares method based on inequality constraint and its computation efficiency optimization［J］.Chinese Journal of Geophysics，2012，55（1）：260-268.

［7］ 劉斌，聶利超，李術才，等.三維電阻率空間結構約束反演成像方法［J］.巖石力學與工程學報，2012，31（11）：2258-2268. LIU Bin，NIE Lichao，LI Shucai，et al.3D electrical resistivity inversion tomography with spatial structural constraint［J］.Chinese Journal of Rock Mechanics and Engineering，2012，31（11）：2258-2268.

［8］ LI S C，NIE L C，LIU B，et al.3D electrical resistivity inversion using prior spatial shape constrains［J］.Applied Geophysics，2013，10（4）：361-372

［9］ 吳小平，徐果明.利用共軛梯度法的電阻率三維反演研究［J］.地球物理學報，2000，43（3）：420-426. WU Xiaoping，XU Guoming.Study on 3-D resistivity inversion using conjugate gradient method［J］.Chinese Journal of Geophysics，2000，43（3）：420-426.

［10］ 底青云，王妙月.積分法三維電阻率成像［J］.地球物理學報，2001，44（6）：843-852. DI Qingyun，WANG Miaoyue.3D resistivity tomography by integral method［J］.Chinese Journal of Geophysics，2001，44（6）：843-852.

［11］ ZHANG J，MACKIE R L，MADDEN T R.Three-dimensional resistivity forward modeling and inversion using conjugate gradients［J］.Geophysics，1995，60（5）：1313-1325.

［12］ 屈有恒，張貴賓，趙連鋒，等.井地有限線源三維電阻率反演研究［J］.地球物理學進展，2007，22（5）：1393-1402. QU Youheng，ZHANG Guibin，ZHAO Lianfeng，et al. Study on 3D resistivity inversion for finite surface-borehole line current source［J］.Progress in Geophysics，2007，22（5）：1393-1402.

［13］ 相鵬.一種改進的二維MT預條件非線性共軛梯度反演方法［J］.中國石油大學學報：自然科學版，2014，38（4）：42-49. XIANG Peng.Two-dimensional MT inversion method based on an improved preconditioned nonlinear conjugate gradient algorithm［J］.Journal of China University of Petroleum（Edition of Natural Science），2014，38（4）：42-49.

［14］ 程勃，底青云.基于遺傳算法和統(tǒng)計學的電阻率測深二維反演研究［J］.地球物理學進展，2012，27（2）：788-795. CHENG Bo，DI Qingyun.2D inversion of resistivity sounding base on GA and statistic method［J］.Progress in Geophysics，2012，27（2）：788-795.

［15］ 程勃，底青云.復雜地電結構條件下統(tǒng)計學建模法電阻率測深二維反演［J］.地球物理學報，2014，57（3）：961-967. CHENG Bo，DI Qingyun.Statistic modeling 2D VES inversion in complex geo-electric structures［J］.Chinese Journal of Geophysics，2014，57（3）：961-967.

［16］ 徐海浪，吳小平.電阻率二維神經(jīng)網(wǎng)絡反演［J］.地球物理學報，2006，49（2）：584-589. XU Hailang，WU Xiaoping.2D resistivity inversion using the neural network method［J］.Chinese Journal of Geophysics，2006，49（2）：584-589.

［17］ 劉斌，王傳武，楊為民，等.基于并行改進遺傳算法的三維電阻率反演方法［J］.巖土工程學報，2014，36（7）：1252-1261. LIU Bin，WANG Chuanwu，YANG Weimin，et al.3D resistivity inversion using an improved parallel genetic algorithm［J］.Chinese Journal of Geotechnical Engineering，2014，36（7）：1252-1261.

［18］ LIU B，LI S C，NIE L C，et al.3D resistivity inversion using an improved genetic algorithm based on control method of mutation direction［J］.Journal of Applied Geophysics，2012，87：1-8.

［19］ RUDIN L I，OSHER S，F(xiàn)ATEMIE E.Nonlinear total variation based noise removal algorithms［J］.Physica D，1992，60（1/4）：259-268.

［20］ 王化祥，唐磊，閆勇.電容層析成像圖像重建的總變差正則化算法［J］.儀器儀表學報，2007，28（11）：2014-2018. WANG Huaxiang，TANG Lei，YAN Yong.Total variation regularization algorithm for electrical capacitance tomography［J］.Chinese Journal of Scientific Instrument，2007，28（11）：2014-2018.

［21］ 劉志文，潘曉露，李一民.L1范數(shù)的TV正則化超分辨率圖像重建［J］.微處理機，2012（3）：37-39. LIU Zhiwen，PAN Xiaolu，LI Yimin.L1 norm of total variation regularization based super resolution reconstruction for images［J］.Microprocessors，2012（3）：37-39.

［22］ 范文茹，王化祥，郝魁紅.基于兩步迭代TV正則化的電阻抗圖像重建算法［J］.儀器儀表學報，2012，33（3）：625-630. FAN Wenru，WANG Huaxiang，HAO Kuihong.Twostep iterative TV regularization algorithm for image reconstruction of electrical impedance tomography［J］. Chinese Journal of Scientific Instrument，2012，33（3）：625-630.

［23］ YUAN S Y，WANG S X，LI G F.Random noise reduction using Bayesian inversion［J］.Journal of Geophysics and Engineering，2012，9（1）：60-68.

［24］ ANAGAW A Y，SACCHI D M.Edge-preserving seismic imaging using the total variation method［J］.Journal of Geophysics and Engineering，2012，9（1）：138-146.

［25］ 韓波，竇以鑫，丁亮.電阻率成像的混合正則化反演算法［J］.地球物理學報，2012，55（3）：970-980. HAN Bo，DOU Yixin，DING Liang.Electrical resistivity tomography by using a hybrid regularization［J］.Chinese Journal of Geophysics，2012，55（3）：970-980.

［26］ 岳建惠.電阻率成像反問題的混合正則化方法研究［D］大連：大連海事大學數(shù)學系，2012. YUE Jianhui.Research of mixture regularization methods for EIT inverse problem［D］.Dalian：Department of Mathematics，Dalian Maritime University，2012.

［27］ 袁義明，孫晨，楊長春.基于迭代再加權最小二乘的地震資料稀疏反演方法［J］.地球物理學進展，2013，28（5）：2536-2546. YUAN Yiming，SUN Chen，YANG Changchun.Seismic data sparse inversion method based on iterative reweighted technique［J］.Progress in Geophysics，2013，28（5）：2536-2546.

［28］ PIDLISECKY A，HABER E，KNIGHT R.RESINVM3D：a 3D resistivity inversion package［J］.Geophysics，2007，72（2）：H1-H10.

（編輯 修榮榮）

Least square method based on hybrid regularization for 3D resistivity inversion imaging

YIN Hongdong，YANG Huaizhang，XUE Yaru，LIU Dejun
（College of Geophysics and Information Engineering in China University of Petroleum，Beijing 102249，China）

The smooth least-squares method based on the classical Tikhonov regularization method is the main method in 3D resistivity inversion.For the geological body whose resistivity is piecewise continuous，this method however fails to distinguish the target area from the background due to the smooth nature of the solution.In order to solve this problem，the classical Tikhonov regularization and the total variation regularization（TV regularization）were combined and introduced into the least squares inversion method.Three inversion results based on classical Tikhonov regularization，total variation regularization and the new hybrid regularization were compared.It is shown that the least square method based on hybrid regularization retains both the stability of the classical Tikhonov regularization and the edge-protected property of the total variation regularization.The improvement in 3D resistivity mapping is immediate，and then the method is further verified in practical engineering problems.

total variation regularization；classical Tikhonov regularization；hybrid regularization；3D electrical resistivity inversion

P 631

A

1673-5005（2015）05-0072-10

10.3969/j.issn.1673-5005.2015.05.010

2014-11-25

國家自然科學基金項目（41374151）

尹洪東（1963-），男，副教授，博士，研究方向為油田信號檢測、勘探技術及儀器設備的研發(fā)。E-mail：yhz1305@163.com。

引用格式：尹洪東，楊懷章，薛亞茹，等.基于混合正則化的最小二乘三維電阻率反演成像［J］.中國石油大學學報：自然科學版，2015，39（5）：72-81.

YIN Hongdong，YANG Huaizhang，XUE Yaru，et al.Least square method based on hybrid regularization for 3D resistivity inversion imaging［J］.Journal of China University of Petroleum（Edition of Natural Science），2015，39（5）：72-81.