一種復(fù)雜分層介質(zhì)中電磁探測(cè)的正演算法研究

葛勤革 韓逍菲 潘威炎 羅 強(qiáng)

(1.海司信息化部，北京 100036；2.中國(guó)電波傳播研究所，山東 青島 266107；3.中國(guó)人民解放軍91919部隊(duì)，湖北 黃岡 438000)

引 言

近年來(lái)，SLF/ELF頻段的電磁波用于地球物理勘探正被廣泛重視。在CSEM的正演計(jì)算中最早及早期的分析模型都是均勻成層模型，電磁波在這些模型中傳播時(shí)在接收點(diǎn)處所產(chǎn)生的場(chǎng)值已有解析表達(dá)式，并有許多算法可以對(duì)這些場(chǎng)值進(jìn)行快速準(zhǔn)確的計(jì)算。但是，實(shí)際的物理模型不是嚴(yán)格平面分層的，對(duì)于這些模型沒(méi)有解析表達(dá)式進(jìn)行計(jì)算。近些年來(lái)，隨著計(jì)算機(jī)性能的提高很多人致力于電磁場(chǎng)在不均勻介質(zhì)中傳播的數(shù)值算法的研究。

常用的電磁場(chǎng)的數(shù)值模擬方法可分為以微分方程為基礎(chǔ)的有限差分法、有限元法和以積分方程為基礎(chǔ)的積分方程法、矩量法、邊界元法兩大類(lèi)。對(duì)于我們所討論的SLF/ELF波的傳播問(wèn)題，由于工作頻率很低，實(shí)際電磁場(chǎng)的空間變化很快，而時(shí)間變化很慢，若采用時(shí)域有限差分法則計(jì)算過(guò)程很難收斂并穩(wěn)定，誤差亦難以控制。若單純采用積分方程法，對(duì)于復(fù)雜模型下的計(jì)算也面臨著矩陣求解以及計(jì)算速度方面的一系列問(wèn)題[1]。對(duì)于目前廣泛采用的有限元方法，雖然其在解決不均勻結(jié)構(gòu)問(wèn)題時(shí)與其他的方法相比有很大的優(yōu)越性，但該方法生成網(wǎng)格的成本高，需要進(jìn)行單元剖分、節(jié)點(diǎn)編號(hào)等繁重的前處理過(guò)程，生成矩陣的存儲(chǔ)量大，在實(shí)際的計(jì)算中這些矩陣往往是病態(tài)的，計(jì)算的精度和速度都很不理想。

對(duì)于探測(cè)正演中所涉及的復(fù)雜分層結(jié)構(gòu)中電偶極子激發(fā)的電磁場(chǎng)的數(shù)值計(jì)算問(wèn)題，本文將復(fù)雜分層結(jié)構(gòu)看成是均勻分層結(jié)構(gòu)中含有目標(biāo)體，如果目標(biāo)體的埋深比較深，并且，導(dǎo)電率與背景介質(zhì)的導(dǎo)電率反差不大，則可以將均勻分層中的目標(biāo)體看成對(duì)于均勻分層的“擾動(dòng)”，將僅含均勻分層的初級(jí)場(chǎng)與由目標(biāo)體的“擾動(dòng)”產(chǎn)生的次級(jí)場(chǎng)之和看成含有目標(biāo)體的不均勻分層介質(zhì)中的總場(chǎng)。在上述地質(zhì)條件下，可以忽略目標(biāo)體內(nèi)部場(chǎng)的相互作用，避免了積分方程法在方程求解方面所存在的一系列問(wèn)題。在計(jì)算過(guò)程中，將充分利用在均勻分層介質(zhì)電波傳播方面的研究成果，結(jié)合格林函數(shù)的物理意義，對(duì)總場(chǎng)值進(jìn)行直接計(jì)算，將該算法稱(chēng)為基于“擾動(dòng)”的格林函數(shù)法。文中給出了該方法在水平電偶源激發(fā)，水平n層分層介質(zhì)中含有目標(biāo)體時(shí)，接收點(diǎn)場(chǎng)強(qiáng)的計(jì)算思路。結(jié)合水平半空間中含有目標(biāo)體這一模型分析了該算法的有效性以及需要進(jìn)一步解決的問(wèn)題。

1 算法原理

考慮海底N層均勻水平分層背景介質(zhì)，在最上層沿著x軸方向放置水平電偶極子發(fā)射源，接收點(diǎn)設(shè)在沿著x軸的一維方向上且位置可變。

當(dāng)在背景介質(zhì)的第i層中存在一個(gè)目標(biāo)層時(shí)，如圖1，可以認(rèn)為整個(gè)區(qū)域的電磁場(chǎng)分布將在原來(lái)的電磁場(chǎng)分布基礎(chǔ)上增加一個(gè)“擾動(dòng)”，即

在這里稱(chēng)背景介質(zhì)中所對(duì)應(yīng)的場(chǎng)Ein和Hin為初級(jí)場(chǎng)，由目標(biāo)體的“擾動(dòng)”所產(chǎn)生的場(chǎng)Ese和Hse為次級(jí)場(chǎng)，初級(jí)場(chǎng)與次級(jí)場(chǎng)之和即為接收點(diǎn)處的總場(chǎng)值Etot和Htot.其中，初級(jí)場(chǎng)Ein和Hin滿足Maxwell方程為

式中：σ0為水平均勻分層介質(zhì)中的電導(dǎo)率分布；r0是激勵(lì)源所在的位置。

圖1 含有目標(biāo)體的水平N層分層模型

當(dāng)存在目標(biāo)層后，在整個(gè)區(qū)域內(nèi)總場(chǎng)分布Etot和Htot應(yīng)滿足

將式(1)～(4)代入(5)、(6)可得整個(gè)區(qū)域內(nèi)的次級(jí)場(chǎng)滿足

(9)

進(jìn)一步表示為

(14)

3) 目標(biāo)區(qū)域Ω內(nèi)各個(gè)點(diǎn)處所對(duì)應(yīng)的x、y、z三個(gè)方向偶極子在接收點(diǎn)處產(chǎn)生的場(chǎng)，即式(14)三階方陣中各元素的計(jì)算。

對(duì)于1)、2)兩種情況下場(chǎng)強(qiáng)的計(jì)算已根據(jù)周永祖[3]和葛德彪[4]等人的推導(dǎo)思路進(jìn)行了具體的推導(dǎo)，并針對(duì)本文的算法得到了適合于編程實(shí)現(xiàn)的具體表達(dá)式，限于篇幅原因?qū)⒘砦陌l(fā)表。對(duì)于3)的計(jì)算可根據(jù)King[5]、潘威炎等人[6-7]給出的積分表達(dá)式直接計(jì)算。

2 算法驗(yàn)證

選用背景分層介質(zhì)為最簡(jiǎn)單的均勻半空間模型，通過(guò)在背景介質(zhì)的區(qū)域2中插入目標(biāo)體Ω的方式來(lái)構(gòu)造不均勻分層結(jié)構(gòu)(如圖2左)。這種地質(zhì)模型下接收點(diǎn)的場(chǎng)強(qiáng)沒(méi)有解析公式對(duì)其進(jìn)行計(jì)算。但是考慮到當(dāng)接收點(diǎn)離發(fā)射源比較近，對(duì)于接收點(diǎn)的場(chǎng)強(qiáng)的計(jì)算結(jié)果應(yīng)接近于水平無(wú)限大四層模型(圖2右)下相應(yīng)接收點(diǎn)場(chǎng)強(qiáng)的計(jì)算結(jié)果。故下文將通過(guò)利用格林函數(shù)法計(jì)算出的圖2左模型的接收點(diǎn)的電場(chǎng)x分量與利用直接積分法計(jì)算出的圖2右模型的接收點(diǎn)的電場(chǎng)x分量作對(duì)比的方式來(lái)說(shuō)明格林函數(shù)算法的有效性。

圖2 均勻半空間含目標(biāo)體模型 與水平無(wú)限大四層模型

對(duì)于圖2模型,式(14)矩陣中的格林函數(shù)表示為

(15)

(16)

(17)

式(15)、(16)、(17)中的φ′為接收點(diǎn)偏離目標(biāo)體中第m個(gè)剖分單元中心節(jié)點(diǎn)的水平相位。將φ′減小π/2可得：Gyx、Gyy、Gyz.

ei(γ1d-γ2z)λ2dλ

(18)

將式(18)中相位φ減小π/2得Gzy.

(19)

式中：J0、J1、J2分別為零階、一階和二階貝塞爾函數(shù)。γj=(kj2-λ2)1/2;kj=ω(μ0(ε0εrj+iσj/ω))1/2；j=1，2;對(duì)于M、N的具體表達(dá)式見(jiàn)參考文獻(xiàn)[5]。

取海水的介質(zhì)參數(shù)為：εr0=80，σ0=3.3 s/m；海底巖層：εr1=10，σ1=1 s/m；目標(biāo)層(四層介質(zhì)中的區(qū)域2)的介質(zhì)參數(shù)：εr2=10，σ2根據(jù)需要取值；海底巖層：εr3=10，σ3=1 s/m.各層中的磁導(dǎo)率均為真空中的磁導(dǎo)率μ0.水平偶極子源的工作頻率f=1 Hz,距離海底巖層的高度d=10 m.

1) 目標(biāo)層導(dǎo)電率σ2=0.1 s/m，將圖2左模型中目標(biāo)體的長(zhǎng)寬均取為5 000 m，埋深l1取為600 m，采用本文所提算法對(duì)于圖2左模型中接收點(diǎn)總電場(chǎng)x分量幅值的計(jì)算結(jié)果與采用直接積分法對(duì)于圖2右模型中接收點(diǎn)總電場(chǎng)x分量幅值的計(jì)算結(jié)果對(duì)比如圖3；將目標(biāo)層(區(qū)域2)的埋深l1設(shè)為1 000 m，其他條件不變所得的結(jié)果對(duì)比見(jiàn)圖4所示。

圖3 目標(biāo)體埋深600 m計(jì)算結(jié)果對(duì)比圖

圖4 目標(biāo)體埋深1 000 m計(jì)算結(jié)果對(duì)比圖

從圖3和圖4可以看出：在收發(fā)距較小的情況下，計(jì)算結(jié)果的相對(duì)誤差均小于10%；當(dāng)接收點(diǎn)接近或超過(guò)目標(biāo)體的水平邊緣時(shí)不能再將圖2左看作圖2右的近似模型，故相對(duì)誤差也明顯增大。另外，通過(guò)圖3和圖4的對(duì)比可以看出：當(dāng)埋深較深時(shí)兩種算法的相對(duì)誤差要小的多。這主要是因?yàn)閷?duì)于無(wú)限大分層介質(zhì)，當(dāng)目標(biāo)層埋深較深時(shí)，由于目標(biāo)體的存在所產(chǎn)生的側(cè)面波和反射波主要是在遠(yuǎn)場(chǎng)區(qū)起作用，因此，當(dāng)收發(fā)距在埋深的2～3倍范圍內(nèi)，兩種模型的近似效果較好，相應(yīng)的計(jì)算誤差也較小，收發(fā)距小于3 000 m時(shí)，相對(duì)誤差小于10%，從而對(duì)算法的有效性進(jìn)行了初步驗(yàn)證。

2) 從算法的推導(dǎo)過(guò)程可以看出：總場(chǎng)中的初級(jí)場(chǎng)由背景均勻介質(zhì)決定，可以認(rèn)為是已知的。可反映目標(biāo)體位置及介質(zhì)特性的是次級(jí)場(chǎng)。圖5、圖6為目標(biāo)體的長(zhǎng)寬均取為5 000 m，埋深1 000 m，導(dǎo)電率取不同的值時(shí)接收點(diǎn)處次級(jí)電場(chǎng)x分量實(shí)部、虛部隨發(fā)-收距的變化。通過(guò)圖5、圖6可以看出：在目標(biāo)體所在區(qū)域附近，次級(jí)場(chǎng)對(duì)于總場(chǎng)的作用較強(qiáng)。

圖5 次級(jí)電場(chǎng)x分量實(shí)部

圖6 次級(jí)電場(chǎng)x分量虛部

圖7為長(zhǎng)寬均為5 000 m的目標(biāo)體位于發(fā)射源正下方，垂直埋深1 000 m時(shí)次級(jí)電場(chǎng)x分量幅值(dB)在xoy面內(nèi)的分布圖。圖8為長(zhǎng)寬均為1 500 m的目標(biāo)體，位于發(fā)射源正下方左側(cè)500 m，垂直埋深1 000 m時(shí)次級(jí)場(chǎng)的水平分布圖。

通過(guò)圖7、圖8可以看出隨著目標(biāo)體的大小和埋藏位置的不同，次級(jí)場(chǎng)的主要作用區(qū)域(圖中深紅色等值線的分布區(qū)域)的大小和位置相應(yīng)地發(fā)生了改變，這為反演計(jì)算提供了依據(jù)。

3 結(jié) 論

針對(duì)水平分層介質(zhì)中含有目標(biāo)體，目標(biāo)體埋深較深、導(dǎo)電率與背景層的導(dǎo)電率反差不大的地質(zhì)模型，給出了基于“擾動(dòng)”的格林函數(shù)算法思路，并且選取簡(jiǎn)單的分層模型對(duì)算法的有效性進(jìn)行了驗(yàn)證。通過(guò)驗(yàn)證得出：在該模型下基于“擾動(dòng)”的格林函數(shù)算法對(duì)接收點(diǎn)總場(chǎng)進(jìn)行計(jì)算，在避免進(jìn)行求解積分方程的同時(shí)，計(jì)算精度可接受。由于采用該算法時(shí)主要涉及到含有貝塞爾函數(shù)的積分計(jì)算，在提高計(jì)算速度方面需要進(jìn)一步研究。

[1] 何展翔, 孫衛(wèi)斌, 孔繁恕, 等. 海洋電磁法[J]. 石油地球物理, 2006, 41(4): 451-457.

HE Zhanxiang, SUN Weibin, KONG Fanshu, et al. Marine electromagnetic approach[J]. Oil Geophysical Prospecting, 2006, 41(4): 451-457. (in Chinese)

[2] Chen-To-Tai. Dyadic Green Functions in Electromagnetic Theory[M]. 2nd. New York: The IEEE Press, 1993: 55-59.

[3] 周永祖. 非均勻介質(zhì)中的場(chǎng)與波[M]. 聶在平, 柳清伙, 譯. 北京: 電子工業(yè)出版社, 1992: 48-61.

[4] 葛德彪. 電磁波理論[M]. 西安: 西安電子科技大學(xué)出版社, 1997: 38-48.

[5] KING R W P, OWENS M, WU T T. Lateral Electromagnetic Wave[M]. New York: Springer-verlag Press, 1992: 151-170.

[6] 潘威炎. 長(zhǎng)波超長(zhǎng)波極長(zhǎng)波傳播[M]. 成都: 電子科技大學(xué)出版社, 2004: 165-170.

[7] 李 凱. 分層介質(zhì)中的電磁場(chǎng)和電磁波[M]. 杭州: 浙江大學(xué)出版社, 2010: 53-72.

[8] WEISS C J, CONSTABLE S. Mapping thin resistors and hydrocarbons with marine EM methods, part II: modeling and analysis in 3D[J]. Geophysics, 2006, 71(6): 321-332.

[9] HE Z X, KURT S, YU G, et al. On reservoir boundary detection with marine CSEM[J]. Applied Geophyiscs, 2008, 5(3): 181-188.

[10] 彭懷云, 陶 偉, 潘威炎, 等. 極低頻垂直偶極子在地-電離層中場(chǎng)的數(shù)值積分算法[J]. 電波科學(xué)學(xué)報(bào), 2012, 27(2): 333-338.

PENG Huaiyun, TAO Wei, PAN Weiyan, et al. Numerical integral method for ELF fields excited by vertical electric dipole in asymmetric earth-ionosphere cavity[J]. Chinese Journal of Radio Science, 2012, 27(2): 333-338. (in Chinese)

[11] ZHAO Luanxiao,GENG Jianhua,ZHANG Shengye,et al. 1-D controlled source electromagnetic forward modeling for marine gas hydrates studies[J]. Applied Geophysics,2008,5(2): 121-126.

[12] BADEA E A. Finite-element analysis of controlled source electromagnetic induction using Coulomb gauged potentials[J]. Geophysics, 2001, 66(3): 786-799.