謝 慧 萬(wàn) 露 岳智彬

（海軍工程大學(xué)電子工程學(xué)院 武漢 430033）

1 引言

水平低架天線以其架設(shè)簡(jiǎn)單、易于維修等優(yōu)點(diǎn)廣泛應(yīng)用于長(zhǎng)波超長(zhǎng)波極長(zhǎng)波無(wú)線電工程中［1～5］，著名的Beverage 天線就是電流為行波分布的水平低架天線［6］。水平低架天線的輻射強(qiáng)度與天線布設(shè)場(chǎng)地的電阻率的開方成正比，天線架設(shè)場(chǎng)地電阻率越大，電磁波輻射強(qiáng)度越大，天線輻射效率越高［7］。根據(jù)前人的研究表明天線架設(shè)場(chǎng)地電阻率對(duì)天線的發(fā)射效率有著重大的影響，但當(dāng)前實(shí)際針對(duì)各向異性大地，以及復(fù)雜多維架設(shè)大地對(duì)水平天線電磁場(chǎng)的影響研究仍然很少。

目前我國(guó)僅在大地電磁勘探中使用可控源電磁法正演方法對(duì)各向異性大地以及復(fù)雜多維地形中發(fā)射源輻射出的電磁場(chǎng)進(jìn)行模擬計(jì)算?？煽卦匆纛l大地電磁領(lǐng)域的研究大多以一維二維的各向異性介質(zhì)為主［8～12］，其中Everett 首先展開了對(duì)一維橫向各向異性海洋可控源電磁的響應(yīng)研究［8］，L?seth 等提出了一種適用于各向異性層狀介質(zhì)中偶極子源電磁場(chǎng)的計(jì)算方法［9］，羅鳴等在Loseth 的基礎(chǔ)上，利用正余弦變換實(shí)現(xiàn)了二維傅里葉變換，由此將波數(shù)域電磁場(chǎng)轉(zhuǎn)換至空間域，實(shí)現(xiàn)了海洋可控源電磁在任意各向異性層狀介質(zhì)中響應(yīng)的計(jì)算［10］，Wiik等采用積分方程法實(shí)現(xiàn)了各向異性介質(zhì)中的2.5 維正演［11］。直到近些年，電導(dǎo)率各向異性三維正演開始受到國(guó)內(nèi)外學(xué)者的重視［13～16］。Wang 等推導(dǎo)并實(shí)現(xiàn)了一種可控音頻大地電磁電阻率軸向各向異性的三維（3D）交錯(cuò)網(wǎng)格有限差分?jǐn)?shù)值模擬方法［13］，邱長(zhǎng)凱等采用非結(jié)構(gòu)四面體網(wǎng)格進(jìn)行剖分的矢量有限元法實(shí)現(xiàn)了可控源音頻大地電磁三維非結(jié)構(gòu)正演［14］，He等基于分解的二次場(chǎng)，采用矢量有限元法實(shí)現(xiàn)了三維任意各向異性正演［15］。

水平低架天線所采用的頻率范圍為3Hz～30kHz，大地電阻率范圍為1Ω?m～5000Ω?m，對(duì)位移電流與傳導(dǎo)電流的比值進(jìn)行估算，其中ω為角頻率，ρ為大地電阻率，f為頻率，位移電流與傳導(dǎo)電流的最大比值7.5×10-3較大，因此不能忽略位移電流。由于可控源電磁法采用的發(fā)射頻率非常低，可忽略位移電流。因此，本文根據(jù)可控源電磁法的正演算法提出了一種新的方法，將位移電流計(jì)算在內(nèi)，使用基于規(guī)則四面體網(wǎng)格的矢量有限元法，實(shí)現(xiàn)了架設(shè)場(chǎng)地各向異性，且呈現(xiàn)多維復(fù)雜結(jié)構(gòu)時(shí)，空間任意點(diǎn)的電磁響應(yīng)數(shù)值模擬。

2 三維各向異性矢量有限元法

2.1 基本理論

從麥克斯韋方程組的微分形式出發(fā)，通過(guò)傅里葉變換將時(shí)變電磁場(chǎng)轉(zhuǎn)變?yōu)橹C變電磁場(chǎng)，選取e-iωt作為時(shí)諧因子，則麥克斯韋方程組可寫為

其中E是電場(chǎng)強(qiáng)度矢量，H是磁場(chǎng)強(qiáng)度矢量，Js是外電流密度矢量，J是傳導(dǎo)電流密度矢量，JD是位移電流密度矢量，、ω、ε分別表示電導(dǎo)率張量，角頻率和介電常數(shù)。

對(duì)式（1）兩邊同時(shí)進(jìn)行旋度運(yùn)算，將式（3）代入后得到電場(chǎng)強(qiáng)度矢量的雙旋度方程表達(dá)式：

將電場(chǎng)分解為背景場(chǎng)（一次場(chǎng)）和散射場(chǎng)（二次場(chǎng)）：

其中EP為一次場(chǎng)電場(chǎng)強(qiáng)度矢量，ES為二次場(chǎng)電場(chǎng)強(qiáng)度矢量。

由總電場(chǎng)強(qiáng)度矢量的雙旋度方程式可以推導(dǎo)出一次場(chǎng)電場(chǎng)強(qiáng)度所滿足的雙旋度方程表達(dá)式：

將式（5）與式（6）相減得到二次場(chǎng)電場(chǎng)強(qiáng)度矢量滿足的方程表達(dá)式：

結(jié)合狄利克雷邊界條件（第一類邊界條件），對(duì)應(yīng)的邊值問(wèn)題表達(dá)式如下：

其中Γ 為Ω 的邊界。

采用加權(quán)余量法中的伽里金法得到式（8）對(duì)應(yīng)的加權(quán)余量方程：

2.2 網(wǎng)格剖分

圖1 為所構(gòu)建的三維張量的可控源音頻大地電磁法模型圖，該模型的研究區(qū)域由空氣和大地兩大部分組成。如圖所示，在地面的某位置放置了沿X 方向有限長(zhǎng)導(dǎo)線源，地面下的地底層中存在一個(gè)三維異常結(jié)構(gòu)體。

圖1 三維張量的可控源音頻大地電磁法模型圖

圖2 是對(duì)三維張量的CSAMT 模型進(jìn)行網(wǎng)格剖分的示意圖，采用采用連續(xù)、規(guī)則的四面體網(wǎng)格剖分單元對(duì)模型剖分，便于構(gòu)建結(jié)構(gòu)相對(duì)簡(jiǎn)單的插值函數(shù)。

圖2 三維張量的可控源音頻大地電磁法模型網(wǎng)格剖分示意圖

2.3 矢量插值函數(shù)

圖3是從圖2整體模型中取出的一個(gè)六面體矢量子單元。通過(guò)矢量有限元法把單元內(nèi)的各個(gè)方向的電場(chǎng)強(qiáng)度矢量（磁場(chǎng)強(qiáng)度矢量）等效到與坐標(biāo)軸同向的棱邊上，電場(chǎng)強(qiáng)度矢量（磁場(chǎng)強(qiáng)度矢量）的正方向也就是坐標(biāo)軸或者棱邊的正方向。

圖3 矢量有限單元的六面體子單元

采用Whitney 矢量基函數(shù)線性插值［17］，假設(shè)規(guī)則六面體子單元沿則X，Y，Z 軸方向的棱邊邊長(zhǎng)分別為lX、lY、lZ，中心點(diǎn)坐標(biāo)為(x0，y0，z0)則該六面體子單元內(nèi)坐標(biāo)軸方向上總的電場(chǎng)強(qiáng)度為

將上式相加得六面體子單元中總的電場(chǎng)強(qiáng)度矢量：

2.4 單元分析和總體合成

將式（13）按照一次場(chǎng)和二次場(chǎng)拆分可得：

對(duì)網(wǎng)格剖分后的每個(gè)六面體子單元內(nèi)積分可得：

其中Ωi代表當(dāng)前六面體子單元的區(qū)域方位。

然后在整個(gè)區(qū)域Ω 上對(duì)所有的六面體子單元進(jìn)行求和得：

經(jīng)化簡(jiǎn)可得：

其中ne為模型整體區(qū)域Ω 中所包含六面體子單元數(shù)量，K1e、K2e、K3e是12*12的矩陣［17］。

其中K1、K2、K3是經(jīng)過(guò)求和后所有六面體子單元所包含棱邊信息的矩陣，K1-K2就是總體系數(shù)矩陣。

根據(jù)計(jì)算出來(lái)的電場(chǎng)強(qiáng)度矢量值，結(jié)合式（1）可以推導(dǎo)出磁場(chǎng)強(qiáng)度矢量的表達(dá)式：

3 實(shí)驗(yàn)結(jié)果與分析

3.1 準(zhǔn)確性驗(yàn)證

如圖4 所示，建立水平層狀電阻率各向異性模型。大地第一層為各向同性地層，電阻率為1000 Ω ?m，層厚度為100m；第二層為各向異性地層，x、y和z 方向電阻率分別為2000 Ω ?m、2000 Ω ?m 和500 Ω ?m。第三層為各向同性地層，電阻率為3000 Ω ?m。沿y 軸放置長(zhǎng)度為1m 的發(fā)射源，中心點(diǎn)位于空氣層（0，0，-1），發(fā)射頻率為3kHz。

圖4 水平層狀各向異性模型

如圖5 所示，將計(jì)算電場(chǎng)強(qiáng)度值結(jié)果與Yin（2003）的層狀各向異性解析解比較，本文計(jì)算結(jié)果（Ex和Ez，Ey=0）與Yin（2003）解析解最大誤差均在4%以內(nèi)，驗(yàn)證了本文提出的算法的準(zhǔn)確性。

圖5 本文算法計(jì)算3D模型的準(zhǔn)確性驗(yàn)證

3.2 異常體模型算例

為了探求天下架設(shè)場(chǎng)地存在異常體時(shí)，對(duì)水平線天線的電場(chǎng)強(qiáng)度的影響，設(shè)計(jì)如圖6 所示的三維大地模型。沿x 軸放置長(zhǎng)度為100m，發(fā)射電流為1A，架設(shè)高度為10m 的水平電偶極子，其上表面中心點(diǎn)為（0，0，10），發(fā)射頻率為3kHz。假設(shè)背景電阻率為1000Ω?m，異常體是位于地表的河流、藻地、裸巖，故設(shè)置異常體大小為10m×50m×30m 的各向同性均勻介質(zhì)，電阻率為5Ω?m（河流）、500Ω?m（藻地）、2500Ω?m（低阻裸巖）、4500Ω?m（高阻裸巖）。接收點(diǎn)為（50000，0，-1）

圖6 三維大地模型

將異常體沿天線徑向移動(dòng)，曲線1000Ω?m表示無(wú)異常體，天線架設(shè)場(chǎng)地均勻，電阻率為1000Ω?m，其電場(chǎng)響應(yīng)為背景電場(chǎng)Es。從圖7 異常體位置與電場(chǎng)強(qiáng)度的變化關(guān)系圖可以看出，當(dāng)天線架設(shè)場(chǎng)地含有異常體，且異常體電阻率小于大地背景電阻率1000Ω?m 時(shí)，x方向電場(chǎng)強(qiáng)度Ex均小于背景電場(chǎng)Es，且異常體電阻率越小，所激發(fā)的電場(chǎng)值越小；當(dāng)異常體位于原點(diǎn)，即正對(duì)天線中心點(diǎn)時(shí)，電場(chǎng)強(qiáng)度最大。當(dāng)異常體電阻率大于大地背景電阻率1000Ω?m 時(shí)，x方向電場(chǎng)強(qiáng)度Ex均大于背景電場(chǎng)Es，且異常體電阻率越大，所激發(fā)的電場(chǎng)值越大；當(dāng)異常體位于原點(diǎn)，即正對(duì)天線中心點(diǎn)時(shí)，電場(chǎng)強(qiáng)度最小。沿y方向電場(chǎng)強(qiáng)度Ey相較于Ex和z方向電場(chǎng)強(qiáng)度Ez，故在此不做分析。Ez總電場(chǎng)強(qiáng)度E 變化規(guī)律與Ex相同。

圖7 異常體位置與電場(chǎng)強(qiáng)度的變化關(guān)系

4 結(jié)語(yǔ)

本文提出的新方法，通過(guò)復(fù)電導(dǎo)率將傳導(dǎo)電流和位移電流合并考慮，成功將位移電流計(jì)算在內(nèi)。從以上研究分析得出，水平線天線輻射電場(chǎng)強(qiáng)度受異常體影響較大。異常體電導(dǎo)率與背景電導(dǎo)率差值越大，影響越大；異常體電阻率大于背景電阻率，越接近天線兩端，輻射能力越強(qiáng)；異常體電阻率小于背景電阻率，接近天線中心，可有效減少低阻異常體降低天線輻射能力的影響。