基于有限差分的2.5維直流電阻率法的改進(jìn)

蔣 亮， 鄧居智， 陳 輝， 祝福榮， 孟小杰

(東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，江西 南昌 330013)

直流電阻率法根據(jù)地下巖石的導(dǎo)電性差異，通過(guò)記錄觀測(cè)和研究激發(fā)人工電流場(chǎng)所得到的數(shù)據(jù)規(guī)律來(lái)反演地下結(jié)構(gòu)。在高精度正演計(jì)算中，應(yīng)用有限差分法進(jìn)行計(jì)算時(shí)，受傅里葉系數(shù)和邊界條件的影響較大，會(huì)導(dǎo)致計(jì)算精度不夠高(張華等，2012)。Dey等 (1979)開(kāi)發(fā)的一款2.5D模型算法，在這種方法中，地質(zhì)描述為2D，點(diǎn)電源為三維源，而不必明確3D模型，這是基于在很多地質(zhì)情況下，可以假設(shè)地下的導(dǎo)電性在一維結(jié)構(gòu)是不變的這一理論(吳信民等，2013)。根據(jù)這種模型理論，在應(yīng)用有限差分的基礎(chǔ)上，對(duì)公式中的傅里葉系數(shù)進(jìn)行改進(jìn);采用邊界校正，對(duì)其邊界條件進(jìn)行優(yōu)化，使正演效果收邊界條件的影響更小。

1 2.5維直流電阻率的理論推導(dǎo)

設(shè)點(diǎn)電源為水平地面上的A點(diǎn)，電流強(qiáng)度為I，用j表示電流的密度矢量，E為電場(chǎng)強(qiáng)度，u表示電位，ρ為巖石的電阻率，σ=1/ρ表示介質(zhì)電導(dǎo)率，由穩(wěn)定電流場(chǎng)的特性可以得到:

根據(jù)(1)(2)式得:

又由電荷守恒定律導(dǎo)出與j和q相關(guān)的連續(xù)性方程:

將(3)式代入(4)式得到U與q的分布式:

引入狄拉克函數(shù)δ，導(dǎo)出電位u與電流強(qiáng)度I的方程，將點(diǎn)電荷e放在(x0，y0，z0)處，最后由δ和 q得出:

通過(guò)(5)(6)可以得到在合適的邊界條件下，電導(dǎo)公式可以寫(xiě)成(楊金鳳，2012):

式(7)中φ是電勢(shì)場(chǎng)，r+r－是電流源的正負(fù)極位置。

假設(shè)(?/?y)σ(x，y，z)=0，那么方程(7)可以變?yōu)?

這是2.5D的電導(dǎo)率公式。

為了進(jìn)一步求解需要進(jìn)行傅里葉變換，這里傅里葉余弦變換的正向和逆向如下:

其中，ky是波數(shù)。

將上述變換代入(8)式可以得到:

將上式用矩陣式表示為:

D和G分別表示為散度和梯度的二維矩陣，S(σ)是電導(dǎo)率，~u是轉(zhuǎn)換后的電勢(shì)矢量，A(σ，k2y)是前移算子矩陣，q是含有正負(fù)電流源位置的矢量。

當(dāng)給定一個(gè)電導(dǎo)模型和波數(shù)，可以將上式改寫(xiě)為:

將公式(13)和(10)聯(lián)合求解變形得到關(guān)于u的方程，這正是所要求的:

其中n是波數(shù)，gn是每個(gè)波數(shù)的總和(Moucha，2003)。

2 算法改進(jìn)

2.1 傅里葉系數(shù)k和g的優(yōu)化

要以最小的波數(shù)得到最準(zhǔn)確的求解，需要一種優(yōu)化技術(shù)，選擇合適的k和g，優(yōu)化k和g空間的點(diǎn)源，使其解決u，k和g的優(yōu)化過(guò)程如下:

在一個(gè)均勻半空間的偶極源中可以得到電位:

r所處的坐標(biāo)是(x，z)，rim是該源的邊界位置，所以將(15)式帶入到(9)式變?yōu)?

積分求解得:

K0是Bessel函數(shù)零階修正的第二種，通過(guò)上式與(14)式進(jìn)行離散反變換可得:

令Mg=v，g是包含n的向量，v是包含vi值的向量，可得到如下方程:

為了讓誤差達(dá)到最小，選擇k和g使目標(biāo)函數(shù)最小化，可以得出

I是一個(gè)向量，(20)式與k和g相關(guān)，而且與g成線性關(guān)系，M的出現(xiàn)是為了應(yīng)對(duì)r的多變性，利用最小二乘法得到:

聯(lián)合(19)式和(20)式得:

上式是一個(gè)與k相關(guān)的非線性函數(shù)。為了進(jìn)一步優(yōu)化，提出一個(gè)小的阻尼因子，使上式得到進(jìn)一步改進(jìn):

v(ki)是線性化的起點(diǎn)，JT是靈敏度矩陣，K是包含k值得向量，W是單位矩陣是對(duì)模型的加權(quán)，β是盡可能小的阻尼參數(shù)。對(duì)上式進(jìn)行微分，得到關(guān)于K的第i次迭代方程:

定義一個(gè)新的k向量，則第i+1次迭代為:

α是為了減少(23)式迭代次數(shù)而設(shè)立的參數(shù)，當(dāng)(23)式經(jīng)過(guò)迭代到達(dá)理想值，就用該值來(lái)計(jì)算(20)式就可以得到 g(Pidlisecky et al.，2008)。

2.2 邊界條件的優(yōu)化

通過(guò)上式計(jì)算出的電壓值會(huì)存在兩個(gè)問(wèn)題，①邊界校正可以有效的解決邊界效應(yīng)，但要減少相關(guān)錯(cuò)誤，需要對(duì)模型的邊界網(wǎng)格進(jìn)行填充。②奇點(diǎn)存在于源位置造成電勢(shì)衰減迅速，故需要使這個(gè)位置得到精細(xì)的離散。為了解決上述問(wèn)題，應(yīng)用一個(gè)校正的源項(xiàng)q，使其降低BC的影響，并提出一個(gè)類(lèi)似的校正因子實(shí)現(xiàn)減少源奇異的效果，并提高邊界條件(Pidlisecky et al.，2008)。所以將(13)式代入(14)式得到:

將其變形可以得到:

σH是穩(wěn)定的電場(chǎng)的電導(dǎo)率，σT是所要求的真電導(dǎo)率，計(jì)算，為了量化誤差，使 qcorr=L－1·uH，向量 qcorr代替 q在(15)中得到:

3 模型算例

這里就溫納裝置建立模型進(jìn)行算法驗(yàn)證和正演對(duì)比，溫納裝置在測(cè)深過(guò)程中保持三等距，AM=MN=NB(李金銘，2005)。在溫納裝置中，可以測(cè)得最佳深度是AB/6，所得測(cè)點(diǎn)一共有(n－1)(n－2)/6，其中n是電極數(shù)(姚文斌，1989)。故構(gòu)建地下模型并將其剖分成80×30的網(wǎng)格，測(cè)區(qū)的剖分單元為0.5，分別建立2層均勻介質(zhì)模型，高阻模型和低阻模型。為了方便驗(yàn)證算法的準(zhǔn)確性，采用理想化的數(shù)值參數(shù)來(lái)構(gòu)建模型。

3.1 算法驗(yàn)證

圖1 算法驗(yàn)證對(duì)比圖Fig.1 The algorithm validation comparison chart

為了方便驗(yàn)證算法的準(zhǔn)確性，采用理想化的數(shù)值來(lái)構(gòu)建模型，2層均勻介質(zhì)的模型(圖1)，第一層的電導(dǎo)率為σ=1 S/m，第二層的電導(dǎo)率為σ=100 S/m，經(jīng)過(guò)計(jì)算可以得到其地表電壓曲線圖(U0理論曲線，U1只加入有限單元法，U2是在U1基礎(chǔ)上加入傅里葉變換，U3是在U2基礎(chǔ)上加入BC校正)，通過(guò)地表電壓曲線可以看出，U3與U1，U2相比更加精確。最后對(duì)所得結(jié)果進(jìn)行誤差分析，由誤差可以看出U3與U2，U1相比雖然大部分誤差趨于零，但在靠近偶極源位置其誤差更小，這充分說(shuō)明了加入改進(jìn)的傅里葉系數(shù)和邊界條件，可以使計(jì)算精度進(jìn)一步提高。

3.2 低阻模型

構(gòu)建低阻模型(圖2)，模型實(shí)際長(zhǎng)度是30 m長(zhǎng)，10 m深，這里剖分成80×30的網(wǎng)格，正方體剖分，低阻是圖2中的紅色區(qū)域，其大小是2×2，電導(dǎo)率為σ=100 S/m，周?chē)鷧^(qū)域?yàn)棣?1 S/m，經(jīng)正演計(jì)算可以得到以x是長(zhǎng)度，z是深度的視電阻率(圖3)，其中圖3a是以有限差分正常反演的結(jié)果，圖3b是進(jìn)一步加入傅里葉變換所產(chǎn)生的結(jié)果，圖3c是在圖3b的基礎(chǔ)上加入了BC校正的結(jié)果。可以看出，圖3b與圖3a的結(jié)果基本相似略好一點(diǎn)，圖3c與圖3a，圖3b相比層次感更加明顯，反演結(jié)果有一定的提高，更能清晰反映出異常體所在。

圖2 低阻模型Fig.2 Low resistivity model

3.3 高阻模型

構(gòu)建高阻模型(其模型與低阻模型構(gòu)建一樣)，高阻異常其電導(dǎo)率為 σ=0.01 S/m，周?chē)鷧^(qū)域?yàn)棣?1 S/m，經(jīng)計(jì)算反演可以得到其以x是長(zhǎng)度，z是深度，的視電阻率ρ(Ω·m)(圖4)其中圖4a是正常反演的的結(jié)果，圖4b是加入傅里葉變換所產(chǎn)生的結(jié)果，圖4c在圖4b的基礎(chǔ)上加入了BC校正的結(jié)果。同樣圖4c與圖4a，圖4b相比層次感更加明顯，反演結(jié)果有一定的提高，更能反映出異常體的現(xiàn)狀。

4 結(jié)語(yǔ)

使用優(yōu)化的傅里葉系數(shù)，可以有效提高數(shù)值模擬的計(jì)算精度;通過(guò)對(duì)邊界進(jìn)行優(yōu)化，在正演過(guò)程可以使異常區(qū)域更加明顯，結(jié)果受邊界條件的影響更小。正演是反演的基礎(chǔ)，希望上述改進(jìn)可以對(duì)今 后反演工作有所幫助。

圖3 低阻模型對(duì)比圖Fig.3 Low resistance model contrast figure

圖4 高阻模型對(duì)比圖Fig.4 High resistance model contrast figure

李金銘.2005地電場(chǎng)與電法勘探[M].北京:地質(zhì)出版社:246-251.

吳信民，楊海燕，楊亞新，等.2013.論電法勘探的理論探測(cè)深度[J].東華理工大學(xué)學(xué)報(bào):自然科學(xué)版，36(1):60-64.

楊金鳳.2012.基于有限差分的直流電阻率法三維正演研究[D].南昌:東華理工大學(xué).

姚文斌.1989.電測(cè)深數(shù)值計(jì)算和解釋入門(mén)[M].北京:地震出版社:3-18.

張華，共育齡.2012.花崗巖裂隙儲(chǔ)層各向異性高階有限差分?jǐn)?shù)值模擬[J].東華理工大學(xué)學(xué)報(bào):自然科學(xué)版，35(4):416-421.

Dey A，Morrison H F.1979.Resistivity modeling for arbitrarily shaped two-dimensional structures[J].Geophysical Prospecting ，27:106-136.

Pidlisecky A，Knight R.2008.FW2_5D:A MATLAB 2.5D-D electrical resistivity modeling code[J].Computers ＆ Geosciences，34:1645-1654.

Moucha R.2003.Multigrid methods for forward and inverse resistivity problems in geophysics[M].Canada:National Library of Canada:1-63.