時(shí)宗洋，趙一宇，渠曉東，馬力超

（1.北京機(jī)械設(shè)備研究所，北京 100854；2.中國科學(xué)院空天信息創(chuàng)新研究院，北京 100094）

0 引言

典型層狀海洋模型下任意姿態(tài)電性天線電磁場的模擬計(jì)算方法主要應(yīng)用海洋目標(biāo)電磁探測應(yīng)用領(lǐng)域，如海底石油、水合物等探測的海洋可控源電磁方法和海水中大型目標(biāo)如沉船等探測等。

以海洋可控源電磁方法（MCSEM）為例，其通常采用幾百米長的水平電性天線在海水中（距海底幾十米的位置）輻射峰值電流幾百安培至千安培、基頻在n×10-1Hz至n×10 Hz范圍內(nèi)的矩形波電流，通過布置在海底或者拖曳在距水平電性天線固定偏移距的電場或磁場傳感器觀測電場/磁場響應(yīng)信號(hào)[1-2]，然后通過采用相適應(yīng)的數(shù)據(jù)處理手段對(duì)電磁信號(hào)進(jìn)行處理，運(yùn)用預(yù)先建立的層狀海洋模型正反演算法獲得對(duì)實(shí)測電磁信號(hào)的定量反演解釋，從而得到海洋中目標(biāo)電阻率信息，圖1給出了海洋可控源電磁法工作示意圖，圖中傳感器陣列rxi(i＝1…N)位于海底沿電性天線軸向布置。

圖1 海洋可控源電磁法(MCSEM)工作示意圖

現(xiàn)有層狀海洋模型電性天線輻射電磁場的模擬計(jì)算方法主要有2種，一種是將有限長電性天線看作電偶極子天線[3-4]，然后利用層狀海洋模型下電偶極子天線輻射電磁場模擬解釋海洋目標(biāo)及海底電性參數(shù)；另一種方法是將有限長度電性天線進(jìn)行均勻密集分割，將分割后每一段等效為電偶極子，再對(duì)每個(gè)等效電偶極子天線的電磁響應(yīng)求和得到有限長電性天線的輻射電磁場。

以上計(jì)算方法存在以下缺點(diǎn)或不足之處，方法一將有限長電性天線看作電偶極子天線，這種方法認(rèn)為：當(dāng)海洋目標(biāo)（海底或海水目標(biāo)）和觀測傳感器距離電性天線很遠(yuǎn)（相比于電性天線長度）時(shí)電性天線可以等效為電偶極子天線。然而，實(shí)際MCSEM中的電性天線長幾百米（一般100～300 m），距離海底幾十米距離，航行作業(yè)時(shí)，電性天線從靠近海底的傳感器陣列到遠(yuǎn)離傳感器陣列的過程中，無法始終保證電性天線距離海底和傳感器的距離遠(yuǎn)大于電性天線的長度（如5倍電性天線長度），從而導(dǎo)致該方法在計(jì)算小偏移距的電磁場時(shí)產(chǎn)生較大的計(jì)算誤差。方法二均勻密集分割有限長度電性天線，然后把分割后的每段等效為電偶極子，并對(duì)每個(gè)等效電偶極子天線的電磁響應(yīng)求和得到有限長電性天線的輻射電磁場。該方法在保證計(jì)算精度的同時(shí)會(huì)極大地降低計(jì)算效率，不利于工程應(yīng)用中數(shù)據(jù)資料的快速解釋[5]。

本文針對(duì)層狀海洋模型電性天線輻射電磁場的模擬計(jì)算現(xiàn)有方法無法兼顧計(jì)算精度和計(jì)算效率的不足，提出了一種典型層狀海洋模型中任意姿態(tài)電性天線輻射電磁場的非均勻稀疏分割計(jì)算方法，采用切比雪夫多項(xiàng)式的零點(diǎn)作為分點(diǎn)，對(duì)原有的均勻分割方案進(jìn)行了以切比雪夫多項(xiàng)式零點(diǎn)為分點(diǎn)的非均勻稀疏分割方案的改進(jìn)，實(shí)現(xiàn)了有限長電性天線的精確快速計(jì)算，兼顧了計(jì)算的精度和效率。

本文首先介紹典型層狀海洋模型下電性天線輻射電磁場的工程應(yīng)用及現(xiàn)階段主要計(jì)算方法和存在的不足，引出本文建立的非均勻稀疏分割積分方法及其優(yōu)勢；然后構(gòu)建了層狀海洋模型電磁天線輻射電磁場模型，并對(duì)電磁場積分計(jì)算式進(jìn)行了理論推導(dǎo)，并且建立基于高斯切比雪夫分點(diǎn)的非均勻稀疏分割積分方法，給出理論計(jì)算式和算法流程；接下來通過模擬仿真算例，對(duì)電偶極子等效計(jì)算方法、密集均勻分割積分方法和本文提出的非均勻稀疏分割積分這3種方法在電磁場計(jì)算精度和效率方面進(jìn)行對(duì)比分析，論證了本文提出方法在兼顧計(jì)算精度和效率方面的優(yōu)勢，能夠?qū)崿F(xiàn)工程應(yīng)用中數(shù)據(jù)解釋的效率和準(zhǔn)確性的提升；最后對(duì)全文進(jìn)行總結(jié)，給出結(jié)論。

1 層狀海洋電性天線輻射電磁場理論

1.1 典型層狀海洋模型

以典型3層海洋模型為例，即空氣層、海水層和海底層。其中，空氣層設(shè)定為半無限大均勻空間，海底層可以為半無限大均勻空間或?qū)訝羁臻g，海水層深度為d1，設(shè)定分解面互相平行，水平方向上無限延伸，空氣、海水和海底（可以為多層海底）的介電參數(shù)分別為σi、εi（i＝0,1,2,…,n,n≥2），其中i＝0表示空氣層，i＝1表示海水層。真空磁導(dǎo)率為μ0，相對(duì)磁導(dǎo)率μr設(shè)置為1。AB表示電性天線首尾，天線長度為L，天線電極矩為P＝IL，天線中點(diǎn)位于原點(diǎn)O正下方，目標(biāo)所在深度為H，距離海底深度為h，各層厚度為d i。笛卡爾坐標(biāo)系正z方向垂直水平面向下，o-xyz滿足右手螺旋定則。假設(shè)任意姿態(tài)電性天線在XOY平面內(nèi)的投影與x軸夾角為θ，天線與XOY平面的夾角為?。

層狀海洋模型電性天線模型如圖2所示。層狀海洋模型下電性天線輻射電磁場的推導(dǎo)基于以下假設(shè)[6-8]：

圖2 層狀海洋模型中的電性天線

1）滿足準(zhǔn)靜態(tài)近似條件，頻率小于100 k Hz；

2）準(zhǔn)靜態(tài)近似下忽略海水和海底環(huán)境的位移電流，波數(shù)為k2i＝-iωμ0σi，i＝1,2,…,n?？諝庵袃H存在位移電流，約定電導(dǎo)率為σ0＝iωε0，波數(shù)為k20＝

3）空氣、海水及同層海底媒介時(shí)各向同性的，參數(shù)與時(shí)間、溫度和壓強(qiáng)無關(guān)；

4）準(zhǔn)靜態(tài)近似下，認(rèn)為媒介參數(shù)與頻率無關(guān)，磁導(dǎo)率和介電常數(shù)采用真空中的參數(shù)。

對(duì)于沒有自由電荷的空間，電性天線在空間中產(chǎn)生的電磁場滿足如下麥克斯韋方程組。

定義電矢量為A，滿足Η＝?×Α。得到電場與電矢位關(guān)系：

式中，U為標(biāo)量位，滿足

1.2 電偶極子天線輻射電磁場

1.2.1 水平電偶極子天線頻域電磁場

層狀海洋模型下水平電偶極子天線如圖3所示，其對(duì)應(yīng)的電磁場具有對(duì)稱性，此時(shí)的電矢量位A僅包含了2個(gè)分量，即沿偶極矩方向的分量A x和沿垂直海水-空氣分解面的分量A z。電矢量位A在笛卡爾坐標(biāo)系下可表示為：

矢量位、電場及磁場滿足的邊界條件為[8-10]：

1）整個(gè)空間中，除發(fā)射源位置外，矢量位A處處為有限值，且在無窮遠(yuǎn)處，矢量位為零，即Α(r→∞)→0；

2）各層分解面上，電場和磁場的切向分量連續(xù)。

根據(jù)邊界條件約束，電場E、磁場H與矢量為A和標(biāo)量位U之間的關(guān)系，采用分離變量法可得到各層中矢量位各分量的解。

空氣中電場和磁場各分量頻域計(jì)算式見式(3)和式(4)。

式中，?為觀測點(diǎn)與電偶極源位置連線在XOY平面投影與x軸的夾角，sin?＝y(tǒng)/r，cos?＝x/r。

海水中電場和磁場各分量頻域計(jì)算式見式(5)和式(6)。

式中，Idl表示電偶極源的電矩，m為積分變量，r＝(x2＋y2)1/2為收發(fā)距，mi＝(m2-k2i)1/2，i＝0,1,2,…,N。積分核中的Ji(mr),i＝0,1代表第i階貝塞爾函數(shù)。參數(shù)Ci、B i、D i和E i為待定系數(shù)，通過理論推導(dǎo)計(jì)算結(jié)果如下：

式中，N01、N21、P01和P21計(jì)算式如下：

式中，R*2和Q*2的遞推計(jì)算式為：

1.2.2 垂直電偶極子天線頻域電磁場

層狀海洋模型下垂直電偶極子天線示意圖如圖4所示，其電矢位A僅包含沿偶極矩方向的分量A z[11]。此時(shí)的電矢量位A在笛卡爾坐標(biāo)系下表示為：

圖4 層狀海洋模型下垂直電偶極子天線

根據(jù)邊界條件約束，及電場E、磁場H與矢量為A和標(biāo)量位U之間的關(guān)系，采用分離變量法可得到各層中矢量位各分量的解。

空氣中電場和磁場各分量頻域計(jì)算式見式(10)和式(11)。

海水中電場和磁場各分量頻域計(jì)算式見式(12)和式(13)。

式中，參數(shù)Ci和D i為待定系數(shù)，通過理論推導(dǎo)計(jì)算結(jié)果如下：

式中，P01和P21計(jì)算式如下：

式中，R*2的遞推計(jì)算式如下：

1.2.3 任意姿態(tài)電偶極子天線頻域電磁場

海水中任意姿態(tài)電偶極子天線輻射電磁場計(jì)算時(shí)，將偶極矩分別向x,y和z軸投影，然后分別計(jì)算三分量投影電偶極子Px、Py和Pz的電磁場，最后將各方向電偶極子天線的電磁場求和即可[12-13]。層狀海洋模型下任意姿態(tài)電偶極子天線如圖5所示。

圖5 層狀海洋模型下任意姿態(tài)電偶極子天線

式中，沿y方向的電矩為Py的電偶極子天線產(chǎn)生的電磁場，可以通過計(jì)算沿x方向的電矩為Py的電偶極子天線的電磁場然后通過坐標(biāo)變換得到。

1.2.4 電偶極子天線電磁場計(jì)算方法

層狀海洋模型下電偶極子天線輻射電磁場的計(jì)算方法流程如圖6所示。該算法包含了輻射電磁場的頻域和時(shí)域模擬計(jì)算方法，總體思路是首先根據(jù)電偶極子天線輻射電磁場的頻域計(jì)算式通過快速Hankel變換數(shù)值濾波方法實(shí)現(xiàn)頻域模擬計(jì)算，獲得頻域電磁結(jié)果；然后通過GS變換數(shù)值濾波方法實(shí)現(xiàn)快速的頻-時(shí)變換，得到時(shí)域電磁計(jì)算結(jié)果[14-15]。

圖6 層狀海洋模型電偶極子天線輻射電磁場計(jì)算流程

1.3 有限長電性天線輻射電磁場

1.3.1 任意姿態(tài)電性天線輻射電磁場的積分計(jì)算式

理論上，典型海洋模型下任意姿態(tài)有限長電性天線輻射電磁場的計(jì)算需要通過對(duì)天線長度進(jìn)行積分獲得。

假設(shè)電性天線l位置的電偶極子在觀測位置r處產(chǎn)生的響應(yīng)（電場或磁場）表示為f(p,θ,?,l,r)，其 中f(·)可 以表 示 輻 射 電 磁 場E x、E y、E z、Bx、B y、B z的計(jì)算式，θ和?為已知參數(shù)，p＝Idl，p表示天線l位置處的電偶極子極矩，I為天線中電流強(qiáng)度，dl表示l位置的偶極子單位長度，l和r分別為電偶極子天線和觀測點(diǎn)的相對(duì)于坐標(biāo)原點(diǎn)的位置矢量。因此，有限長電性天線的電磁場計(jì)算式為：

1.3.2 現(xiàn)有電性天線輻射電磁場的積分計(jì)算方法

現(xiàn)有的電性天線輻射電磁場的積分計(jì)算方法主要采用電偶極子天線等效方法和均勻密集分割積分方法。對(duì)于電偶極子天線等效方法，是將有限長電性天線看作電偶極子天線，此時(shí)積分式(18)變?yōu)椋?/p>

式中，l p表示有限長電性天線中心位置的位置矢量。

對(duì)于均勻密集分割積分方法是將積分式(18)均勻離散化為眾多的電偶極子天線響應(yīng)之和的形式。

式中，N為均勻分割的分點(diǎn)數(shù)，一般為保證計(jì)算精度，對(duì)于長度100 m以上的天線，N取值不小100。

2 層狀海洋電性天線輻射電磁場的非均勻稀疏分割計(jì)算方法

為解決現(xiàn)有方案在計(jì)算效率和精度無法兼顧的不足，本節(jié)將基于切比雪夫多項(xiàng)式零點(diǎn)建立適于加速計(jì)算電性天線輻射電磁場的非均勻稀疏分點(diǎn)的高斯-切比雪夫積分方法，以保證了計(jì)算精度的同時(shí)極大的減小了積分節(jié)點(diǎn)數(shù)，提升了計(jì)算效率。

2.1 基于高斯—切比雪夫積分的非均勻稀疏分割計(jì)算方法

在高斯勒讓德數(shù)值積分方法在地面有限長電性天線輻射電磁場加速計(jì)算方面的應(yīng)用發(fā)展基礎(chǔ)上[9,16]，本論文提出基于高斯—切比雪夫積分的非均勻稀疏分割方法，利用具備帶權(quán)正交性的切比雪夫多項(xiàng)式的零點(diǎn)作為高斯積分點(diǎn)，采用帶權(quán)值的拉格朗日插值多項(xiàng)式計(jì)算積分系數(shù)，并通過對(duì)積分核函數(shù)的構(gòu)造實(shí)現(xiàn)快速精確的計(jì)算。下面給出基于高斯-切比雪夫積分的非均勻稀疏分割計(jì)算方法的流程。

1）求n＋1次切比雪夫多項(xiàng)式的n＋1個(gè)零點(diǎn)，用于后續(xù)產(chǎn)生積分節(jié)點(diǎn)。

由于切比雪夫多項(xiàng)式在區(qū)間[-1,1]內(nèi)帶權(quán)值ρ(x)＝(1-x2)-1/2,x∈[-1,1]正交性，n＋1次切比雪夫多項(xiàng)式的n＋1個(gè)零點(diǎn)也是高斯點(diǎn)，對(duì)應(yīng)的高斯-切比雪夫積分的代數(shù)精度是2n＋1次的。切比雪夫多項(xiàng)式滿足下面的遞推關(guān)系式：

其對(duì)應(yīng)的n＋1個(gè)零點(diǎn)為：

2）計(jì)算各節(jié)點(diǎn)對(duì)應(yīng)的帶權(quán)值ρ(x)＝(1-x2)-1/2,x∈[-1,1]的拉格朗日插值基函數(shù)l k(x)，用于計(jì)算各積分節(jié)點(diǎn)對(duì)應(yīng)的積分系數(shù)。x k對(duì)應(yīng)的拉格朗日插值基函數(shù)l k(x)為：

3）求各積分節(jié)點(diǎn)對(duì)應(yīng)的積分系數(shù)A k，計(jì)算式為：

4）積分區(qū)間變換，調(diào)整積分區(qū)間與有限長電性天線長度范圍匹配。

積分節(jié)點(diǎn)x k及其對(duì)應(yīng)積分系數(shù)A k與積分區(qū)間無關(guān)，僅與分點(diǎn)數(shù)（積分階數(shù)）有關(guān)。因此，算法實(shí)現(xiàn)時(shí)可以將常用階數(shù)的積分節(jié)點(diǎn)和積分系數(shù)存儲(chǔ)以備隨時(shí)調(diào)用，以此獲得額外的加速計(jì)算效率。積分區(qū)間變換式如下：

式中，a和b分別為電性天線長度范圍對(duì)應(yīng)的坐標(biāo)。

當(dāng)計(jì)算天線的X軸分量時(shí)，b＝L x/2，a＝-L x/2，計(jì)算得到的為天線X軸分量各分點(diǎn)坐標(biāo)l xk；

當(dāng)計(jì)算天線的Y軸分量時(shí)，b＝L y/2，a＝-L y/2，計(jì)算得到的為天線Y軸分量各分點(diǎn)坐標(biāo)l yk；

當(dāng)計(jì)算天線的Z軸分量時(shí)，b＝L z/2，a＝-L z/2，計(jì)算得到的為天線的Z軸分量各分點(diǎn)坐標(biāo)l zk。求解過程中的x k保持不變。L x、L y和L z通過投影變化得到：

5）代入式（25）中計(jì)算任意姿態(tài)有限長電性天線的輻射電磁場。

本文建立的基于高斯-切比雪夫積分方法進(jìn)行典型層狀海洋模型任意姿態(tài)有限長電性天線電磁場的快速計(jì)算方法，是以n＋1次切比雪夫多項(xiàng)式的零點(diǎn)作為積分節(jié)點(diǎn)構(gòu)建代數(shù)精度為2n＋1次的積分求解方法，實(shí)現(xiàn)了非均勻稀疏分割積分方法，以保證計(jì)算效率和精度。

2.2 模擬計(jì)算結(jié)果與性能分析

為驗(yàn)證本文提出方法的效果，表1給出下列模擬仿真計(jì)算案例，以對(duì)比本文提出方法相比現(xiàn)有計(jì)算方法的性能提升。

表1 模擬仿真案例參數(shù)列表

圖7 —14給出了海底沿x軸和y軸方向0～1 km范圍內(nèi)磁感應(yīng)強(qiáng)度總場、各分量及電場總場和各分量的幅度和相位分布曲線。從幅度分布曲線的計(jì)算結(jié)果可知，電偶極子天線近似計(jì)算結(jié)果誤差最大，甚至短偏移距時(shí)的分布曲線形態(tài)已無法反應(yīng)真實(shí)分布曲線形態(tài)，如圖8和圖12給出的電場總場及各分量幅度分布曲線；從相位分布曲線計(jì)算結(jié)果可知，電偶極子天線近似計(jì)算結(jié)果的相位分布在短偏移距時(shí)存在較大誤差，如圖10和圖14所示。對(duì)比均勻密集分割頻域計(jì)算結(jié)果與本論文提出方法的頻域計(jì)算結(jié)果可知，本文方法的計(jì)算結(jié)果與均勻密集分割計(jì)算結(jié)果精度相當(dāng)，幅度和相位分布曲線具有很好的一致性。

圖7 海底x方向0～1 km范圍內(nèi)磁感應(yīng)總場及各分量的幅度分布曲線

圖8 海底x方向0～1 km范圍內(nèi)電場總場及各分量的幅度分布曲線

圖10 海底x方向0～1 km范圍內(nèi)電場各分量的相位分布曲線

圖12 海底y方向0～1 km范圍內(nèi)電場總場及各分量的幅度分布曲線

圖14 海底y方向0～1 km范圍內(nèi)電場各分量的相位分布曲線

圖15和圖16出了海底（0,500 m）位置的磁感應(yīng)總場、各磁場分量、電場總場和各電場分量的負(fù)階躍響應(yīng)曲線。其中，實(shí)線為電偶極子天線近似結(jié)果，劃線為均勻精細(xì)分割計(jì)算結(jié)果，點(diǎn)線為本文提出的基于高斯-切比雪夫積分方法的非均勻稀疏分割積分計(jì)算結(jié)果。

圖9 海底x方向0～1 km范圍內(nèi)各磁感應(yīng)強(qiáng)度分量的相位分布曲線

圖15 海底(0,500 m)位置的磁感應(yīng)總場及各分量的負(fù)階躍響應(yīng)

圖16 海底(0,500 m)位置的電場總場及各分量的負(fù)階躍響應(yīng)

在計(jì)算機(jī)配置為Win7系統(tǒng)，六核Intel i5-8400，主頻2.80 GHz，8 GB RAM條件下，仿真軟件采用MatlabR2018a-64位版本，模擬計(jì)算單頻點(diǎn)，200個(gè)觀測點(diǎn)6個(gè)分量的仿真計(jì)算耗時(shí)間如表2所示。推薦使用的非均勻稀疏分割計(jì)算參數(shù)為積分節(jié)點(diǎn)數(shù)7～13，誤差限1e-8～1e-12。

表2 計(jì)算耗時(shí)統(tǒng)計(jì) s

圖11 海底y方向0～1 km范圍內(nèi)磁感應(yīng)總場及各分量的幅度分布曲線

對(duì)比本文方法與電偶極子天線近似計(jì)算方法、均勻密集分割計(jì)算方法的計(jì)算結(jié)果精度和計(jì)算效率可知，在保證計(jì)算精度的情況下，本文提出的方法對(duì)層狀海洋有限長電性天線輻射電磁場的時(shí)域和頻域計(jì)算效率的提升達(dá)10倍以上。

圖13 海底y方向0～1 km范圍內(nèi)各磁感應(yīng)強(qiáng)度分量的相位分布曲線

3 結(jié)束語

本文針對(duì)當(dāng)前層狀海洋下電性天線輻射電磁場時(shí)頻域計(jì)算精度和效率無法兼顧的問題，基于高斯-切比雪夫積分方法，構(gòu)建了適用于層狀海洋模型電性天線輻射電磁場的非均勻稀疏分割快速計(jì)算方法，通過采用n＋1次切比雪夫多項(xiàng)式的零點(diǎn)作為分點(diǎn)，實(shí)現(xiàn)將原有的密集均勻分割點(diǎn)數(shù)極大地減小到以切比雪夫多項(xiàng)式的n＋1個(gè)零點(diǎn)為分點(diǎn)的非均勻稀疏分割計(jì)算方案，在保證計(jì)算精度的前提下，實(shí)現(xiàn)了有限長電性天線輻射電磁場的快速計(jì)算，兼顧了計(jì)算的精度和效率。這種非均勻稀疏分割積分方法能夠?qū)崿F(xiàn)工程應(yīng)用中數(shù)據(jù)解釋的效率和準(zhǔn)確性的提升。■