謝 慧 萬 露 岳智彬

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

1 引言

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

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

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

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

2.1 基本理論

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

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

對式（1）兩邊同時進行旋度運算，將式（3）代入后得到電場強度矢量的雙旋度方程表達式：

將電場分解為背景場（一次場）和散射場（二次場）：

其中EP為一次場電場強度矢量，ES為二次場電場強度矢量。

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

將式（5）與式（6）相減得到二次場電場強度矢量滿足的方程表達式：

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

其中Γ 為Ω 的邊界。

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

2.2 網(wǎng)格剖分

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

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

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

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

2.3 矢量插值函數(shù)

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

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

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

將上式相加得六面體子單元中總的電場強度矢量：

2.4 單元分析和總體合成

將式（13）按照一次場和二次場拆分可得：

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

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

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

經(jīng)化簡可得：

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

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

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

3 實驗結(jié)果與分析

3.1 準確性驗證

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

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

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

圖5 本文算法計算3D模型的準確性驗證

3.2 異常體模型算例

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

圖6 三維大地模型

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

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

4 結(jié)語

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