王智，王程，方思南

(1.長(zhǎng)江大學(xué) 電子信息學(xué)院，湖北 荊州 434023；2.中煤科工集團(tuán) 西安研究院有限公司，陜西 西安 710077；3.長(zhǎng)江大學(xué) 地球物理與石油資源學(xué)院，湖北 武漢 430100)

0 引言

電磁法勘探是以地殼中不同巖石之間的電性差異為物理基礎(chǔ)(如導(dǎo)電性、導(dǎo)磁性、介電性以及電化學(xué)性質(zhì)的差異)，通過觀測(cè)研究人工或天然形成電磁場(chǎng)的時(shí)間和空間分布規(guī)律，從而達(dá)到尋找地下有用礦產(chǎn)資源、查明地下地質(zhì)構(gòu)造以及解決地質(zhì)問題的目的[1]。電磁法勘探方法主要分為基于幾何測(cè)深原理的直流電法、基于頻率域測(cè)深原理的大地電磁法(MT/AMT/CSAMT)與基于時(shí)間域測(cè)深原理的瞬變電磁法(TEM)[2]等。直流電阻率法是電法勘探的最經(jīng)典方法之一，在金屬礦產(chǎn)、非金屬礦產(chǎn)、煤田、油氣及地下水等資源勘探中獲得廣泛而有效的應(yīng)用，并已拓展到水文、工程、環(huán)境、考古等與國(guó)民經(jīng)濟(jì)建設(shè)、人類社會(huì)生活密切相關(guān)的探測(cè)領(lǐng)域[3]。直流電阻率法有多種靈活的觀測(cè)方式，如單極、偶極及多極裝置的電剖面法、電測(cè)深法等，井地電阻率法是在井中供電，在地面接收電磁場(chǎng)的一類電法，供電裝置總是放入鉆孔深部使其接近于被探測(cè)的對(duì)象，從而增大電流強(qiáng)度或被接收的異常響應(yīng)[4-5]，井地電阻率法主要用于金屬礦礦山的二次資源勘查及油藏邊界的預(yù)測(cè)[6]等，通常采集裝置沿著具有不同電極陣列的平行線在網(wǎng)格上進(jìn)行數(shù)據(jù)采集，并使用3D反演算法進(jìn)行解釋。隨著計(jì)算機(jī)及數(shù)值計(jì)算技術(shù)的發(fā)展，電阻率三維正反演算法分別在網(wǎng)格的剖分(結(jié)構(gòu)化[7-12]、非結(jié)構(gòu)化[13-21])，數(shù)值方法(有限差分法[22-27]、有限元法[28-29])，目標(biāo)函數(shù)的求解(高斯—牛頓法(GN)[3, 19, 30-37]、擬牛頓法(QN)[38-54]、非線性共軛梯度法(NLCG)[45, 51, 55-63]等)等方面取得很大進(jìn)展。

由于結(jié)構(gòu)化單元的網(wǎng)格不能適應(yīng)復(fù)雜地形下的電阻率三維反演，因此非結(jié)構(gòu)有限元方法在復(fù)雜地形電阻率三維數(shù)值模擬中取得極大成功，非結(jié)構(gòu)化網(wǎng)格具有單元質(zhì)量可控、允許局部加密、能夠模擬復(fù)雜幾何模型等優(yōu)點(diǎn)，使得三維非結(jié)構(gòu)有限單元求解效率大幅提高。同等計(jì)算精度情況下，非結(jié)構(gòu)化網(wǎng)格相對(duì)于結(jié)構(gòu)化規(guī)則網(wǎng)格，計(jì)算時(shí)間和存儲(chǔ)量均可下降約一個(gè)數(shù)量級(jí)[3]。在電阻率三維反演中，由于反演參數(shù)多且數(shù)據(jù)量大，雅可比矩陣(偏導(dǎo)數(shù)矩陣)對(duì)巨大的計(jì)算量和存儲(chǔ)需求難于克服，許多學(xué)者研究可避開計(jì)算雅可比矩陣的反演算法，Zhang等[26]、吳小平等[36]引入共軛梯度方法，解決了電阻率三維正演效率及三維反演中雅可比矩陣求取、存儲(chǔ)兩方面問題，實(shí)現(xiàn)比較快速有效地電阻率三維反演。在三維數(shù)據(jù)反演中采用的最優(yōu)化方法主要有非線性共軛梯度法(NLCG)、高斯—牛頓法(GN)和擬牛頓法(QN)。NLCG和QN只需要目標(biāo)函數(shù)的梯度信息，無需顯式計(jì)算靈敏度矩陣，由于僅僅使用了目標(biāo)函數(shù)的一階導(dǎo)數(shù)(梯度)信息，只有線性收斂性，模型需要反復(fù)更新很多次才能收斂到理想的值，而GN法同時(shí)使用了目標(biāo)函數(shù)的一階導(dǎo)數(shù)(Gradient)和二階導(dǎo)數(shù)(Hessian)的信息，表現(xiàn)出近似二次收斂性，所需的模型更新次數(shù)遠(yuǎn)遠(yuǎn)少于NLCG和QN法，每次高斯牛頓迭代的法方程使用預(yù)條件共軛梯度法(PCG)求解，不僅可以避免顯示求解和存儲(chǔ)靈敏度矩陣，還能減少PCG的迭代次數(shù)來減小反演的計(jì)算時(shí)間。由于反演的非線性特征，這種非精確求解對(duì)于線性化迭代反演不僅減小了計(jì)算量，而且還避免了每次GN迭代中對(duì)法方程的過度求解，因此在反演搜索過程中有機(jī)會(huì)跳出局部極小值[33]。

壓制電阻率三維反演多解性的有效方法是發(fā)展約束反演，將已有的先驗(yàn)信息通過數(shù)學(xué)手段加入到反演過程中，提高反演精度。Sasaki Y[21]提出了光滑約束，使得相鄰網(wǎng)格之間的電阻率差異極小，光滑過渡，這種約束具有天然的合理性，已經(jīng)成為電阻率反演中常用的一種約束形式；Kim等[64]按照光滑約束的施加形式，定義了一個(gè)新的參數(shù)向量來表征模型參數(shù)，將不等式約束施加到反演方程中，獲得了不錯(cuò)的效果；Li等[65]利用約束最優(yōu)化里的內(nèi)點(diǎn)罰函數(shù)法將不等式約束條件加入到目標(biāo)函數(shù)里，提高了反演成像的精度；黃俊革等[66]和宛新林等[67]分別對(duì)粗糙度矩陣進(jìn)行改進(jìn)，提高了三維反演的效果；劉斌等[68]在反演目標(biāo)函數(shù)中引入自適應(yīng)加權(quán)光滑約束，控制深部分辨率并在其文獻(xiàn)中[69-70]討論了不等式和先驗(yàn)結(jié)構(gòu)約束條件的加入方法，去除了反演成像中的假異常和多余構(gòu)造，明顯地壓制了電阻率探測(cè)反演的多解性問題。

電磁探測(cè)反演是典型的不適定問題，易造成反演結(jié)果的多解性，不適定性是反演自身固有的特征，沒有求解的附加信息這一本質(zhì)困難是很難克服的[71]，如何利用相關(guān)技術(shù)引入先驗(yàn)信息實(shí)現(xiàn)地下結(jié)構(gòu)的高精度成像是未來電磁法重要的發(fā)展方向之一[72]。電阻率的變化范圍是一個(gè)很重要的先驗(yàn)信息，如何將不等式約束施加到反演方程中是一個(gè)關(guān)鍵性的難題。Kim等[64]采取的方法是一種經(jīng)驗(yàn)式施加方式，在反演過程中該約束項(xiàng)繁瑣、不容易理解；Li等[65]采用內(nèi)點(diǎn)罰函數(shù)法，在目標(biāo)函數(shù)中加入對(duì)數(shù)形式的障礙函數(shù)，這種施加形式簡(jiǎn)單，但是內(nèi)點(diǎn)法要求初始值的選取必須在可行域內(nèi)，并且初始點(diǎn)應(yīng)選擇一個(gè)離約束邊界較遠(yuǎn)的可行點(diǎn)，太靠近某一約束邊界會(huì)造成求解的困難。如果約束條件較為復(fù)雜時(shí)，選擇初始點(diǎn)有一定難度，并且內(nèi)點(diǎn)法只能解決約束為不等式情形，外點(diǎn)法可以解決約束為等式和不等式混合的情形，外點(diǎn)法對(duì)初始點(diǎn)也沒有要求。結(jié)合前人研究與上述反演存在的問題，本文將表征模型參數(shù)變化范圍的不等式約束作為先驗(yàn)信息以外點(diǎn)罰函數(shù)的方式引入到目標(biāo)函數(shù)中[73]，選擇井地觀測(cè)方式、非結(jié)構(gòu)化四面體的有限單元法與高斯牛頓法(GN)相結(jié)合方式實(shí)現(xiàn)了井地電阻率法的三維約束反演算法，通過典型三維理論模型合成觀測(cè)數(shù)據(jù)來驗(yàn)證所開發(fā)的算法的穩(wěn)定性及有效性。

1 三維正演理論

1.1 控制方程

三維直流電阻率法點(diǎn)電流源場(chǎng)總電位滿足的偏微分方程及其邊界條件為[9, 16, 74]:

(1)

其中：σ是巖石電導(dǎo)率；u是電位；I是電流強(qiáng)度；δ是狄拉克函數(shù)；rA是源點(diǎn)位置；ω是源點(diǎn)對(duì)地下區(qū)域Ω的張開角，若源點(diǎn)在地面，則ωA=2π；若源點(diǎn)在地下，ωA=4π，n是邊界Γ∞邊界上的外法向向量；r是點(diǎn)電源到邊界點(diǎn)的距離；Γs是自然邊界條件；?！奘菬o窮遠(yuǎn)邊界；θ是無窮遠(yuǎn)邊界上點(diǎn)的外法向向量n與場(chǎng)源到邊界點(diǎn)距離r之間的夾角。采用加權(quán)余量法可導(dǎo)出式(1)對(duì)應(yīng)的變分問題的積分方程:

(2)

計(jì)算區(qū)域采用四面體剖分與線性差值，最終形成的大型稀疏對(duì)稱線性方程組，矩陣表達(dá)形式如下：

Ku=P

，

(3)

其中：K是一個(gè)n×n的對(duì)稱矩陣；u為n×1列向量,表示三維網(wǎng)格節(jié)點(diǎn)上的電位向量；P是含有場(chǎng)源信息的列向量；P=(0,…,I,…,0)，為了節(jié)約內(nèi)存采用行壓縮存儲(chǔ)(compressed sparse row or compressed row storage)處理系數(shù)矩陣K，使用對(duì)稱超松弛預(yù)條件共軛梯度算法[11, 20](SSORPCG)求解方程組(3)。

1.2 算法驗(yàn)證

為了驗(yàn)證本文正演程序的正確性，本文設(shè)計(jì)了一個(gè)均勻半空間中埋藏球體的模型，如圖1a所示，該模型為均勻半空間中嵌入一個(gè)球形異常體[16]，球體的半徑為R=2.25 m，球體中心坐標(biāo)為(0,0,-4.5),球體電阻率值ρ0=1 Ω·m，地層圍巖電阻率值為ρ1=10 Ω·m，點(diǎn)電流源A坐標(biāo)在(-5,0,0)處，電流大小為1 A,測(cè)量電極M沿x方向，間距為0.25 m。整個(gè)計(jì)算區(qū)域大小為500 m×500 m×300 m，目標(biāo)區(qū)域大小為100 m×100 m×100 m，網(wǎng)格節(jié)點(diǎn)總數(shù)為273 582，網(wǎng)格單元總數(shù)為1 625 838，網(wǎng)格剖分效果見圖1b～1d，分別在源點(diǎn)、測(cè)點(diǎn)與球體處進(jìn)行了加密處理。采用pole-pole裝置，由本文的有限元正演程序計(jì)算得到測(cè)點(diǎn)處電位值，并與解析解進(jìn)行了比較，如圖2所示，本文的計(jì)算結(jié)果與解析解有很好的一致性，最大誤差小于0.34%(見圖3)。

a—球體模型示意；b—球體模型網(wǎng)格剖分放大效果；c—點(diǎn)與測(cè)點(diǎn)網(wǎng)格加密效果；d—網(wǎng)格剖分示意

圖2 地下球體的解析解與數(shù)值解對(duì)比

圖3 解析解與數(shù)值解的相對(duì)誤差

2 三維反演理論

2.1 目標(biāo)函數(shù)

依據(jù)Tikhonov正則化理論，最小二乘意義下的反演目標(biāo)函數(shù)為[3, 19, 31]:

Φ(m)=Φd(m)+λΦm(m)=

(4)

其中：F(m)為正演響應(yīng)函數(shù)；m為模型參數(shù)(mi,i=1,2,…,M)，dobs(dj,j=1,2,…,N)為觀測(cè)數(shù)據(jù)，Wd為N×N階觀測(cè)數(shù)據(jù)誤差加權(quán)矩陣，Wd=diag(1/σ1,1/σ2,…,1/σn)，σi為第i個(gè)觀測(cè)數(shù)據(jù)的標(biāo)準(zhǔn)差；Wm為M×M階模型加權(quán)矩陣，通常用模型單元的離散差分算子來定義，一般取一階正則化約束。λ為正則化因子，用來平衡數(shù)據(jù)擬合和模型光滑度的權(quán)重；mref為含有模型參數(shù)先驗(yàn)信息的參考模型。反演的觀測(cè)參數(shù)dobs為pole-pole的電位值，模型參數(shù)為單元的電導(dǎo)率值σ，由于變化范圍很大，為了反演的穩(wěn)定性，通常采用對(duì)數(shù)來標(biāo)定觀測(cè)數(shù)據(jù)及模型參數(shù)，即d=lnφobs,m=lnσ。光滑約束的施加請(qǐng)參考文獻(xiàn)[62]，不等式約束施加的原理詳見文獻(xiàn)[62-63]，常規(guī)電阻率反演目標(biāo)函數(shù)中構(gòu)造出帶有不等式約束的反演目標(biāo)函數(shù)見式(5)：

，

(5)

式中：mmin，mmax為每個(gè)單元網(wǎng)格電阻率的上下限。式(5)即為本文帶不等式約束條件的三維電阻率反演的目標(biāo)函數(shù)，與常規(guī)三維電阻率反演目標(biāo)函數(shù)相比，該方法增加了不等式約束項(xiàng)，該表達(dá)式簡(jiǎn)單、容易理解、易于實(shí)現(xiàn)。目前常用于三維電阻率反演的最優(yōu)化方法有擬牛頓法(QN),高斯牛頓法(GN)、非線性共軛梯度法(NLCG)以及這些方法的變種[33,75]。對(duì)目標(biāo)函數(shù)(5)進(jìn)行二階Taylor展開，忽略高階項(xiàng)得到高斯—牛頓法(Gauss-Newton，GN)的法方程(normal equation)為：

H(m)△(m)=-g(m)

，

(6)

其中：g(m)為目標(biāo)函數(shù)的一階導(dǎo)數(shù)(梯度矩陣，Gradient)，H(m)為目標(biāo)函數(shù)的二階近似導(dǎo)數(shù)(海森矩陣,Hessian)，分別表示如下:

，

(7)

其中：J(m)為靈敏度矩陣(正演響應(yīng)函數(shù)F(m)的一階導(dǎo)數(shù)),Wc為不等式約束矩陣：其元素Wc(i,j)取值如下:

(8)

2.2 模型更新

利用CG法求得法方程(6)的解Δm后，需要進(jìn)行線性搜索來獲得每次GN迭代的模型更新量，更新公式為：

mk+1=mk+αΔm

，

(9)

其中：α為模型的更新步長(zhǎng)，步長(zhǎng)α的精確線搜索往往需要大量的計(jì)算，一般采用非精確一維線性搜索方法得到步長(zhǎng)，常用的策略是采用Wofe-Powell準(zhǔn)則來進(jìn)行非精確線性搜索得到更新步長(zhǎng)α，Wofe-Powell準(zhǔn)則包括充分下降條件與曲率條件，公式為:

(10)

其中：Φ為正演算子，mk為第k次反演迭代的模型參數(shù)；αk為迭代步長(zhǎng)；pk為搜索方向；c1與c2為常數(shù)，滿足0

2.5 反演流程

GN-CG法反演流程如下[73, 76]：

GN-CG反演算法1) 給定初始模型m0與參考模型mref,設(shè)置收斂系數(shù)ε,正則化因子λ和懲罰因子μ,給定光滑度矩陣Wm;2) 開始GN-CG 迭代, 設(shè)k=1,2,…,NmaxGN;3) 計(jì)算Δdk=dobs-F(mk),并計(jì)算r0=-gk;gk=-JTkWTdWdΔdk+λWTmWm(m-mref)+μ[min(m-mmin,0)-min(mmax-m,0)];4) 開始CG 迭代, 令Δm=0,p0=r0,設(shè)i=1,2,…,Nmaxcg; 4-1)計(jì)算Hk=JTkWTdWdJk+λWTmWm+μWc, 得到Hkpi; 4-2)計(jì)算ti=rTiripTiHpi; 4-3)計(jì)算Δmi+1=Δmi+tipi; 4-4)計(jì)算ri+1=ri-tiHkpi; 4-5)計(jì)算βFRi=rTi+1ri+1rTiri; 4-6)計(jì)算pi+1=ri+1+βipi;5) CG迭代結(jié)束, 求得Δm,根據(jù)Wofe-Powell條件(式12)進(jìn)行一維搜索得到最優(yōu)步長(zhǎng)αk,更新模型:mk+1=mk+αkΔmk;6) 如果RMS<ε 則停止GN-CG迭代,輸出mk,否則令k=k+1 轉(zhuǎn)到步驟3);

3 合成數(shù)據(jù)反演

為了驗(yàn)證開發(fā)的基于不等式約束井地電阻率三維反演算法的正確性與有效性，本文設(shè)計(jì)了2個(gè)三維模型來分別模擬井地電阻率法的測(cè)量情況，并對(duì)這2個(gè)三維模型的合成數(shù)據(jù)進(jìn)行三維反演測(cè)試。本文的所有計(jì)算均在小型計(jì)算機(jī)上完成，其配置CPU為Intel I7-4712MQ，主頻為2.3 GHz，內(nèi)存16 G，正反演程序均采用Inter Fortran編譯運(yùn)行，非結(jié)構(gòu)化四面體網(wǎng)格剖分采用Gmsh-4.8.4[77]，非結(jié)構(gòu)化四面體網(wǎng)格的顯示采用Paraview-5.6.0[78]。

3.1 平坦地形下單個(gè)異常體

長(zhǎng)方體模型如圖4a、4b所示，圍巖電阻率ρ0=100 Ω·m，低阻體電阻率ρ1=10 Ω·m，低阻長(zhǎng)方體大小為10 m×10 m×5 m，頂面埋深h=5 m，底面埋深為10 m。藍(lán)色三角形代表鉆井井口位置(18,25,0)與(32,25,0)，異常體離兩邊鉆井的水平距離d=2 m。本次反演采用的測(cè)線范圍為10～40 m，測(cè)線之間的間距為2 m，每條測(cè)線上16個(gè)測(cè)點(diǎn)，測(cè)點(diǎn)M(測(cè)量電極)沿x方向的位置坐標(biāo)為x=10～40 m，y=10～40 m，間距為2 m；源點(diǎn)A(供電電極)移動(dòng)范圍在井下-4～-20 m，間隔為4 m，一共設(shè)置10個(gè)點(diǎn)源；地表測(cè)量裝置位置見圖4c。為了測(cè)試反演算法的穩(wěn)定性，避免反演偏差，反演網(wǎng)格采用不同于合成數(shù)據(jù)所用的正演網(wǎng)格：正演網(wǎng)格分別在測(cè)點(diǎn)、源點(diǎn)與異常體附近進(jìn)行網(wǎng)格加密；反演網(wǎng)格在源點(diǎn)與測(cè)點(diǎn)處進(jìn)行加密，不僅提高了正演的精度，還減少了源點(diǎn)附近的數(shù)值模擬誤差，測(cè)點(diǎn)局部加密示意圖見圖4d。圖4e、圖4f分別為正演與反演所使用的非結(jié)構(gòu)化四面體網(wǎng)格，正演網(wǎng)格單元總數(shù)為1 288 357。反演網(wǎng)格中網(wǎng)格單元總數(shù)為1 090 207，通過三維有限元正演程序共得到2 560個(gè)一次場(chǎng)的”電位實(shí)測(cè)數(shù)據(jù)”，正演一次約130 s，反演一次約14 400 s。反演參數(shù)選擇如下：正則化參數(shù)λ=0.05，懲罰因子μ=1.0，在反演過程中始終保持不變，反演終止的收斂系數(shù)ε=1.0×10-5，反演的迭代次數(shù)10次，反演過程中RMS變化如圖5所示。從圖6可以看出，數(shù)據(jù)擬合差穩(wěn)定下降，說明本文的帶約束的GN法對(duì)于三維井地電阻率法的數(shù)據(jù)收斂性較好，圖5分別給出了傳統(tǒng)正則化反演與施加不等式約束反演結(jié)果的xoy與xoz剖面，兩種反演方法得到的模型與真實(shí)模型吻合較好，地下低阻異常體的位置與電阻率值基本一致，驗(yàn)證了本文反演方法的有效性。

a—xoy水平截面; b—xoz垂直斷面;c—三維異常體模型；d—測(cè)點(diǎn)加密示意;e—正演網(wǎng)格;f——反演網(wǎng)格

圖5 迭代次數(shù)與RMS的關(guān)系

a—傳統(tǒng)電阻率反演結(jié)果z=-7 m處xoy切片;b—傳統(tǒng)電阻率反演結(jié)果y=25 m處xoz切片;c—施加不等式約束的反演結(jié)果z=-7 m處xoy切片;d—施加不等式約束的反演結(jié)果y=25 m處xoz切片;e—傳統(tǒng)電阻率反演結(jié)果的三維剖面(電阻率小于50 Ω·m);f—施加不等式約束反演結(jié)果的三維剖面(電阻率小于50 Ω·m)

3.2 五個(gè)異常體模型

為了測(cè)試所開發(fā)的反演算法對(duì)于復(fù)雜模型的有效性，考慮淺地表存在的電阻率異常干擾，在背景為1 000 Ω·m的均勻半空間中放置5個(gè)電阻率異常體[79-80]，其中3個(gè)長(zhǎng)方體頂界面位于淺地表以下10 m(S1,S2,S3)，用于模擬淺地表的電阻率異常干擾，另外兩個(gè)頂界面埋深分別為60 m和95 m的長(zhǎng)方體(B1，B2)模擬探測(cè)目標(biāo)體，各項(xiàng)參數(shù)見表1，三維模型見圖7a。整個(gè)模型的研究區(qū)域(目標(biāo)區(qū)域與邊界區(qū)域)x、y、z這3個(gè)方向的剖分范圍分別是(-1 000，2 000)、(-1 000，2 000)和(-1 500，0)，目標(biāo)區(qū)域x、y、z這3個(gè)方向的剖分范圍分別是(0，1 000)、(0，1 000)和(-500，0)。采用的測(cè)線范圍為50～950 m，測(cè)線之間的間距為50 m，每條測(cè)線上19個(gè)測(cè)點(diǎn)，測(cè)點(diǎn)M(測(cè)量電極)沿x方向的位置坐標(biāo)為x=50～950 m，y=50～950 m，間距為50 m；井口坐標(biāo)分別為(500，200，0)與(500，800，0)，源點(diǎn)A(供電電極)移動(dòng)范圍在井下-40～-360 m，間隔為80 m，一共設(shè)置10個(gè)點(diǎn)源。地表測(cè)量裝置與點(diǎn)源位置分別見圖7b與7c。正演網(wǎng)格單元總數(shù)為2 511 215，反演網(wǎng)格中網(wǎng)格單元總數(shù)為1 037 303，通過三維有限元正演程序共得到3 610個(gè)一次場(chǎng)的”電位實(shí)測(cè)數(shù)據(jù)”，正演一次約275 s，反演一次約17 039 s。

a—三維異常體模型位置; b—地表數(shù)據(jù)觀測(cè)系統(tǒng);c—垂直剖面yoz; d—非結(jié)構(gòu)化網(wǎng)格加密示意

表1 長(zhǎng)方體參數(shù)

淺層地震勘探、地質(zhì)雷達(dá)探測(cè)法對(duì)異常體界面識(shí)別和定位方面具有較好的效果，在電阻率探測(cè)時(shí)利用這些方法獲得的異常體位置與形態(tài)信息是一種極為重要的先驗(yàn)構(gòu)造信息[62, 63, 69]，因此在三維井地電阻率反演中將構(gòu)造位置的先驗(yàn)信息作為局部約束加入。由其他地球物理探測(cè)方法假設(shè)可知異常體的構(gòu)造位置與形態(tài)信息，由探區(qū)礦石標(biāo)本物性參數(shù)可得：分別在①250 m≤x≤350 m，250 m≤y≤350 m，-50 m≤z≤-10m處異常體的電阻率值在220 Ω·m以內(nèi)，將此處網(wǎng)格單元的電阻率上限設(shè)為220 Ω·m，下限設(shè)為0；②300 m≤x≤700 m、400 m≤y≤500 m、-50 m≤z≤-10 m處異常體的電阻率值在120 Ω·m以內(nèi)，將此處網(wǎng)格單元的電阻率上限設(shè)為120 Ω·m，下限設(shè)為0 Ω·m；③350 m≤x≤700 m、650 m≤y≤700 m、-50 m≤z≤-10 m處異常體的電阻率值在2 200 Ω·m以內(nèi)，將此處網(wǎng)格單元的電阻率上限設(shè)為2 200 Ω·m，下限設(shè)為1 800 Ω·m；④300 m≤x≤650 m、300 m≤y≤400 m、-260 m≤z≤-60 m處異常體的電阻率值在2 200 Ω·m以內(nèi)，將此處網(wǎng)格單元的電阻率上限設(shè)為2 200 Ω·m，下限設(shè)為1 800 Ω·m；⑤400 m≤x≤500 m、500 m≤y≤650 m、-275 m≤z≤-95 m處異常體的電阻率值在120 Ω·m以內(nèi)，下限設(shè)為0 Ω·m。各項(xiàng)反演參數(shù)的設(shè)置與模型1的反演參數(shù)設(shè)置一樣。圖8c與圖8d為施加不等式約束后的電阻率反演結(jié)果，由反演結(jié)果可以看出，5個(gè)異常體水平和橫向的分辨率較傳統(tǒng)反演結(jié)果(圖8a，圖8b)有了明顯提高，與原模型在電阻率值、規(guī)模、位置、形狀等特征方面均基本一致，帶不等式約束的正則化反演能夠得到更為精確的三維反演成像結(jié)果。

a—傳統(tǒng)正則化反演x=400 m處yoz剖面; b—傳統(tǒng)正則化反演x=600 m處yoz剖面;c—約束反演x=400 m處yoz剖面; d—約束反演x=600 m處yoz剖面

4 結(jié)論

本文實(shí)現(xiàn)了一種以外點(diǎn)懲罰函數(shù)法的方式將帶有電阻率取值范圍的不等式約束施加到井地電阻率三維反演方程中的GN-CG反演算法，正演采用有限元法與非結(jié)構(gòu)化網(wǎng)格對(duì)控制方程離散化并使用SSORPCG法求解離散所得的大型線性方程；反演中采用高斯牛頓最優(yōu)化算法，在目標(biāo)函數(shù)中加入不等式與構(gòu)造約束來降低反演的多解性。以上這些技術(shù)的綜合使用，使得本文帶約束的GN-CG算法穩(wěn)定且準(zhǔn)確，適用于陸地與井中的電阻率法勘探情況。通過對(duì)理論模型的反演試算，驗(yàn)證了本文反演算法的可靠性與有效性，帶約束的三維井地電阻率反演方法能得到更精確的反演結(jié)果；懲罰因子的取值與不等式約束的范圍給定是兩個(gè)難點(diǎn)，本文懲罰因子的取值采用的是實(shí)驗(yàn)法，過大或者過小的懲罰因子都會(huì)影響反演的效果，有時(shí)甚至?xí)霈F(xiàn)不收斂的情況，可以通過結(jié)合其他物探方法給出異常體較為準(zhǔn)確的邊界范圍以及電阻率值的上下限，給出的范圍越精確，反演結(jié)果就越準(zhǔn)確。