用數(shù)值模式匹配算法高效仿真軸對(duì)稱(chēng)型散射體海洋可控源電磁響應(yīng)?

林藺 焦利光 陳博 康莊莊 馬玉剛 汪宏年

(吉林大學(xué)物理學(xué)院,長(zhǎng)春 130012)

用數(shù)值模式匹配算法高效仿真軸對(duì)稱(chēng)型散射體海洋可控源電磁響應(yīng)?

林藺 焦利光 陳博 康莊莊 馬玉剛 汪宏年?

(吉林大學(xué)物理學(xué)院,長(zhǎng)春 130012)

(2016年11月17日收到;2017年4月10日收到修改稿)

圓盤(pán)、球體以及球冠狀體是地球物理研究中非常重要的一類(lèi)散射類(lèi)型.在海洋環(huán)境中,圓盤(pán)可以用于描述玄武巖基巖以及油氣圈閉構(gòu)造等電阻率異常體,而球冠可以近似描述某些基巖隆起或起伏地形等.這類(lèi)散射體的一個(gè)重要特征是其電阻率空間分布具有軸對(duì)稱(chēng)性.如果能夠針對(duì)這類(lèi)形狀的散射體研究建立一套有效的海洋可控源電磁數(shù)值模擬方法,對(duì)于認(rèn)識(shí)復(fù)雜地層條件下海洋電磁響應(yīng)的變化特征、研究建立相關(guān)的資料處理和解釋方法具有非常重要的意義.本文根據(jù)電導(dǎo)率軸對(duì)稱(chēng)分布特征,設(shè)法用一個(gè)或多個(gè)不同半徑、不同厚度的水平同心圓盤(pán)逼近這類(lèi)軸對(duì)稱(chēng)電導(dǎo)率散射體,并將這些同心圓盤(pán)與海洋環(huán)境中的空氣、海水、沉積層和基巖等背景介質(zhì)結(jié)合,形成一個(gè)在水平方向電導(dǎo)率具有軸對(duì)稱(chēng)分布、在垂直方向又具有分層特征的水平層狀非均質(zhì)模型.在此基礎(chǔ)上,應(yīng)用數(shù)值模式匹配法研究水平電偶極子天線電磁場(chǎng)的數(shù)值模擬方法,給出位于對(duì)稱(chēng)軸上的水平發(fā)射天線電磁場(chǎng)在層狀非均質(zhì)地層中的半解析解,建立海洋可控源電磁響應(yīng)高效算法.最后通過(guò)數(shù)值模擬結(jié)果對(duì)該算法進(jìn)行檢驗(yàn)并考察海洋可控源三維電磁響應(yīng)特征.

海洋可控源電磁,軸對(duì)稱(chēng)散射體,數(shù)值模式匹配法,半解析解

1 引 言

自20世紀(jì)80年代以來(lái),海洋可控源電磁法(marine controlled source electromagnetic system, MCSEM)已成為研究海洋地形及海底構(gòu)造等重要的地球物理方法[1].由于MCSEM對(duì)海底高阻目標(biāo)體的有效探測(cè)能力,海洋可控源電磁法已逐漸發(fā)展成為海洋油氣資源勘探的有效手段[2-4],關(guān)于海洋可控源電磁數(shù)值模擬與響應(yīng)特征的研究已成為當(dāng)前重要的一項(xiàng)研究課題[5,6].經(jīng)過(guò)幾十年的發(fā)展,一維水平層狀地層中電磁響應(yīng)的解析法[7],二維和三維非均質(zhì)地層中的有限元法[8,9]、有限差分法[10,11]、有限體積法[12]和積分方程法[13,14]等數(shù)值模擬技術(shù)在海洋可控源電磁響應(yīng)研究中均得到研究與應(yīng)用.從理論上說(shuō),三維有限元和三維有限體積法為解決任意復(fù)雜海洋地電條件下的可控源電磁響應(yīng)數(shù)值模擬提供了有力手段,但這類(lèi)數(shù)值模擬的效率往往較低,對(duì)網(wǎng)格劃分也要求嚴(yán)格,特別是對(duì)于高頻、大源距(6000 m以上)情況,由于其電磁強(qiáng)度往往很弱,導(dǎo)致計(jì)算精度明顯下降[15].三維積分方程法是求解局部散射體電磁響應(yīng)比較理想的選擇,但由于涉及滿(mǎn)元矩陣方程的求解,當(dāng)散射體較大時(shí)需要求解的未知參數(shù)很多,給數(shù)值結(jié)果的準(zhǔn)確計(jì)算帶來(lái)極大困難.實(shí)際上,海洋可控源往往采用多頻多源距測(cè)量方式,其在復(fù)雜情況下的精確模擬仍然面臨著嚴(yán)峻挑戰(zhàn)[11],特別是當(dāng)模型中存在起伏地形和起伏邊界時(shí),由于起伏邊界本身的不確定性,大大增加了可控源電磁響應(yīng)的計(jì)算難度.

近期的相關(guān)研究表明,海洋電磁勘探中遇到的典型高阻探測(cè)目標(biāo)是玄武巖、碳酸鹽巖以及油氣圈閉,這些目標(biāo)體可以近似表示為水平盤(pán)狀散射體[5,7].同樣地,某些基巖隆起或海底起伏也可以用球冠狀散射體加以描述.由于水平圓盤(pán)和球冠等散射體的電導(dǎo)率空間分布具有軸對(duì)稱(chēng)性,完全可以用一個(gè)或多個(gè)半徑與厚度不同的水平同心圓盤(pán)加以逼近,如果將這些同心圓盤(pán)與海洋環(huán)境中的空氣、海水、沉積層與基巖等背景介質(zhì)結(jié)合在一起,將形成一個(gè)在水平方向具有軸對(duì)稱(chēng)分布、在垂直方向又具有分層特征的水平層狀非均質(zhì)電導(dǎo)率模型.本文以這種水平層狀非均質(zhì)地層模型為基礎(chǔ),設(shè)法將井中地球物理中廣泛使用的數(shù)值模式匹配算法(NMM)[16-22]推廣應(yīng)用到海洋可控源電磁響應(yīng)數(shù)值模擬中,建立海洋可控源三維電磁響應(yīng)的一種新型高效算法,給出位于對(duì)稱(chēng)軸上水平電偶極子電磁場(chǎng)的半解析解,并用這種算法研究考察水平圓盤(pán)、球形電阻率散射體在基巖或海底具有球冠狀起伏表面時(shí)海洋可控源電磁響應(yīng)的變化特征.

2 基本原理

2.1 海洋地電模型與等效水平層狀非均質(zhì)地層

圖1(a)是本文要首先研究考察的一種海洋可控源地電模型,模型中包含有空氣層、海水、沉積層、基巖以及水平圓盤(pán)散射體和球冠型基巖隆起,其中,空氣和海水電導(dǎo)率分別用σair和σsea表示,沉積層和基巖電導(dǎo)率分別為σsed和σbs,海水和沉積層厚度分別為hsea和hsed.沉積層中水平圓盤(pán)散射體的半徑和厚度分別為adsk和hdsk,其頂面離海底距離為ddsk;而球冠狀基巖隆起的高度為hcrn,其對(duì)應(yīng)的球體半徑為acrn.圖1(b)是由空氣、海水等組成的等效水平層狀非均質(zhì)模擬,水平層界面位置用dn(n=1,2,···,N+5)表示.其中,圓盤(pán)所在地層是非均質(zhì)的,其電導(dǎo)率在柱坐標(biāo)系下可以表示為如下的階躍函數(shù)[16,19]：

球冠狀基巖隆起部分被劃分成N個(gè)等厚度的薄圓盤(pán),各圓盤(pán)水平界面位置為d4+j= hsed-hcrn[1-(j-1)/N](j=1,2,···,N),各圓盤(pán)的半徑根據(jù)球冠邊界形狀從上到下按公式逐漸增加：

圖1 高阻圓盤(pán)散射體和球冠型基巖隆起形成的具有軸對(duì)稱(chēng)電導(dǎo)率分布的海洋地電模型(a)及其等價(jià)水平層狀非均質(zhì)地層(b)Fig.1.M arine geoelectricmodelwith axisymm etric conductivity d istribu tion includ ing both resistive d isk and spherical cap type bedrock up lift(a)and its equivalentmodel consisting of horizontal layered inhom ogeneous beds(b).

而每個(gè)圓盤(pán)所在的水平層狀地層也是非均質(zhì)的,其電導(dǎo)率也具有如下的階躍函數(shù)形式：

而其他地層電導(dǎo)率是常數(shù),可以表示為

因此,整個(gè)模型上電導(dǎo)率是軸對(duì)稱(chēng)分片常數(shù)函數(shù),其空間分布為

2.2 M axw ell方程的分解

海洋可控源電磁響應(yīng)的正演模擬實(shí)質(zhì)上是求解如下M axwell方程[12]：

其中,σ=σ(ρ,z)是軸對(duì)稱(chēng)電導(dǎo)率函數(shù),磁導(dǎo)率μ假定是常數(shù),ω為發(fā)射信號(hào)的角頻率且時(shí)間變化關(guān)系假定為eiωt.由于發(fā)射源工作頻率很低,忽略了位移電流.為簡(jiǎn)單起見(jiàn),假定發(fā)射源Js是一個(gè)單位強(qiáng)度的水平電偶極子發(fā)射天線且位于對(duì)稱(chēng)軸z上.利用直角坐標(biāo)系與柱坐標(biāo)系下單位向量間的變換公式和歐拉公式,可以將單位強(qiáng)度水平電流源展開(kāi)成柱坐標(biāo)系下的Fourier級(jí)數(shù)形式[16]：

其中,

是柱坐標(biāo)系下單位水平電偶極子源的k階諧分量, ρs是足夠小的常數(shù)(例如10-4m),保證在柱坐標(biāo)軸附近電磁場(chǎng)仍然有界[16].

同樣地,將方程(5)中的電磁場(chǎng)E = (σEρ,Eθ,Ez)T和H=(Hρ,Hθ,Hz)T展開(kāi)成類(lèi)似形式的Fourier級(jí)數(shù)：

將(6)和(7)式代入(4)式,并對(duì)方程進(jìn)行分解,得到水平磁場(chǎng)的k階諧分量HS,k=(Hρ,k,Hθ,k)T滿(mǎn)足的方程[17,18,22]：

而右側(cè)項(xiàng)是柱坐標(biāo)系中的二維水平向量,其表達(dá)式為

為保證(8)式的解在z軸附近仍然有意義,引入如下第二類(lèi)齊次邊界條件[16,20]：

這里的ρMN是足夠小的數(shù)(例如10-5m).由于存在趨膚效應(yīng),遠(yuǎn)處的電磁場(chǎng)快速衰減,所以在足夠大外邊界ρMX(例如105m)上,引入截?cái)噙吔鐥l件：

此外,水平電場(chǎng)分量ES,k=(σEρ,k,Eθ,k)T與水平磁場(chǎng)分量間具有如下關(guān)系[17]：

垂直電磁場(chǎng)分量可通過(guò)如下方程得到[16,17]：

2.3 海洋可控源電磁響應(yīng)的半解析解

附錄A中,通過(guò)數(shù)值模式匹配算法給出了完全柱狀介質(zhì)中方程(8)的求解過(guò)程以及相應(yīng)的半解析解,同時(shí)通過(guò)方程(12)和(13)給出了計(jì)算其他電磁分量的方法與相應(yīng)的半解析解,這些不含水平地層界面的半解析解稱(chēng)之為電磁場(chǎng)的基本解.為了得到圖1(b)中水平層狀非均質(zhì)地層中的電磁場(chǎng),需要分析海水中k階水平磁場(chǎng)的基本解(入射波)的傳播過(guò)程.當(dāng)入射波傳播到界面d1和d2時(shí),將產(chǎn)生反射和透射,海水中界面d1附近的總場(chǎng)是入射場(chǎng)和上下界面反射場(chǎng)的疊加,因此可表示為[17,18,22]

等于下行波在界面d2上的反射,即

求解方程(16a)和(16b)可以得到

將(17)式代入(15)式中并經(jīng)適當(dāng)整理,得到發(fā)射源所在的海水中水平磁場(chǎng)k階諧分量的半解析解：

其中,

利用各個(gè)地層界面上電磁場(chǎng)的連續(xù)性條件,可以推導(dǎo)出地層n中的水平磁場(chǎng)k階諧分量的半解析表達(dá)式[16,17]：

利用水平磁場(chǎng)分量的半解析解(19),(21)以及(12)式和附錄A中的(A 10)式,可以得到各個(gè)地層上水平電場(chǎng)k階諧分量的表達(dá)式[16,17].其中,發(fā)射源所在的海水中的電場(chǎng)可表示為

而在其他地層中電場(chǎng)具有如下類(lèi)似的表達(dá)形式：

其中,

是地層n中水平電場(chǎng)k階諧分量的傳播矩陣.同樣地,可以推導(dǎo)出電磁場(chǎng)垂直分量的半解析解[16,17],其中,在發(fā)射源所在的海層中：

而在其他地層n中：

將上面的關(guān)于電磁場(chǎng)水平分量與垂直分量結(jié)合起來(lái),最后得到各個(gè)地層n中電磁場(chǎng)k階諧分量的半解析表達(dá)式：

其中,上式右端出現(xiàn)的上標(biāo)“±”分別表示接收點(diǎn)位于發(fā)射源的上部或下部.

進(jìn)一步,將方程(29)和(30)代入方程(7)就可以計(jì)算出單位水平電流源在各個(gè)地層中任意位置產(chǎn)生的電磁場(chǎng),換言之,通過(guò)解決兩次二維問(wèn)題,就可以得到層狀非均質(zhì)地層中水平電流源產(chǎn)生的電磁場(chǎng),從而大幅提高了數(shù)值模擬的計(jì)算效率.

3 數(shù)值模擬結(jié)果

本節(jié)首先利用圖1中地電模型,分別設(shè)計(jì)出水平層狀地層和含有圓盤(pán)的水平層狀地層兩個(gè)不同的簡(jiǎn)化地電模型,并用傳輸線法(TLM)[7]以及三維有限體積法(3D FV)[12]分別計(jì)算這兩個(gè)地電模型上可控源電磁響應(yīng),通過(guò)與上述NMM的數(shù)值結(jié)果進(jìn)行對(duì)比,檢驗(yàn)NMM的計(jì)算精度與效率.然后,通過(guò)NMM研究高阻圓盤(pán)與基巖隆起地電模型上海洋可控源電磁響應(yīng),在本節(jié)的最后部分,利用NMM進(jìn)一步考察球形散射體與基巖凹陷以及球形散射體、基巖凹陷和海底隆起這兩個(gè)地電模型的海洋電磁響應(yīng).

圖2 水平層狀海洋地電模型NMM和TLM法計(jì)算得到的主測(cè)線上電磁水平分量(Ex和Hy)的振幅與相位曲線對(duì)比(a)Ex振幅曲線;(b)Ex相位曲線;(c)Hy振幅曲線;(d)Hy相位曲線Fig.2.Com parison of the in line magnitudes(M VO)and phases(PVO)versus off set of electrom agnetic(EM) horizontal com ponents(Ex and Hy)ob tained by two d iff erent algorithm s of NMM and TLM in an horizontal layered m arine geoelectric model:(a)MVO of Ex;(b)PVO of Ex;(c)MVO of Hy;(d)PVO of Hy.

在整個(gè)數(shù)值模擬過(guò)程中,徑向求解區(qū)間的變化范圍為ρMN=10-4m和ρMX=5×104m,總節(jié)點(diǎn)數(shù)M+1=181,且所有節(jié)點(diǎn)ρα(α= 1,2,···,M+1)均按對(duì)數(shù)等間距增加.空氣層、海水、沉積層電導(dǎo)率分別為σair=10-6,σsea=3.33和σsed=0.667 S/m;基巖和圓盤(pán)(或球體)的電導(dǎo)率分別為σbs=0.05和σdsk=0.01 S/m;海水層厚度為hsea=1000m.單位強(qiáng)度水平電流源位于z軸上且在海底上方50m處即b=-50m,發(fā)射源工作頻率假定為0.25和0.75 Hz,而接收器均勻分布在海床面上,水平間距為200 m.收發(fā)距的變化范圍為-6400m到6400m.

3.1 算法檢驗(yàn)

首先,假定圖1(a)中圓盤(pán)厚度和頂部離海底距離的分別是hdsk=150 m和ddsk=925 m,且沉積層厚度為hsed=2000 m.為了進(jìn)行TLM與NMM數(shù)值結(jié)果的對(duì)比,進(jìn)一步假定圖1(a)中圓盤(pán)的半徑adsk趨近于無(wú)限大,并且基巖與沉積層電導(dǎo)率相同,即σsed= σbs=0.667 S/m,這時(shí)圖1(a)簡(jiǎn)化為由空氣、海水、沉積層、高阻薄層以及沉積層組成的5層水平層狀模型,因此可以用TLM進(jìn)行正演計(jì)算,圖2是該簡(jiǎn)化模型上由TLM與NMM計(jì)算得到數(shù)值模擬結(jié)果的對(duì)比.圖2(a)和圖2(c)分別是主測(cè)線上電磁場(chǎng)水平分量Ex,Hy的振幅與偏離距(magnitude versus off set,MVO)曲線,圖2(b)和圖2(d)是主測(cè)線上Ex,Hy的相位與偏離距(phase versus off set,PVO)曲線,其中實(shí)線是NMM數(shù)值結(jié)果,而離散符號(hào)是TLM數(shù)值結(jié)果.從圖2可以看出,兩種方法計(jì)算得到的電磁場(chǎng)數(shù)值結(jié)果均符合得很好.電磁場(chǎng)振幅的最大相對(duì)誤差僅為2.65%,NMM的計(jì)算精度達(dá)到了TLM的水平.

圖3 含圓盤(pán)狀散射體的水平層狀海洋地電模型中由NMM和3D FV計(jì)算得到的主測(cè)線上電磁場(chǎng)水平分量(Ex和Hy)振幅與相位曲線對(duì)比 (a)Ex振幅曲線;(b)Ex相位曲線;(c)Hy振幅曲線;(d)Hy相位曲線Fig.3.Com parison of the in line MVO and PVO of EM horizontal com ponents(Ex and Hy)obtained by two d iff erent algorithm s of NMM and 3D FV in the horizontal layered m arine geoelectric model with a resistive d isk: (a)MVO of Ex;(b)PVO of Ex;(c)M VO of Hy;(d)PVO of Hy.

為了進(jìn)一步考察非均質(zhì)層狀海洋地電模型中NMM算法的計(jì)算精度,假定圖1(a)中圓盤(pán)半徑adsk=2000m,基巖與沉積層電導(dǎo)率仍然相同,形成一個(gè)具有不同厚度的5層水平層狀非均質(zhì)地電模型.圖3是該模型分別用NMM法與3D FV法得到的主測(cè)線上Ex,Hy的MVO與PVO曲線對(duì)比結(jié)果,其中,實(shí)線與虛線是NMM計(jì)算結(jié)果,而離散的符號(hào)是3D FV數(shù)值結(jié)果.兩種不同方法計(jì)算得到的數(shù)值結(jié)果同樣符合得較好.圖3中兩個(gè)頻率、65個(gè)接收點(diǎn)的數(shù)值結(jié)果,在HP Z820(Intel?CPU E5-2697 V 2 2.7GHz,256GB RAM)工作站上,NMM和3D FV所用CPU時(shí)間分別為5m in 57 s和65m in 3 s, NMM比3D FV所用時(shí)間少10倍以上,且NMM只占用0.536 Gb內(nèi)存而3D FV需要的內(nèi)存則達(dá)到35 Gb.因此,NMM在普通的PC機(jī)就可以運(yùn)行,從而大大降低了計(jì)算成本.

3.2 含基巖隆起的水平圓盤(pán)響應(yīng)

對(duì)于含有三維起伏界面的海洋地電模型可控源電磁響應(yīng)的計(jì)算目前仍然是一項(xiàng)非常困難的工作,往往需要使用復(fù)雜的有限元技術(shù)才能加以解決.下面利用NMM法研究考察具有軸對(duì)稱(chēng)起伏地層邊界情況下的海洋可控源電磁響應(yīng).

假定圖1(a)中基巖電導(dǎo)率為σbs=0.05 S/m、基巖球冠隆起的高度hcrn=600 m、相應(yīng)的球體半徑acrn=3000 m,而圓盤(pán)參數(shù)、沉積層等參數(shù)與圖3相同,從而形成一個(gè)含有圓盤(pán)與基巖隆起的海洋地電模型.數(shù)值模擬過(guò)程中,將基巖隆起劃分成N=25個(gè)等厚度薄水平圓盤(pán),圖4是兩個(gè)不同工作頻率下計(jì)算得到的Ex和Hy的MVO與PVO曲線.對(duì)比圖2、圖3以及圖4中的數(shù)值結(jié)果可以看到一個(gè)十分有趣的現(xiàn)象：Ex和Hy的MVO曲線均隨源距增加而快速衰減、頻率越高振幅衰減越快;由于基巖隆起的影響,在遠(yuǎn)場(chǎng)處,圖4中Ex和Hy的振幅比圖2和圖3中的振幅更小.此外,不同地電模型上Ex和Hy的PVO曲線變化特征也相差較大.

圖4 包含水平圓盤(pán)與基巖隆起模型的海洋地電模型上由NMM法計(jì)算得到的主測(cè)線上電磁場(chǎng)水平分量(Ex和Hy)振幅與相位曲線 (a)Ex振幅曲線;(b)Ex相位曲線;(c)Hy振幅曲線;(d)Hy相位曲線Fig.4.The in line M VO and PVO of both Ex and Hy ob tained by NMM in the horizontal layered m arine geoelectric model including both resistive disk and spherical cap type bed rock up lift:(a)MVO of Ex; (b)PVO of Ex;(c)M VO of Hy;(d)PVO of Hy.

為進(jìn)一步了解可控源海洋電磁場(chǎng)的空間分布特征,圖5(a)和圖5(b)分別是xOz鉛垂面上圓盤(pán)與基巖隆起周?chē)碾妶?chǎng)與電流密度實(shí)分量的向量圖以及振幅灰度圖.結(jié)果顯示,電場(chǎng)與電流密度均隨著源距增加而快速衰減,但在高阻圓盤(pán)內(nèi)部與其邊界周?chē)约案咦杌鶐r隆起邊界附近,由于積累面電荷影響,電場(chǎng)強(qiáng)度的振幅明顯增大.但電流密度振幅與電場(chǎng)強(qiáng)度變化特征正好相反,由于積累面電荷對(duì)電流的排斥,導(dǎo)致高阻層中電流密度明顯小于低阻層中的電流密度.此外,在高阻圓盤(pán)和高阻基巖隆起的邊界附近,電場(chǎng)強(qiáng)度與電流密度的方向也產(chǎn)生了非常明顯的變化.

圖5 (網(wǎng)刊彩色)基巖隆起與水平圓盤(pán)周?chē)妶?chǎng)強(qiáng)度與電流密度的實(shí)分量向量圖與振幅灰度圖 (a)電場(chǎng)(f=0.25 Hz);(b)電流密度(f=0.25 Hz)Fig.5.(color on line)The vector diagram and am p litude grayscale of real com ponents of electrical intensity(f=0.25 Hz)(a)and electrical current intensity (f=0.25 Hz)(b)in the horizontal layered m arine geoelectricmodel including both resistive d isk and spherical cap type bed rock up lift.

3.3 含基巖凹陷的高阻球體響應(yīng)

圖6是由球形散射體和球冠型基巖凹陷形成的海洋地電模型與其等價(jià)模型示意圖,高阻球形散射體的半徑和電導(dǎo)率分別為asph=700 m和σ=0.01 S/m,其中心到海床面距離dsph=750m,基巖凹陷對(duì)應(yīng)的球冠高度和球體半徑分別是hcrn=500 m和acrn=3000 m.模型中沉積層與基巖電導(dǎo)率與圖1中模型完全相同,其值分別為σsed=0.667 S/m和σbs=0.05 S/m.在數(shù)值模擬中,將整個(gè)球形散射體和球冠基巖凹陷所在區(qū)域分別化分成N1=42和N2=25個(gè)等厚度薄層,并根據(jù)其邊界位置用不同半徑的薄水平圓盤(pán)加以逼近.圖7是兩個(gè)不同工作頻率下Ex和Hy的MVO與PVO曲線.對(duì)比圖7與圖4的計(jì)算結(jié)果不難看出,兩者差異很大.由于球形散射體的體積比圓盤(pán)大得多在加之存在基巖凹陷,導(dǎo)致Ex和Hy振幅隨源距增加衰減更快,工作頻率變化對(duì)振幅與相位的影響也更明顯.

圖6 由球形散射體和球冠基巖凹陷形成的海洋地電模型(a)與其等價(jià)模型示意圖(b)Fig.6.M arine geo-electric model with a resistive spherical scatter and spherical cap type bed rock dep ression(a) and its equivalent model consisting of horizontal layered inhom ogeneous beds(b).

圖7 球形散射體和球冠基巖凹陷模型上由NMM法得到的主測(cè)線上Ex,Hy的MVO與PVO曲線 (a)Ex振幅曲線; (b)Ex相位曲線;(c)Hy振幅曲線;(d)Hy相位曲線Fig.7.The in line MVO and PVO of both Ex and Hy ob tained by NMM in the horizontal layered m arine geoelectric model including both a resistive spherical scatter and spherical cap type bed rock dep ression:(a)MVO of Ex; (b)PVO of Ex;(c)MVO of Hy;(d)PVO of Hy.

圖8 (網(wǎng)刊彩色)球形散射體和球冠基巖凹陷周?chē)妶?chǎng)強(qiáng)度與電流密度的實(shí)分量向量圖與振幅灰度圖 (a)電場(chǎng)(f=0.25 Hz); (b)電流密度(f=0.25 Hz)Fig.8. (color on line)The vector diagram and am p litude grayscale of real com ponents of electrical intensity(f = 0.25 Hz)(a)and electrical cu rrent intensity(f=0.25 Hz) (b)in the horizontal layered m arine geoelectric model including both a resistive spherical scatter and spherical cap type bed rock dep ression.

圖8(a)和圖8(b)分別是xOz鉛垂面上球形散射體和基巖凹陷周?chē)碾妶?chǎng)與電流密度實(shí)分量的向量圖以及振幅灰度圖.同樣可以看到電阻率邊界上積累面電荷對(duì)電場(chǎng)和電流密度空間分布的影響.在整個(gè)高阻球形散射體內(nèi)部及其邊界周?chē)?電場(chǎng)強(qiáng)度明顯增大,但在高阻球體內(nèi)部電流密度明顯減小,而在其邊界周?chē)浇娏髅芏让黠@增大.此外在基巖凹陷界面附近,也可以看到積累面對(duì)周?chē)妶?chǎng)的影響,但由于受高阻球形散射體的屏蔽作用以及離發(fā)射源距離更遠(yuǎn)因素等影響,與球體散射體表面附近相比,其對(duì)周?chē)妶?chǎng)的影響程度要小得多.

3.4 含海底隆起、基巖凹陷的高阻球體響應(yīng)

圖9 由球狀散射體、球冠型海底隆起與基巖球冠凹陷組成的海洋地電模型(a)及其等價(jià)的水平層狀非均質(zhì)模型(b)Fig.9.M arine geo-electricmodelwith a resistive sphericalscatter,topography with spherical cap up lift and spherical cap bed rock dep ression(a)and its equivalent model consisting of horizontal layered inhom ogeneous beds(b).

圖10 球狀散射體、球冠型海底隆起與基巖凹陷組成的海洋地電模型上由NMM法得到的主測(cè)線上水平電磁場(chǎng)(Ex和Hy)的MVO與PVO曲線 (a)Ex振幅曲線;(b)Ex相位曲線;(c)Hy振幅曲線;(d)Hy相位曲線.Fig.10.The inline MVO and PVO of both Ex and Hy obtained by NMM in them arine geoelectric model including both a resistive spherical scatter,topography with spherical cap up lift and spherical cap bed rock dep ression: (a)MVO of Ex;(b)PVO of Ex;(c)M VO of Hy;(d)PVO of Hy.

為進(jìn)一步考察海底地形中海洋電磁的響應(yīng)特征,在圖6地電模型的基礎(chǔ)上,在海底面上增加一個(gè)高度hcrn1=300 m的球冠隆起,形成一個(gè)同時(shí)包含起伏海床面、高阻球形散射體、基巖凹陷等更加復(fù)雜的海洋地電模型,并且接收器按等水平距離放置在海床面上(見(jiàn)圖9),而球冠對(duì)應(yīng)的球體半徑仍然為acrn=3000 m.此外,圖9中的高阻球形散射體、基巖凹陷、沉積層和基巖電導(dǎo)率等參數(shù)與圖6中的地電模型完全相同.在進(jìn)行數(shù)值模擬時(shí),將沉積層隆起、球形散射體和基巖凹陷分別劃分成N1=25,N2=42和N3=25個(gè)等厚度薄層,形成一系列不同半徑的薄水平圓盤(pán).圖10是兩個(gè)不同工作頻率下Ex和Hy的MVO與PVO曲線,不難看出由于起伏地形的影響,Ex和Hy的振幅與相位曲線與圖7的結(jié)果差異十分明顯,沉積層隆起導(dǎo)致Ex和Hy的MVO與PVO曲線變化范圍明顯減小.

圖11(a)和圖11(b)分別是xOz鉛垂面上海底隆起、球形散射體和基巖凹陷周?chē)碾妶?chǎng)與電流密度實(shí)分量的向量圖以及振幅灰度圖.從圖11同樣可以看到海底隆起、球形散射體和球冠基巖凹陷邊界上積累面電荷對(duì)整個(gè)電場(chǎng)強(qiáng)度與電流密度空間分布的影響.

圖11 (網(wǎng)刊彩色)球狀散射體、球冠型海底隆起與基巖凹陷組成的海洋地電模型上電場(chǎng)強(qiáng)度(a)與電流密度(b)的實(shí)分量向量圖與振幅灰度圖 (a)電場(chǎng)(f=0.25 Hz);(b)電流密度(f=0.25 Hz)Fig.11.(color on line)The vector diagram and am p litude grayscale of real com ponents of electrical intensity (f=0.25 Hz)(a)and electrical current intensity(f=0.25 Hz)(b)in them arine geoelectric model including a resistive spherical scatter,topography with spherical cap up lift and spherical cap bed rock dep ression.

4 結(jié) 論

針對(duì)海洋環(huán)境中廣泛存在水平圓盤(pán)、球體、球冠等具有軸對(duì)稱(chēng)電導(dǎo)率分布的散射體,可用多個(gè)不同半徑與厚度的水平薄圓盤(pán)加以逼近,將復(fù)雜的三維海洋地電模型轉(zhuǎn)化為具有軸對(duì)稱(chēng)的水平層狀非均質(zhì)地電模型.利用電阻率空間分布的軸對(duì)稱(chēng)特征,可以將位于對(duì)稱(chēng)軸上的水平電偶極子天線電磁場(chǎng)正演模擬問(wèn)題轉(zhuǎn)化為兩個(gè)軸對(duì)稱(chēng)問(wèn)題加以求解,并通過(guò)數(shù)值模式匹配法可以得到海洋可控源電磁場(chǎng)的半解析解,實(shí)現(xiàn)海洋可控源三維電磁響應(yīng)以及全空間電磁場(chǎng)的快速計(jì)算.

數(shù)值計(jì)算結(jié)果顯示,對(duì)于水平層狀地層模型,數(shù)值模式匹配法的計(jì)算精度可以得到傳輸線的水平;對(duì)于含有單一水平圓盤(pán)散射體模型,數(shù)值模式匹配法的數(shù)值模擬結(jié)果與3D FV算法符合得也很好,但數(shù)值模式匹配法算法效率和計(jì)算精度比3D FV算法更高,且占用的內(nèi)存也更少.此外,圓盤(pán)、球體、海底隆起以及基巖隆起和凹陷等對(duì)可控源電磁勘探中的電場(chǎng)與電流密度空間分布均有非常明顯影響,不同地電模型上Ex和Hy的MVO與PVO曲線往往差異巨大.在高阻散射體(圓盤(pán)、球體)內(nèi)部與其邊界周?chē)约斑吔缏∑鸹虬枷萁绺浇?由于積累面電荷影響,電場(chǎng)明顯增強(qiáng);但電流密度振幅與電場(chǎng)強(qiáng)度變化特征正好相反,積累面電荷對(duì)電流的排斥作用,導(dǎo)致高阻層內(nèi)部電流密度明顯減小.此外,在高阻散射體(圓盤(pán)、球體)與起伏地層邊界上,電場(chǎng)強(qiáng)度與電流密度的方向也產(chǎn)生了非常明顯的變化.

附錄A無(wú)限厚層柱狀非均質(zhì)地層中電磁場(chǎng)的基本解

數(shù)值模式匹配算法的關(guān)鍵是需要事先確定每個(gè)地層(無(wú)限大均勻或柱狀地層)中電磁場(chǎng)的基本解.為此,用有限長(zhǎng)區(qū)間[ρMN,ρMX]逼近半無(wú)限區(qū)間,并將其剖分成M(偶數(shù))個(gè)小區(qū)間,剖分節(jié)點(diǎn)位置用ρα(α=1,2,···,M+1)表示.選用二次New ton插值函數(shù)作為徑向節(jié)點(diǎn)ρα上的基函數(shù)φα(ρ)(α=1,···,M)[16],根據(jù)方程(10)的第二類(lèi)齊次邊界條件以及方程(11)的截?cái)鄺l件,選擇M個(gè)基函數(shù)組成列向量S(ρ)=(φ1(ρ),φ2(ρ),···,φM-1(ρ),φM(ρ))T,利用Petrov-Galerkin法將方程(6)中的±1階水平磁場(chǎng)分量展開(kāi)為如下形式[]:

求解(A 2)式可以得到2M個(gè)本征值κ2k,α和相應(yīng)的本征向量Hk,α(α=1,2,···,2M).將本征值代入(A 3)式中并進(jìn)行求解得:

進(jìn) 一 步 求 解 方 程(A 4)可 確 定 向 量Ck=,其中,是本征向量矩陣.最后,我們得到無(wú)限大均勻或柱狀地層中k階諧變水平磁場(chǎng)分量的基本解[16]:

根據(jù)每個(gè)地層上電導(dǎo)率大小和空間分布,利用上面的方法可以計(jì)算出各個(gè)地層中本征值和本征向量矩陣和,從而確定每個(gè)地層的基本解.

此外,利用方程(12)和(14)可以得到水平電場(chǎng)分量的基本解

以及垂直電磁場(chǎng)分量的基本解

[1]Edwards N 2005 Surv.Geophys.26 675

[2]Constab le S 2010 Geophysics 75 75A 67

[3]Yuan J,Edwards N 2000 Geophys.Res.Lett.27 2397

[4]W eiss C J,Constab le S 2006 Geophysics 71 G 321

[5]Constab le S C,W eiss C J 2006 Geophysics 71 G 43

[6]Hoversten G M,Newm an G A,Geier A,Flanagan G 2006 Geophysics 71 G 239

[7]W ang J X,W ang H N,Zhou J M,Yang SW,Liu X J, Y in C C 2013 Acta Phys.Sin.62 224101(in Chinese) [汪建勛,汪宏年,周建美,楊守文,劉曉軍,殷長(zhǎng)春2013物理學(xué)報(bào)62 224101]

[8]Li Y G,Dai S K 2011 Geophys.J.In t.185 622

[9]Xu Z F,Wu X P 2010 Chinese J.Geophys.53 1931(in Chinese)[徐志鋒,吳小平2010地球物理學(xué)報(bào)53 1931]

[10]Shen J S 2003 Chin.J.Geophys.46 280(in Chinese) [沈金松2003地球物理學(xué)報(bào)46 280]

[11]Streich R 2009 Geophysics 74 F95

[12]Zhou J M,Zhang Y,W ang H N,Yang SW,Y in C C 2014 Acta Phys.Sin.63 159101(in Chinese)[周建美,張燁,汪宏年,楊守文,殷長(zhǎng)春2014物理學(xué)報(bào)63 159101]

[13]Chen G B,W ang H N,Yao J J,Han Z Y,Yang S W 2009 Acta Phys.Sin.58 1608(in Chinese)[陳桂波,汪宏年,姚敬金,韓子夜,楊守文2009物理學(xué)報(bào)58 1608]

[14]Chen G B,Wang H N,Yao J J,Han Z Y 2009 Acta Phys.Sin.58 3848(in Chinese)[陳桂波,汪宏年,姚敬金,韓子夜2009物理學(xué)報(bào)58 3848]

[15]Kong F N,Johnstad S E,R?sten T,W esterdah l H 2008 Geophysics 73 F9

[16]W ang H N,Tao H G,Yao J J,Zhang Y 2012 IEEE Trans.Geosci.Rem ote Sens.50 3383

[17]W ang H N,Tao H G,Yang SW 2008 Chin.J.Geophys. 51 1591

[18]Liu Q H,Chew W C 1992 Radio Sci.27 569

[19]W ang H N 2011 IEEE Trans.Geosci.Rem ote Sens.49 4483

[20]W ang H N,So P M,Yang SW,Hoefer W J R,Du H L 2008 IEEE Trans.Geosci.Rem ote Sens.46 1134

[21]Zhu T Z,Yang SW,Bai Y,Chen T,W ang H N 2017 Chin.J.Geophys.60 1221(in Chinese)[朱天竹,楊守文,白彥,陳濤,汪宏年2017地球物理學(xué)報(bào)60 1221]

[22]Chew W C 1990 W aves and Fields in Inhom ogeneous M edia(New York:van Nostrand Reinhold)

(Received 17 Novem ber 2016;revised manuscript received 10 April 2017)

Efficient simulation of marine controlled source electromagnetic responses for axisymmetric scatter by using numerical mode matching approach?

Lin Lin Jiao Li-Guang Chen Bo Kang Zhuang-Zhuang Ma Yu-Gang Wang Hong-Nian?

(College of Physics,Jilin University,Changchun 130012,China)

Horizontal disk,sphere,and spherical crown are a very im portant type of scatter in geophysics research.In the marine environm ent,a disk-like scatter can be used to describe several resistive targets,e.g.,basaltic sills and stratigraphic hyd rocarbon reservoirs while spherical crown can be used to approximately depict the topography of interface for basem ent rock.This type of scatter has characteristics of axisymm etrical distribution of the conductivity. If som e app roaches can be established to efficiently simulate the m arine controlled source electrom agnetic(MCSEM) response to this scatter,it w ill be meaningful to investigate the nature of MCSEM responses in com p lex formation and to build app ropriate method of processing and explaining MCSEM data.In this paper,the resistive scatters are approximated by one or several horizontal concentric diskswith different radiiand thickness values,based on the axially symm etrical spatial distribution of conductivity.Then,a combination of these concentric disks with air,sea water and surrounding bedsw ill construct a horizontally stratified inhom ogeneous form ation with comm on axis-center,whose spatial distribution of conductivity is layered in the verticaldirection and axisymmetric in the horizontal direction.Based on the approxim ationsm entioned above,the com putation of MCSEM response excited by horizontal electrical dipole (HED)located at the z-axis is entirely transform ed into two axially symmetrical prob lem s for the Fourier harmonic com ponents of the electromagnetic(EM)fields.The differential operators about the horizontalmagnetic com ponents and transform ation of horizontal electrical com ponents and other EM com ponents from horizontalm agnetic com ponents are derived.Then,the num ericalm ode m atching approach is extended to the simulation of the EM field and threedimensional(3D)MCSEM responses excited by the HED in the formation.The p rocedure for solving the EM field is presented.The sem i-analytic solution of EM field in the whole space is obtained to efficiently and num erically model MCSEM response in the com p lex formation.Finally,the efficiency and accuracy of the presentmethod are demonstrated num erically.The characteristics of 3D MCSEM responses in three different cases are further investigated.

marine controlled-source electromagnetic method,axisymmetric scatter,numerical mode matching approach,sem i-analytic solution

PACS:91.25.Qi,02.30.Zz,41.20.-q DO I:10.7498/aps.66.139102

?國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào):41574110)和國(guó)家高技術(shù)研究發(fā)展計(jì)劃重大項(xiàng)目(批準(zhǔn)號(hào):2012AA 09A 20103)資助的課題.

?通信作者.E-m ail:wanghn@jlu.edu.cn

PACS:91.25.Qi,02.30.Zz,41.20.-q DO I:10.7498/aps.66.139102

*Project supported by the National Natural Science Foundation of China(G rant No.41574110)and the National High-tech R&D Program,M ajor Pro ject,China(G rant No.2012AA 09A 20103).

?Corresponding author.E-m ail:wanghn@jlu.edu.cn