侯宇健，李桐林，張镕哲，陳漢波

吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院，長春 130026

0 引言

激發(fā)極化法是以巖、礦石的激電效應(yīng)差異為物質(zhì)基礎(chǔ)，通過觀測和研究大地激電效應(yīng)，探查地下地質(zhì)情況的一種電法分支方法。該方法在尋找金屬礦、能源和工程地質(zhì)等不同領(lǐng)域均有廣泛的應(yīng)用，尤其在尋找浸染狀礦產(chǎn)資源方面效果顯著[1]。

20世紀60、70年代，國內(nèi)學(xué)者開始了解激發(fā)極化法并應(yīng)用[2]。最初的反演方式多是用擬斷面圖進行定性解釋或用一維方法來解釋。阮百堯等[3]在二維電阻率/激發(fā)極化法反演程序的基礎(chǔ)上進行了改進，使用算子線性化反演方法在國內(nèi)首次實現(xiàn)了二維電阻率/極化率的反演。吳小平[4]提出了用共軛梯度方法的激發(fā)極化三維快速反演，實現(xiàn)了極化率的三維反演。

地球物理反演問題存在多解性[5]，為減少多解性的影響，一般期望電阻率與極化率的結(jié)構(gòu)相同。Gallardo[6]最早提出了多參數(shù)聯(lián)合反演方法，實現(xiàn)了兩種物性在結(jié)構(gòu)上互相約束的反演方式。李兆祥等[7]實現(xiàn)了基于交叉梯度約束的激發(fā)極化二維反演，首次進行了極化率電阻率的聯(lián)合反演。閆政文等[8]實現(xiàn)了多參數(shù)的三維聯(lián)合反演。

筆者在前人工作的基礎(chǔ)上，基于等效電阻率公式，實現(xiàn)了三維電阻率/極化率正反演和聯(lián)合反演。電阻率正演使用有限體積法進行數(shù)值模擬，反演過程中使用不精確牛頓法進行求解。使用共軛梯度法減少了需要占用存儲空間，使三維反演速度和效率更高。推導(dǎo)了三維規(guī)則網(wǎng)格下的交叉梯度項及其梯度公式，提出了一種間接實現(xiàn)電阻率極化率聯(lián)合反演的方式，并對各反演結(jié)果進行了分析。

1 正演理論

1.1 電阻率法正演

本文采用單元中心控制體積法[9]進行正演，對三維空間進行剖分，將地下空間離散為規(guī)則六面體單元，每一塊的控制體積的基本方程為：

▽·(σ▽U)=-s

(1)

對其在控制體積Vijk內(nèi)積分：

(2)

由高斯定理得：

(3)

等式右邊為控制體積表面積分，等于各個面積分和：

(4)

經(jīng)推導(dǎo)可得出控制方程：

λi+1, j, kUi+1, j, k+λi-1, j, kUi-1, j, k+λi, j+1, kUi, j+1, k+

λi, j-1, kUi,j-1, k+λi, j, k+1Ui, j, k+1+λi, j, k-1Ui, j, k-1-λi, j, kUi, j, k=Q

(5)

式中：

λi, j, k=-(λi+1, j, k+λi-1, j, k+λi, j+1, k+λi, j-1, k+λi, j, k+1+λi, j, k-1)

(6)

以此可以得到正演系數(shù)矩陣K，于是正演問題變成解決：

KU=q

(7)

線性方程組問題，q是源項向量。U為地下各點電位。

1.2 極化率正演方法

同樣，用真實電阻率模型做一次正演，算出不帶有激發(fā)極化效應(yīng)的正演響應(yīng)結(jié)果ρs。

利用等效電阻率公式，即可求出極化率的正演響應(yīng):

(8)

2 反演理論

2.1 電阻率聯(lián)合反演方法

使用高斯-牛頓法建立目標函數(shù)，交叉梯度約束反演的目標函數(shù)為：

(9)

式中：

(10)

式中：D矩陣是數(shù)據(jù)誤差矩陣；其中S是實測數(shù)據(jù)的誤差；ε的存在是為保證分母不為0。β為正則化因子，是控制模型光滑程度的系數(shù)。W是模型光滑度矩陣，mref是期待的參考模型，d(m)為正演數(shù)據(jù)列向量，dobs為實測數(shù)據(jù)列向量。t為交叉梯度項，定義式為：

t(x,y,z)=▽mp(x,y,z)×▽mr(x,y,z)=(tx,ty,tz)

(11)

式中：mp，mr是兩種物性的反演結(jié)果，如果反演初始模型是均勻的，那么t0是0。λ為交叉梯度項的權(quán)重。

對式(9)數(shù)據(jù)項泰勒展開，并對m求導(dǎo)，令其等于0，可得：

(JTDTDJ+βWTW+λBTB)·δm=-(JTDTD(d(mk-1)-dobs)+βWTW(m-mref)+λBT(t-t0))

(12)

式中：J是數(shù)據(jù)靈敏度矩陣；B是交叉梯度靈敏度矩陣；k是迭代次數(shù)。

令：

H=JTDTDJ+βWTW+λBTB

(13)

g=JTDTD(d(mk-1)-dobs)+βWTW(m-mref)+λBT(t-t0)

(14)

式(12)可寫為：

H·δm=-g

(15)

要求解該式，如果用直接方法求解需要大量計算時間和儲存空間，本文使用Dembo et al.[11]提出的不精確高斯牛頓法(IGN),解式(15)的具體算法為：

(a)在計算g時，計算出一個高精度的g(不用直接計算存儲J，而是計算J和任意向量x的乘積Jx和JTx，詳見下節(jié))，保證下降方向正確。

(b)用預(yù)處理共軛梯度法(PCG)求解式(15)。

(c)在使用PCG求解時，以一個較低的精度計算J以得到H。

(d)在PCG中進行不超過五次的迭代，得到一個近似的模型改正量(每次迭代都要重新生成J)。

(e)得到的模型改正量δm可以對反演模型進行更新：

mi+1=mi+∝δm

(16)

式中：∝是線性搜索參數(shù)，用來控制模型修正量的大小。

以此可以求得極化率或其他參數(shù)約束電阻率反演的模型更新向量。

2.2 極化率聯(lián)合反演方法

筆者在極化率單獨反演中，利用等效電阻率公式將兩次電阻率反演結(jié)果(等效電阻率反演結(jié)果和真電阻率反演結(jié)果)轉(zhuǎn)換成極化率，由于此方法沒有直接的極化率目標函數(shù)，前文描述中使用的聯(lián)合反演理論無法直接算出極化率聯(lián)合反演結(jié)果。筆者提出了一種間接實現(xiàn)極化率聯(lián)合反演的方法。在不使用極化率直接反演的情況下成功實現(xiàn)了電阻率極化率的聯(lián)合反演。具體過程為：

在已知極化率單獨反演結(jié)果的基礎(chǔ)上，進行由極化率約束的真電阻率反演，目標函數(shù)為：

(17)

用上一步得到的真電阻率反演結(jié)果ρ來約束等效電阻率的反演，目標函數(shù)為：

(18)

構(gòu)造一個新的的目標函數(shù)：

(19)

通過等效電阻率公式：

(20)

求得的結(jié)果η即為電阻率約束的極化率聯(lián)合反演結(jié)果。

2.3 目標函數(shù)梯度項的計算

(21)

對上式兩邊對于電導(dǎo)率σ求導(dǎo)，得到：

(22)

由第二節(jié)可知，正演方程組可寫為：

KU=q

(23)

兩邊同時對電導(dǎo)率求導(dǎo)可得：

(24)

把式(24)代入式(22)中，得到：

(25)

(26)

令：

(27)

則：

(28)

把式(28)右端看成正演方程組中的源項，按照正演的方法得到向量φ，代入式(26)，求得偏導(dǎo)數(shù)矩陣與任意向量x的乘積Jx。

偏導(dǎo)數(shù)矩陣的轉(zhuǎn)制JT與任意向量y的乘積JTy的求解方法相似，在此不再贅述。

帶有J的梯度項可以通過這種方式計算，以下為帶有B和t的梯度項的計算方式。

交叉梯度函數(shù)t定義式為：

t(x,y,z)=▽mp(x,y,z)×▽mr(x,y,z)

=(tx,ty,tz)

(29)

彭淼[12]將式中tx，ty，tz進行了離散化，具體形式為：

(30)

(31)

(32)

式中：m的下角標p，r代表極化率和電阻率，上角標代表網(wǎng)格離散化過程中網(wǎng)格位置(圖1)。

圖1 交叉梯度項網(wǎng)格離散示意圖Fig.1 Grid discrete diagram of cross gradient term

2.4 交叉梯度網(wǎng)格離散化邊界的處理方法

筆者采用中心差分方式對交叉梯度項進行離散化。閆政文[8]推導(dǎo)了詳細的三維交叉梯度函數(shù)公式，提出了在中心差分網(wǎng)格情況下的交叉梯度函數(shù)公式，將地下網(wǎng)格分為4類：內(nèi)部網(wǎng)格、邊界面網(wǎng)格、邊界線網(wǎng)格和邊界點網(wǎng)格。每一種網(wǎng)格的交叉梯度項公式都不同。所有梯度項一共是24個公式。筆者提出了一種將地下網(wǎng)格增加假想層的方法。只需要使用上文提到的3個公式及其導(dǎo)數(shù)即可表示出所有位置網(wǎng)格的交叉梯度項。具體方法是在反演網(wǎng)格外面假想地增加一層。如圖2所示(透明網(wǎng)格為假想層)，前述反演網(wǎng)格變成了內(nèi)部網(wǎng)格，可以用內(nèi)部網(wǎng)格公式進行計算。右側(cè)三個圖分別為tx，ty，tz的示意圖。假想層的物性參數(shù)與反演初值相同，保證交叉梯度項在邊界附近的穩(wěn)定性(均為0)，而又不影響感興趣區(qū)的交叉梯度結(jié)果，提高了編程效率。

圖2 交叉梯度項假想層示意圖Fig.2 Imaginary layer diagram of cross gradient term

2.5 交叉梯度項計算可靠性驗證

為了驗證本程序提出的計算交叉梯度項及其導(dǎo)數(shù)的算法的可靠性，本文設(shè)計了如下理論模型：電阻率模型隨深度漸變，極化率模型是一個6 m*4 m*4 m的高極化異常體埋藏在背景極化率為0的區(qū)域(圖3，圖4)，均為yz剖面。

圖3 電阻率理論模型(可靠性驗證)Fig.3 Theoretical model of resistivity (test)

圖4 極化率理論模型(可靠性驗證)Fig.4 Theoretical model of polarizability(test)

用理論模型計算交叉梯度的tx，ty項，得到的結(jié)果(圖5，圖6)。

圖5 交叉梯度項tx在YZ方向上的剖面圖Fig.5 Cross gradient term tx in YZ direction

圖6 交叉梯度項ty在XZ方向上的剖面圖Fig.6 Cross gradient term ty in XZ direction

由于X方向和Y方向網(wǎng)格劃分不同，導(dǎo)致XZ方向上的交叉梯度值看起來更寬。另外由于電阻率理論模型在XY方向上均勻，交叉梯度項tz的值為0。

觀察上述交叉梯度結(jié)果可知，在物性變化的邊界交叉梯度值不為0，在其他區(qū)域交叉梯度項均為0。由此可證明交叉梯度項t計算的可靠性。

3 模型試算

本文用理論模型進行試算，采用井下三極裝置，網(wǎng)格剖分情況為48 m*24 m*22 m。異常體情況如圖所示(圖7，圖8)，異常體A(左側(cè))和異常體B(右側(cè))大小形狀相同。電阻率相同，右側(cè)異常體B有極化率異常，異常體A沒有極化率異常。

圖7 電阻率正演理論模型Fig.7 Theoretical model of resistivity forward modeling

圖8 極化率正演理論模型Fig.8 Forward model of polarizability

圖9為電阻率單獨反演結(jié)果，圖10為極化率單獨反演結(jié)果，圖11為電阻率聯(lián)合反演結(jié)果，圖12為極化率聯(lián)合反演結(jié)果。可以看出，無論是電阻率結(jié)果還是極化率結(jié)果，在兩種異常同結(jié)構(gòu)的情況下，電阻率/極化率聯(lián)合反演相較于單獨反演在數(shù)值恢復(fù)上效果更好，在兩種異常結(jié)構(gòu)不同的情況下，交叉梯度聯(lián)合反演的效果和單獨反演效果相同。這是因為交叉梯度值只在兩種物性同時變化時才有異常，在任意一種物性均勻的情況下均為0，這時交叉梯度項對反演結(jié)果沒有影響。

圖9 電阻率單獨反演結(jié)果Fig.9 Single inversion results of resistivity

圖10 極化率單獨反演結(jié)果Fig.10 Single inversion results of polarizability

圖11 電阻率聯(lián)合反演結(jié)果Fig.11 Joint inversion results of resistivity

圖12 極化率聯(lián)合反演結(jié)果Fig.12 Joint inversion results of polarizability

4 結(jié)論

(1)交叉梯度約束反演比單獨反演對同結(jié)構(gòu)異常體的異常值恢復(fù)效果更好，符合預(yù)期成果。

(2)不同結(jié)構(gòu)異常體的聯(lián)合反演效果與單獨反演相同，交叉梯度項只對同結(jié)構(gòu)的異常體起作用。