時域有限差分方法分析水下目標散射電場分布

劉 健， 郭元龍， 張慶一， 魏寶君， 劉學鋒， 王殿生

（中國石油大學（華東）理學院，山東 青島 266580）

0 引言

隨著水下電磁波通信和海洋勘探技術的發展，水下目標的電磁探測逐漸成為水下電磁場研究的重點方向。海水的高電導率，中高頻電磁波在海水中傳播具有極強的衰減效應［1］，研究低頻電磁波在水下的傳播和電磁散射問題愈發重要。目前解決水下電磁散射的計算方法主要有時域有限差分法（Finite Difference Time Domain，FDTD）［2］、矩量法［3］、體積分方程法［4］和有限元法［5］等。其中FDTD 方法憑借其計算全空間、全時段和大尺寸目標散射電場分布的快速以及能夠通過傅里葉變換求解全頻域散射結果，已發展成為研究水下電磁散射的常用方法［6］。FDTD 方法是通過將時域麥克斯韋旋度方程進行差分，得到電場強度和磁場強度滿足的時空耦合離散的差分方程，通過迭代求解計算得到全空間電磁場分布［7-8］。在海水電磁場計算中，海水是有耗的色散介質，FDTD方法可在時域下求解色散介質瞬態電磁問題，故FDTD 方法在處理色散介質的電磁衰減并計算水下目標全空間瞬時散射電磁場方面具有優勢。

國內外學者做了大量水下電磁波傳播與散射電磁場分析計算的理論研究。在水下電磁波傳播方面，Xia等［9］利用FDTD方法分析了湖水中低頻電磁波的傳播規律，驗證了FDTD 方法計算水下電磁波傳播的可行性。在電磁散射方面，Luebbers等［10］通過卷積遞歸法差分計算真空中頻率相關材料的三維電磁散射場，為海水等色散介質中FDTD 計算提供了支持；鄭奎松等［11］提出大比例網格變換的總場-散射場源時域FDTD方法，在二維空間計算水下油、氣、水合物等高阻和金屬等低阻物體模型的電磁散射，為利用FDTD方法計算水下電磁散射給予了參考。

由于海水和水下目標是色散介質，入射電磁波頻率f和目標電導率σ 都會對實際水下探測產生影響，需對入射電磁波頻率和水下目標材料電導率進行系統性對比研究?；谏鲜隹紤]，在前人研究的基礎上，構建三維海水色散介質的時域有限差分電磁計算模型。計算電磁波在海水中的傳播，分析低頻電磁波在海水中的衰減規律，并通過與解析解對比，驗證模型在計算海水中電磁波傳播的正確性。進一步建立橢球形水下目標模型，基于三維海水色散介質FDTD 模型計算水下總場和散射電場分布，探究總場和散射電場隨電磁波頻率f的關系。討論電磁波頻率f和目標電導率σ對散射電場分布的影響，為實際水下探測提供參考。

1 三維色散介質中FDTD計算模型

針對建立三維海水色散介質計算模型，利用麥克斯韋旋度方程進行FDTD計算

式中：E、H、D 分別為海水中的電場強度，磁場強度和電位移矢量；μ為海水磁導率。電位移矢量D 和電場強度E在頻域下滿足

式中：ε0為真空介電常數；為復介電常數；ω 為角頻率。海水為有耗散的色散介質，頻域下在色散介質中與頻率有關，

式中：σ為介質電導率；εr為相對介電常數。通過改變式（4）中εr以及σ 可得到某一頻率下介質材料的復介電常數。將式（4）代入式（3）后得到頻域下有耗散的色散介質中電位移矢量：

由于FDTD方法需要把時間差分，得到各時間點的電磁場，依照時間步求解電磁場瞬態變化。利用傅里葉積分變換可把式（5）變為關于時間t的積分，即

基于式（1）～（6），通過FDTD方法便可得到差分方程。結合劃分的電磁場強度空間和時間三維差分網格，即可計算電磁場全空間全時段分布。圖1 所示為海水深度x方向上電場和磁場在空間和時間坐標上的中心差分采樣方法。三維色散介質的FDTD模型中，y和z方向的差分網格與x方向一致。

圖1 計算區域海水深度x方向上空間和時間的中心差分網格

基于圖1，以海水深度x方向為例，取時間離散間隔為Δt，空間離散間隔為Δx，將式（1）、（2）的電場強度E和磁場強度H在時間和空間上進行中心差分，得差分方程為

式中：上標n為計算的時間步長；k為計算的空間步長。

在差分中電場Ez在空間和時間的整數步長取值，磁場Hy在空間和時間的半整數步長取值。計算時間t＝n·Δt，n+1 為向前運行一個時間步，n+1／2 為向前運行半個時間步。在深度x方向計算的空間距離x＝k·Δx，k+1 為運行一個空間步，k+1／2 為運行半個空間步。

式（7）、（8）給出了E和H的差分關系，由于計算區域為海水色散介質，需對時域下的式（6）進行差分采樣。將式（6）中積分近似為對時間步長Δt的求和，給出D和E的關系

從n＝0 時間步開始，引入正弦電場（k），結合式（9）電位移矢量的差分計算式，得（k）。代入式（7）、（8）便可依次求出（k+2）…，即可利用第k層空間步的電場值計算第k+1 層的值，從初始條件開始逐層計算。在n＝1 時間步，同樣引入初始條件（k）和上一時間步的全空間電場強度代入式（7）、（8）便可依次求出（k+2）…。在時間上進行重復迭代便可得到計算區域內各個時間步的全空間水下電場分布。

為在計算區域邊界處電磁波不發生反射，需在邊界處引入完全匹配層（Perfect Matched Layer，PML）［12］。PML基本原理為在計算區域邊界設置一種和計算區域介質完全匹配的介質層，使得電磁波進入PML后不會反射并且完全衰減。PML 和計算區域設置及具體參數如圖2 所示。

圖2 吸收邊界和計算區域示意圖圖中x為海水深度方向，y、z為海平面方向

計算過程中，電磁波由介質1 傳播至介質2 時，電磁波的反射系數Γ 由兩種介質的內在波阻抗決定。為使電磁波在PML 內不發生反射，需在PML 邊界處應保持Γ不變。電磁波在兩介質處反射系數

式中，η1、η2分別為介質1 和介質2 的波阻抗。波阻抗

式中：μ為介質的磁導率；εr為介質的相對介電常數。因此PML層設置磁導率和相對介電常數參數，需滿足阻抗匹配條件［13］：

式中：n為法向；t為切向。

由于電磁波在PML內衰減，式（11）中的相對介電常數和相對磁導率須考慮損耗，即存在虛部

式中：σ（ρ）為PML內的電導率；ρ為到PML內邊界的距離；σm為用導磁率表示的電導率；在阻抗匹配條件式（11）下與σ（ρ）等價。σ（ρ）通常寫成冪指數形式：

式中：σmax＝；d為PML 的厚度；δ 為空間步長；m為常數。本文中將m設為4，確保PML 具有較好的吸收效果。

式（12）在計算中設為邊界兩側的切向電磁參數相同

綜合式（1）、（2）、（13）～（15）可得PML內麥克斯韋方程：

式（16）、（17）與計算區域內的麥克斯韋方程形式相同，因而可選用相同的迭代公式。

利用式（7）～（9）、（16）、（17），結合圖1 劃分的差分網格和圖2 設置的吸收邊界和計算區域便可計算得到全空間、全時段水下電磁場傳播的總場與散射場。

2 水下電磁波傳播散射結果分析

論文利用時域有限差分方法，構建了三維海水色散介質模型，計算水下目標電場強度的總場和散射場分布。設置參數：海水的電導率為σ ＝4 S／m，相對介電常數εr＝80，磁導率μ ＝4π×10-7N／A2［14］。劃分的差分網格選取空間步長δ ＝1 cm，時間步長dt＝。其中，c0為真空中光速。在邊界設置8 個空間步長間隔的PML模擬電磁波的吸收，確保圖2 計算區域內電磁波無反射。

計算中采用正弦點源模擬水下電磁探測，將點源置于海平面下30 m處，點源差分

式中：f為正弦電磁波頻率；t為運行的總時間步數。從初始條件正弦電場源開始逐步計算，在時間和空間上進行迭代得到全空間、全時段水下電磁場分布。

2.1 無散射目標FDTD數值解和解析解對比

為驗證構建的三維海水色散介質中FDTD計算模型的正確性，計算無散射目標時水下電磁波傳播和電場幅值衰減結果，并將FDTD 方法數值解與解析解進行對比。

在不考慮偏振的情況下，色散介質中電磁波所滿足的傳播方程可由式（1）、（2）解析求解：

式中：E0為發射源的初始幅度；α、β 分別為電磁波的衰減常數和相位常數；ηc為介質的本征阻抗。其中［15］：

圖3 不同頻率水下電場強度Ez 在海水深度x方向上FDTD數值解與解析解對比

基于建立的三維色散介質模型，系統計算每個時刻下計算區域內的電場強度及其幅值。在圖3、4 給出頻率分別為5、20 和100 Hz的低頻電磁波在海水中傳播1.25 個周期時，電場強度及其幅值的FDTD數值解和式（19）的解析解對比。

由圖3 可見，基于FDTD 方法得到電場強度數值解與解析解吻合很好，說明FDTD 方法計算數值結果的正確性和可靠性。進一步分析可得，隨頻率的降低，電磁波的波長越長，電磁波在海水色散介質中的衰減越弱，傳播距離越遠，有利于對海水中目標的探測。

由圖3 可見，5 Hz電場強度FDTD 數值解與解析解在深度x＞300 m處有微弱偏差，誤差極值為15.3 mV／m。該偏差的主要影響因素為電場的色散耗散誤差kr，其與傳播時間和劃分的空間網格步長有關［16］。當頻率很低時，電磁波在海水色散介質中的傳播距離很長，傳播時間也就長，造成色散耗散誤差kr在深度x＞300 m處越大。

圖4中給出了3 個頻率下電場強度幅值數值解和解析解對比。由圖中可見，隨著頻率f降低，電磁波的衰減變慢，傳播深度增大。這是由于海水是色散介質，海水中自由電子在高頻強電場的驅動下會在海水中形成電流密度，電場能量轉變為焦耳熱，造成電磁波在海水中傳播的衰減。根據圖4 中電場幅值變化曲線，可得海水中頻率分別為5、20 和100 Hz 的電磁波，其電場強度幅值衰減到初始幅值的1／e時的傳播距離分別為112.54、56.27 和25.16 m。對比圖4 中各頻率電場強度幅值衰減可見，電磁波在海水中的傳播深度與式（21）中α 相關，即與為反比關系［17］。結合圖3、4可得，當頻率低于5 Hz 時差分計算有一定的幅值誤差，會對FDTD模擬水下電磁探測產生影響。而當頻率高于20 Hz時，電場幅值衰減過快不利于水下遠距離探測。低頻電磁波的幅值在海水中衰減較慢，且在探測區域中的電磁場顯著大于洋流、潮波等海水運動切割地磁場產生的感應電磁場［18］，因此可用于水下目標的探測。超低頻電磁波具有傳播穩定、抗干擾能力強等優點，對水下大尺寸目標進行遠距離電磁探測時應選用頻率較低的電磁波，通過探測水下散射電場幅值分布探測水下目標的大小和位置。

圖4 不同頻率水下電場幅值在海水深度x方向上FDTD數值解與解析解對比

2.2 計算水下大尺寸目標散射電場分布

實際探測中，水下目標如潛艇等呈流線型，將水下目標形狀設為旋轉橢球形，研究海水中橢球形目標散射電磁場的分布。旋轉橢球目標長軸為120 m，短軸長為20 m，設置水下目標材料為鋼鐵，電導率σ ＝5 MS／m，相對介電常數εr＝8，相對磁導率μr＝1 000［19］。水下散射目標在計算區域位置如圖2 所示，位于計算區域中心，橢球中心距離海面150 m，構建水下目標模型如圖5 所示。

圖5 水下橢球形目標模型

根據2.1 中分析，選擇頻率為10、20 Hz的正弦點源，構建好模型后利用FDTD方法進行差分，計算旋轉橢球水下電磁場分布，并將x、y切面的電場強度幅值展示在圖6 中。由于全空間電場分布具有較大的數量級差異，引入對數單位dB，其計算式為dB ＝20lg（A／A0），其中A為電場強度或幅值，A0＝1 V／m。

圖6 不同入射電磁波頻率水下目標的電場總場幅值分布

由圖6（a）可知，由于金屬目標散射的影響，不同深度的電場總場及散射場的幅值有顯著差異。在深度30 ～100 m范圍內即目標上方水下總場幅值呈e 指數規律衰減，目標散射電場的空間分布先增大后減小并趨近于0。其原因為水下目標產生的感應電場和散射電場在100 m處相互抵消，目標上方范圍內散射電場起主導作用。在深度100 ～140 m 范圍由于橢球目標的電導率較大，目標對水下電磁波的傳播產生影響，散射電場迅速增大，與入射電場相互抵消，進而導致電場強度總場的幅值迅速衰減。在深度140 ～160 m范圍，即水下目標內部，電磁波在目標表面位置被吸收，幾乎無法傳入目標內部，電場近似為0。這是由于入射到目標表面的電磁波發生了趨膚效應，電磁波在鋼中無法傳播，只能在其表面傳播，在金屬殼體表面產生了感應電流，進而在水下目標的下端產生散射電場，對電磁波傳播起引導作用。當深度大于160 m即目標下方區域，由于目標對電磁波進行了阻擋，總場幅值逐步衰減，感應電流產生的散射電場起主導作用。散射電場隨著深度的增加迅速減小逐漸趨近于0。由圖6 可見，目標上側散射電場比下側強，其原因主要為低頻場源離水下目標的上側較近，電磁波目標對低頻電磁波具有阻擋作用，散射電場主要分布于目標上側。且水下目標的下端產生的散射電場來自于感應電流，其量級小于目標上側。

由圖6（a）可見，20 Hz電磁波在深度60 m處散射電場幅值為-48.41 dB，深度250 m 處電場幅值為-52.82 dB，當前水下電場探測量級可達μV／m，即-120 dB［20］，因此兩處的散射電場均可被水下電磁探測器探測。將水下電磁探測器置于上述海水深度，將低頻正弦點源置于圖2 中所示位置在水下持續發射低頻電磁波，通過分析有散射目標和無散射目標的電場總場幅值變化即散射電場便可精確確定目標位置。

根據圖6（a）、（b）中不同入射頻率下，水下電場幅值分布的對比可以發現，頻率越高，在目標周圍及其下方的電場總場幅值越大。在目標下方90 m 處即深度250 m處10 Hz入射源給出的總場幅值比20 Hz大-10.16 dB。說明10 Hz的電磁波相較于20 Hz探測距離更遠，其水下電場的幅值變化更容易被電磁探測器精確探測。

2.3 入射電磁波頻率對總場及散射場分布的影響

為研究入射電磁波頻率對電場分布及水下電磁探測的影響，基于圖5 水下橢球形目標模型和圖2 設置的計算區域，選擇頻率為10、20 和30 Hz 的入射正弦點源，利用FDTD方法差分計算得到水下電場強度總場和散射場分布，結果分別展示在圖7、8 中，分析不同頻率電磁波的水下電磁探測效果。

圖7 不同入射電磁波頻率在計算區域中軸線（y ＝300 m）處電場總場幅值分布

由圖7 可見，總場幅值大小與入射正弦點源的頻率有關，入射電磁波的頻率越低，總場整體幅值越大。當探測電磁波頻率高于30 Hz時，目標下方的總場幅值小于-80 dB。而海水洋流和潮波運動切割地磁場產生的感應電場數值大約在10-4量級，即為-80 dB左右［21］。海水運動產生的感應電場與30 Hz 入射電磁波的水下總場幅值數量級一致，海水背景運動以及背景磁場電場變化均會極大影響水下電磁探測［22］，在水下電磁探測中應選擇頻率低于30 Hz 的入射電磁波。由圖7 可見，水下總場的整體幅值正比于e-■f，這可由電場強度和頻率的關系式（19）解釋。

圖8 給出了不同頻率電磁波水下散射電場在計算區域軸線上的分布?？梢姡S著深度的增加，散射電場幅值在目標上方范圍呈先增大后減小的趨勢。設旋轉橢球體的高度為單位長度，在目標上方3 ～4 個單位長度處可達到散射電場的極值。在目標下方散射電場以e指數關系隨深度迅速衰減至0。頻率影響散射電場分布特征如下：①頻率越低，海水中的整體散射電場幅值越大，越容易探測目標位置；②頻率越低，目標上方散射電場的極值分布距離潛艇越遠，表明可在更遠處探測散射電場變化，探測更深處的潛艇位置。

圖8 不同入射電磁波頻率在計算區域中軸線（y ＝300 m）處散射電場幅值分布

分析圖7、8，頻率影響總場和散射電場的原因主要有：①電磁波在水下的衰減與頻率有關，頻率低電磁波在水下的衰減弱；②金屬表面產生的感應電流與頻率相關，頻率越低，感應電流越強，在目標上、下方產生的散射電場越強，散射電場的極值分布距離目標越遠。例如將水下電磁探測器放置于目標上方95 m處，10 Hz入射頻率的水下散射電場幅值為-44.65 dB，明顯大于20 Hz的散射電場幅值。因此在水下電磁探測中應盡量選擇頻率低的入射電磁波。

綜合以上可得，入射電磁波的頻率越低，水下目標的總場和散射電場的整體幅值越大，散射電場極值分別距離潛艇越遠，電磁探測效果越好。因此在水下探測時應選用頻率較低的電磁波。

2.4 目標電導率對水下電場分布的影響

水下目標的電場分布不僅與入射電磁波的頻率有關，還與目標的電導率有關。為探究目標電導率對水下電場分布的影響，利用FDTD 方法計算入射電磁波頻率為20 Hz時，水下目標電導率σ分別為102、103和105S／m時的散射電場分布，并將結果展示在圖9中。

圖9 不同電導率水下目標在計算區域中軸線（y ＝300 m）處散射電場幅值分布

圖9給出了在計算區域中軸線（y＝300 m）處不同電導率水下目標的散射電場分布?？梢姡S目標的電導率增加，散射電場的整體幅值也增加，而散射電場的極值分布點保持不變，即目標電導率的變化僅影響散射電場幅值大小而不影響散射電場的分布。這從式（19）～（21）可知，目標電導率僅與電場幅值大小有關。目標材料的電導率越低，對電磁波的吸收也就越弱，表面感應電流產生的散射電場也就越小。

圖10 給出了20 Hz入射電磁波頻率下，不同電導率的水下目標在其上方75 m（散射電場極值）處散射電場FDTD計算幅值及其對數擬合曲線?？芍?，該處散射電場幅值在電導率0 ～2 000 S／m 區間內增幅明顯，隨電導率呈對數增加趨勢。對該散射電場極值處的幅值和電導率進行擬合，即

圖10 20 Hz入射電磁波頻率下不同電導率的水下目標在目標上方75 m（散射電場極值位置）處散射電場幅值FDTD計算結果及對數擬合曲線

得到擬合曲線中的參數分別為：a＝0.231 mV／m、b＝1.21 kS／m。圖10 給出了對數擬合結果和散射電場幅值FDTD計算結果吻合得很好，說明散射電場極值處的幅值正比于ln σ。圖中，當目標電導率大于2 kS／m后，散射電場幅值的增幅逐漸變小直至趨于穩定，這是由于當目標電導率很大時，根據趨膚效應，電磁波已無法傳入金屬內部，各種電導率的水下目標對電磁波傳播的影響一致，在目標表面產生的感應電流基本相同，感應電流產生的散射電場變化趨于穩定，導致不同電導率下水下目標的散射電場趨于穩定。

當前水下目標外殼材料更新較快，在選用高強度鋼滿足下潛承受力的同時還會選擇添加復合材料用于耐海水腐蝕和增強電磁信號的傳播［23］。這會導致目標的電導率降低，也就難利用電磁波精確探測到目標的位置。假設水下電磁探測的極限為5 mV／m，則該電磁探測器將無法精確探測電導率小于105S／m的水下目標。結合圖9、10 可見，水下目標的電導率越大，散射電場在水下的整體幅值越大，也就越容易被電磁探測器精確探測。當電導率大于2 kS／m后，增幅變化逐漸趨于穩定，即電導率的變化不會影響目標散射電場的探測。

3 結語

在海水色散介質中，建立三維時域有限差分電磁計算模型。結合構建水下橢球形散射目標，計算全空間不同時刻下水下目標的電場強度總場和散射場幅值分布，得到水下不同深度處總場和散射電場的變化規律。通過對結果分析發現，入射電磁波頻率和目標電導率對水下探測有直接影響。對于頻率f，水下總場和散射電場的整體幅值正比于，且頻率越低，散射電場極值位置距離目標越遠，電磁探測效果越好。對于目標電導率σ，散射電場幅值在目標電導率0 ～2 000 S／m區間內增幅明顯，后變化趨于穩定，散射電場極值正比于ln σ。本文的研究結果，分析并解釋了頻率和目標電導率對水下目標散射電場的影響，可為后續水下電磁探測應用提供理論參考。