層狀海洋電性天線輻射電磁場高效模擬計算

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

（1.北京機械設備研究所，北京 100854；2.中國科學院空天信息創新研究院，北京 100094）

0 引言

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

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

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

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

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

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

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

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

1.1 典型層狀海洋模型

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

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

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

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

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

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

4）準靜態近似下，認為媒介參數與頻率無關，磁導率和介電常數采用真空中的參數。

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

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

式中，U為標量位，滿足

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中，參數Ci和D i為待定系數，通過理論推導計算結果如下：

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

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

1.2.3 任意姿態電偶極子天線頻域電磁場

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

圖5 層狀海洋模型下任意姿態電偶極子天線

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

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

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

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

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

1.3.1 任意姿態電性天線輻射電磁場的積分計算式

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

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

1.3.2 現有電性天線輻射電磁場的積分計算方法

現有的電性天線輻射電磁場的積分計算方法主要采用電偶極子天線等效方法和均勻密集分割積分方法。對于電偶極子天線等效方法，是將有限長電性天線看作電偶極子天線，此時積分式(18)變為：

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

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

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

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

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

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

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

1）求n＋1次切比雪夫多項式的n＋1個零點，用于后續產生積分節點。

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

其對應的n＋1個零點為：

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

3）求各積分節點對應的積分系數A k，計算式為：

4）積分區間變換，調整積分區間與有限長電性天線長度范圍匹配。

積分節點x k及其對應積分系數A k與積分區間無關，僅與分點數（積分階數）有關。因此，算法實現時可以將常用階數的積分節點和積分系數存儲以備隨時調用，以此獲得額外的加速計算效率。積分區間變換式如下：

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

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

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

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

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

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

2.2 模擬計算結果與性能分析

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

表1 模擬仿真案例參數列表

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

圖7 海底x方向0～1 km范圍內磁感應總場及各分量的幅度分布曲線

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

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

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

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

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

圖9 海底x方向0～1 km范圍內各磁感應強度分量的相位分布曲線

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

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

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

表2 計算耗時統計 s

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

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

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

3 結束語

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