徐龍飛 鄭奎松 張世田 李亞靜

(1.西北工業(yè)大學(xué)電子信息學(xué)院,西安 710072;

2.中國(guó)電波傳播研究所 電波環(huán)境特性及?；夹g(shù)重點(diǎn)實(shí)驗(yàn)室,青島 266107)

引 言

隨著水下目標(biāo)探測(cè)和海洋環(huán)境勘探的快速發(fā)展,水下目標(biāo)的探測(cè)研究逐漸成為人們關(guān)注的焦點(diǎn)[1-2],特別是開發(fā)一種處理水下應(yīng)用的準(zhǔn)確算法是多年來一直的挑戰(zhàn). 在解決電磁散射問題的眾多計(jì)算方法中,時(shí)域有限差分(finite-difference time-domain, FDTD)方法由于其直觀的物理過程和簡(jiǎn)單的計(jì)算格式,已經(jīng)發(fā)展成為計(jì)算水下低頻電磁波傳播的有效工具[3-4]. FDTD方法是將微分形式的麥克斯韋旋度方程進(jìn)行差分離散從而得到一組時(shí)域推進(jìn)公式,并在時(shí)間軸上逐步求解空間的電磁場(chǎng). 對(duì)于一個(gè)三維仿真計(jì)算問題,FDTD方法要求將問題空間離散成許多Yee元胞. 假設(shè)Yee元胞的尺寸為δ,則時(shí)間間隔Δt的大小與元胞尺寸有關(guān),受到Courant穩(wěn)定性條件的限制[5].

對(duì)于傳統(tǒng)的FDTD方法,在計(jì)算仿真低頻電磁波時(shí),需要大量的時(shí)間步才能將仿真計(jì)算完成,計(jì)算時(shí)間非常長(zhǎng). 并行散射場(chǎng)源FDTD(tolal-field scattered-field source FDTD, TSS-FDTD)方法的提出為解決該問題提供了一種途徑,但有些情況下依舊需要非常長(zhǎng)的計(jì)算時(shí)間[6]. 文獻(xiàn)[7]提出了一種在太赫茲頻段利用色散介質(zhì)特性提高計(jì)算效率的新思路. 該方法是通過改變色散介質(zhì)的控制參數(shù)來模擬不同媒質(zhì)材料[7]. 例如,設(shè)置模擬空間填充的背景材料為介電常數(shù)為80的海水,時(shí)間間隔的取值可大約設(shè)置為背景材料為空氣時(shí)的9倍,從而計(jì)算效率可大約提高9倍.

為了解決水下低頻電磁波波長(zhǎng)與殼體目標(biāo)厚度尺寸在電長(zhǎng)度方面相差較大的電磁散射計(jì)算問題,本文提出了一種基于移位算子時(shí)域有限差分方法的加速計(jì)算方法.利用該方法計(jì)算了水下橢球殼體目標(biāo)的電磁散射.數(shù)值結(jié)果表明該方法可以提高計(jì)算效率,縮短計(jì)算時(shí)間,為水下殼體目標(biāo)的探測(cè)、識(shí)別和跟蹤提供一定的工程應(yīng)用價(jià)值.

1 色散介質(zhì)參數(shù)的設(shè)置

本文采用的色散介質(zhì)FDTD方法的控制方程為

(1)

(2)

D(ω)=ε0εr(ω)E(ω)．

(3)

式中:式(1)和(2)是時(shí)域麥克斯韋旋度方程;εr(ω)為相對(duì)介電系數(shù),是頻率的函數(shù),故式(3)表現(xiàn)為頻域的形式.

考慮有耗Drude模型的色散介質(zhì),其相對(duì)介電系數(shù)為

(4)

式中:ωp為等離子體頻率;νc為電子平均碰撞頻率. 按照移位算子方法,式(4)可改寫為

(5)

通過改變式(4)中參數(shù)ε∞、ωp和νc的取值,來模擬出某一頻率下不同的介質(zhì)材料.

2 應(yīng)用算例

2.1 低頻驗(yàn)證

考慮整個(gè)計(jì)算空間填充海水媒質(zhì)的計(jì)算模型,海水媒質(zhì)采用Drude色散介質(zhì)模型來模擬. 在工作頻率f=30 Hz、150 Hz時(shí),式(4)中相關(guān)參數(shù)設(shè)置如下:

νc=1×10-20,ε∞=1×108.

(6)

將式(6)代入式(4)可以得到海水在工作頻率f=30 Hz、150 Hz下的介電常數(shù)都為80. 計(jì)算過程中,空間離散間隔為δ=dx=dy=dz=0.05 m. 時(shí)間間隔ΔtN=N·dx/(2c),其中c為真空中的光速. 在這里,變量N取9 000,即Δt9000=9000·dx/(2c). 為了模擬開域空間,設(shè)置整個(gè)計(jì)算空間被厚度為8個(gè)網(wǎng)格的完全匹配層(convolutional perfect matched layer, CPML)包圍[5].

輻射源設(shè)置為x方向極化的磁偶極子源,加源的表達(dá)式為

(7)

式中,磁偶極矩pinc(t)為

(8)

式(8)中:脈寬τ=1.5 ms;頻率ωi=2πi,i=1,2,3;系數(shù)β為初始振幅,在這里β=1.

2.1.1 工作頻率f=30 Hz

整個(gè)計(jì)算空間為長(zhǎng)方體,其物理尺寸取為1000 m×2 m×3 m. 偶極子輻射場(chǎng)中Hx和Ez分量在yoz面上振幅相位的快照如圖1所示.

(a) Hx振幅 (b) Hx相位(a) Ampulitude of Hx (b) Phase of Hx

(c) Ez振幅 (d) Ez相位(c) Ampulitude of Ex (d) Phase of Ex 圖1 f=30 Hz時(shí)yoz面上振幅相位的快照Fig.1 The amplitude and phase distributions in the yoz plane at the working frequency f=30 Hz

觀察圖1中Hx及Ez的相位快照可以看出,在海水中工作頻率f=30 Hz的波長(zhǎng)大約為λ=385 m. 根據(jù)電磁波理論,這個(gè)頻率的電磁波在海水中波長(zhǎng)的理論解大約為333 m. 這兩個(gè)結(jié)果有差別是由于根據(jù)圖1相位變化的方法測(cè)得的波長(zhǎng),觀測(cè)位置離偶極子源很近,沒有位于偶極子輻射場(chǎng)的遠(yuǎn)區(qū)區(qū)域,這就造成與利用平面波特性方法計(jì)算得到的理論解結(jié)果不一致. 如果圖1所示的相位變化法的觀測(cè)位置是離偶極子源更遠(yuǎn)的位置,這兩個(gè)值之間的差別會(huì)逐漸縮小. 這一原因也通過計(jì)算工作頻率f=150 Hz在海水中波長(zhǎng)的算例得到了驗(yàn)證.

為了驗(yàn)證偶極子產(chǎn)生場(chǎng)振幅的變化情況,根據(jù)電磁波理論給出了電偶極子在有耗媒質(zhì)中產(chǎn)生場(chǎng)的解析解為

(9)

(a) 歸一化Hx幅值 (b) 相對(duì)誤差(a) Normalized ampulitude of Hx (b) Relative error 圖2 解析解和計(jì)算解計(jì)算結(jié)果的對(duì)比Fig.2 Comparison between analytical and computational solutions

2.1.2 工作頻率f=150 Hz

本算例中,計(jì)算空間的物理尺寸為450 m×4 m×4 m. 偶極子輻射場(chǎng)中Hx和Ez分量在yoz面上的振幅相位快照如圖3所示.

(a) Hx振幅 (b) Hx相位(a) Ampulitude of Hx (b) Phase of Hx

(c) Ez振幅 (d) Ez相位(c) Ampulitude of Ex (d) Phase of Ex 圖3 f=150 Hz時(shí)yoz面上振幅相位的快照Fig.3 The amplitude and phase distributions in the yoz plane at the working frequency f=150 Hz

同理,根據(jù)電磁波理論,海水中工作頻率f=150 Hz波長(zhǎng)的理論解約為λ=149 m. 根據(jù)相位變化法從圖3中可以讀出波長(zhǎng)的計(jì)算解約為λ=157 m. 通過對(duì)比圖1中的結(jié)果,可以看出隨著工作頻率的升高偶極子產(chǎn)生的場(chǎng)會(huì)更快地到達(dá)遠(yuǎn)區(qū)場(chǎng). 圖4給出yoz面上輻射場(chǎng)沿y軸的變化情況.

將計(jì)算所得結(jié)果進(jìn)行歸一化,與解析解結(jié)果進(jìn)行對(duì)比,結(jié)果及相對(duì)誤差如圖4所示. 通過觀察圖2(b)和圖4(b)可以看出,計(jì)算結(jié)果和理論結(jié)果吻合.

(a) 歸一化Hx幅值 (b) 相對(duì)誤差(a) Normalized ampulitude of Hx (b) Relative error 圖4 解析解與歸一化后計(jì)算結(jié)果的對(duì)比Fig.4 Comparison between analytical and computational solutions

2.1.3 加速性能的測(cè)試

為了測(cè)試本文方法的加速計(jì)算性能,設(shè)置了本算例.計(jì)算環(huán)境為惠普Z620工作站(CPU E5-2620 v2 2.10 GHz,內(nèi)存32 GB). 計(jì)算空間的網(wǎng)格尺寸為30δ×30δ×30δ,時(shí)間間隔為ΔtN=Nδ/(2c),變量N取正整數(shù). 在本算例中N取三個(gè)數(shù)值,即N=1,5,7.加速比Sp和加速效率Ep的定義與文獻(xiàn)[8]中的定義相同,計(jì)算結(jié)果如表1所示.

表1 加速性能的測(cè)試結(jié)果Tab.1 Test results of the accelerated performance

在測(cè)試程序加速性能算例中,要使三次計(jì)算的總時(shí)間保持一致,即三次計(jì)算需滿足M×Δt一樣. 通過觀察表1可知,隨著N的增加,加速性能接近線性變化,加速效率達(dá)到約95%以上.

2.2 金屬殼體目標(biāo)的電磁散射

2.2.1 起伏海面的建模及影響

運(yùn)用Monte Carlo方法建立起伏海面[9],基本原理簡(jiǎn)述如下:假設(shè)x、y方向的長(zhǎng)度分別為L(zhǎng)x、Ly,兩個(gè)方向的采樣點(diǎn)數(shù)分別為M、N,采樣間距分別為Δx和Δy,則起伏面上每一點(diǎn)(xm,yn)(xm=mΔx,yn=nΔy)的高度表示為

F(km,kn)exp[j(kmxm+knyn)].

(10)

式中,

(11)

式(11)中:S(km,kn)為起伏海面的功率譜;km=2πm/Lx;kn=2πn/Ly;二維高斯粗糙表面的功率譜為

(12)

式(12)中:δ代表粗糙面的均方根高度;lx和ly代表沿x軸和y軸的相關(guān)長(zhǎng)度.

考慮海平面為三級(jí)海況時(shí),可令粗糙面的均方根高度δ=0.8 m,相關(guān)長(zhǎng)度為lx=ly=0.2 m,得到海面如圖5(a)所示;保持均方根高度不變,改變其相關(guān)長(zhǎng)度為lx=ly=0.15 m,得到起伏海面如圖5(b)所示.

將圖5所示的起伏海面模型分別疊加到圖6所示的模型中,對(duì)頻率f=150 Hz的低頻場(chǎng)進(jìn)行計(jì)算. 為了查看起伏海面對(duì)低頻電磁波的影響,將上述兩種起伏海面的計(jì)算結(jié)果做差值處理,得到海面對(duì)低頻電磁波的影響,如圖7所示. 觀察圖7海面可以看出起伏海面影響低頻電磁場(chǎng)的空間分布,特別是對(duì)電場(chǎng)分布的影響.

(a) δ=0.8 m,lx=ly=0.2 m

(b) δ=0.8 m,lx=ly=0.15 m圖5 起伏海面示意圖Fig.5 Schematic diagram of rough sea surface

圖6 計(jì)算起伏海面下低頻電磁波傳播的三維模型Fig.6 3D model for calculating the low frequency electromagnetic wave under rough sea surface

(a) Hx振幅差 (b) Ez振幅差(a) Amplitude difference of Hx (b) Amplitude difference of Ex圖7 起伏海面對(duì)低頻電磁波的影響Fig.7 Impact of rough sea surface on the low-frequency electromagnetic waves

2.2.2 水下金屬殼體目標(biāo)

圖8給出了計(jì)算平滑海面下低頻電磁波傳播的三維模型.在z方向依次建立空氣和海水兩種分層媒質(zhì),其中空氣媒質(zhì)厚度為2.5 m,海水媒質(zhì)厚度為6 m. 整個(gè)計(jì)算空間的物理尺寸為21 m×6 m×8.5 m. 建立如圖8所示的金屬橢球狀殼體目標(biāo),其半長(zhǎng)軸為7.5 m,半短軸為1 m. 殼體目標(biāo)中心點(diǎn)距離海面高度為3 m,殼體厚度為0.1 m,殼體為金屬材質(zhì),殼體內(nèi)部填充空氣.

圖8 計(jì)算平滑海面下低頻電磁波傳播的三維模型Fig.8 3D model for calculating the low frequency electromagnetic wave under the smooth sea surface

整個(gè)計(jì)算空間的媒質(zhì)材料都采用有耗Drude模型進(jìn)行模擬. 激勵(lì)源為x方向極化的磁偶極子源,如式(8)所示. 在工作頻率f=150 Hz時(shí),空氣、海水和金屬的相關(guān)設(shè)計(jì)參數(shù)為:

空氣

vc=1×10-18,ε∞=1×1010;

(13)

海水

vc=1×10-18,ε∞=1×1010;

(14)

金屬

vc=1×10-18,ε∞=1×1010.

(15)

經(jīng)過計(jì)算得到的電磁場(chǎng)各分量在yoz面上的振幅快照如圖9所示.

(a) Hx振幅 (b) Hy振幅(a) Ampulitude of Hx (b) Ampulitude of Hy

(c) Hz振幅 (d) Ex振幅(c) Ampulitude of Hz (d) Ampulitude of Ex

(e) Ey振幅 (f) Ez振幅(e) Ampulitude of Ey (f) Ampulitude of Ez圖9 低頻電磁場(chǎng)各分量在yoz面上的振幅快照Fig.9 The amplitude distributions in the yoz plane

觀察圖9可以看出,殼體目標(biāo)影響低頻電磁波的傳播. 由于殼體目標(biāo)是金屬導(dǎo)體,所以對(duì)低頻電磁波的傳播起到了引導(dǎo)的作用. 由于工作頻率為f=150 Hz,通過計(jì)算發(fā)現(xiàn)低頻電磁波透過殼體進(jìn)入殼體內(nèi)部傳播,但是量級(jí)很小. 海面影響低頻電磁波各分量的空間分布,特別是電場(chǎng)分量,這是由于海水的電導(dǎo)率比空氣中的電導(dǎo)率要大. 觀察圖9(f)還可以得到,當(dāng)?shù)皖l電磁波穿過海面到達(dá)空氣中傳播時(shí),其傳播速度要比在海水中的傳播速度快,特別是對(duì)電場(chǎng)的空間分布影響比較大,從圖中可以很明顯地看出海面的空間位置. 由圖9(d)可知,由于水下目標(biāo)的存在導(dǎo)致在海水中存在電場(chǎng)Ex分量.

為了考慮水下目標(biāo)散射場(chǎng)的空間分布,將圖9的計(jì)算結(jié)果重新作了處理,即將存在橢球體目標(biāo)計(jì)算的結(jié)果減去無橢球體目標(biāo)的計(jì)算結(jié)果,具體結(jié)果如圖10所示.

(a) Hx振幅 (b) Hy振幅(a) Ampulitude of Hx (b) Ampulitude of Hy

(c) Hz振幅 (d) Ex振幅(c) Ampulitude of Hz (d) Ampulitude of Ex

(e) Ey振幅 (f) Ez振幅(e) Ampulitude of Ey (f) Ampulitude of Ez圖10 水下目標(biāo)在yoz面上的低頻散射場(chǎng)Fig.10 The amplitude distributions in the yoz plane

觀察圖10可以明顯地看出,在水下殼體目標(biāo)周圍存在散射場(chǎng),且殼體目標(biāo)左側(cè)的散射場(chǎng)要比右側(cè)散射場(chǎng)強(qiáng),這是由于低頻源離殼體目標(biāo)的左側(cè)較近. 由于低頻電磁源的照射,在金屬殼體表面產(chǎn)生了感應(yīng)電流,從而在殼體目標(biāo)的右端點(diǎn)會(huì)產(chǎn)生較強(qiáng)的散射場(chǎng),例如圖10(b)的Hy分量和圖10(e)的Ey分量. 由于海水電導(dǎo)率的存在,導(dǎo)致電場(chǎng)的衰減要比磁場(chǎng)的衰減快,這可以從圖10(a)～(c)與圖10(d)～(f)的對(duì)比中得出. 顯然,海面也影響散射場(chǎng)各分量的空間分布,特別是對(duì)電場(chǎng)各分量的影響比較明顯,而且散射場(chǎng)在海面上方的傳播速度要比在海水中的傳播速度快.

3 結(jié) 論

利用本文方法計(jì)算了水下殼體目標(biāo)的散射場(chǎng)分布,并對(duì)該方法的加速性能進(jìn)行了測(cè)試.計(jì)算結(jié)果表明海面的存在可以影響殼體目標(biāo)散射場(chǎng)的空間分布規(guī)律,特別是對(duì)電場(chǎng)空間分布影響比較大,并且不同的磁分量具有不同的空間分布規(guī)律.本文的計(jì)算結(jié)果可以為水下殼體目標(biāo)的探測(cè)、識(shí)別和跟蹤提供參考.