時宗洋，趙一宇，渠曉東，馬力超

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

0 引言

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

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

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

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

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

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

本文首先介紹典型層狀海洋模型下電性天線輻射電磁場的工程應(yīng)用及現(xiàn)階段主要計算方法和存在的不足，引出本文建立的非均勻稀疏分割積分方法及其優(yōu)勢；然后構(gòu)建了層狀海洋模型電磁天線輻射電磁場模型，并對電磁場積分計算式進行了理論推導(dǎo)，并且建立基于高斯切比雪夫分點的非均勻稀疏分割積分方法，給出理論計算式和算法流程；接下來通過模擬仿真算例，對電偶極子等效計算方法、密集均勻分割積分方法和本文提出的非均勻稀疏分割積分這3種方法在電磁場計算精度和效率方面進行對比分析，論證了本文提出方法在兼顧計算精度和效率方面的優(yōu)勢，能夠?qū)崿F(xiàn)工程應(yīng)用中數(shù)據(jù)解釋的效率和準確性的提升；最后對全文進行總結(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，相對磁導(dǎo)率μr設(shè)置為1。AB表示電性天線首尾，天線長度為L，天線電極矩為P＝IL，天線中點位于原點O正下方，目標所在深度為H，距離海底深度為h，各層厚度為d i。笛卡爾坐標系正z方向垂直水平面向下，o-xyz滿足右手螺旋定則。假設(shè)任意姿態(tài)電性天線在XOY平面內(nèi)的投影與x軸夾角為θ，天線與XOY平面的夾角為?。

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

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

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

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

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

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

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

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

式中，U為標量位，滿足

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

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

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

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

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

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

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

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

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

海水中電場和磁場各分量頻域計算式見式(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)計算結(jié)果如下：

式中，N01、N21、P01和P21計算式如下：

式中，R*2和Q*2的遞推計算式為：

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

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

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

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

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

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

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

式中，P01和P21計算式如下：

式中，R*2的遞推計算式如下：

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

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

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

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

1.2.4 電偶極子天線電磁場計算方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其對應(yīng)的n＋1個零點為：

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

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

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

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

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

當(dāng)計算天線的X軸分量時，b＝L x/2，a＝-L x/2，計算得到的為天線X軸分量各分點坐標l xk；

當(dāng)計算天線的Y軸分量時，b＝L y/2，a＝-L y/2，計算得到的為天線Y軸分量各分點坐標l yk；

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

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

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

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

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

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

圖7 —14給出了海底沿x軸和y軸方向0～1 km范圍內(nèi)磁感應(yīng)強度總場、各分量及電場總場和各分量的幅度和相位分布曲線。從幅度分布曲線的計算結(jié)果可知，電偶極子天線近似計算結(jié)果誤差最大，甚至短偏移距時的分布曲線形態(tài)已無法反應(yīng)真實分布曲線形態(tài)，如圖8和圖12給出的電場總場及各分量幅度分布曲線；從相位分布曲線計算結(jié)果可知，電偶極子天線近似計算結(jié)果的相位分布在短偏移距時存在較大誤差，如圖10和圖14所示。對比均勻密集分割頻域計算結(jié)果與本論文提出方法的頻域計算結(jié)果可知，本文方法的計算結(jié)果與均勻密集分割計算結(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)總場、各磁場分量、電場總場和各電場分量的負階躍響應(yīng)曲線。其中，實線為電偶極子天線近似結(jié)果，劃線為均勻精細分割計算結(jié)果，點線為本文提出的基于高斯-切比雪夫積分方法的非均勻稀疏分割積分計算結(jié)果。

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

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

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

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

表2 計算耗時統(tǒng)計 s

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

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

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

3 結(jié)束語

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