馮世杰

（電子科技大學(xué)電子科學(xué)技術(shù)研究院，四川成都 611731）

0 引 言

臨近空間（距離地面20～100 km 高度的區(qū)域）內(nèi)的高超聲速目標(biāo)（如彈道導(dǎo)彈、航天飛機(jī)）再入或在臨近空間中飛行時會與空氣產(chǎn)生強(qiáng)烈摩擦，導(dǎo)致目標(biāo)周圍空氣發(fā)生電離，形成包裹目標(biāo)的等離子體鞘套及目標(biāo)身后范圍較大、強(qiáng)度較低的等離子尾跡，再入過程中的通信黑障便是由此產(chǎn)生［1-2］。洲際導(dǎo)彈等再入武器的等離子體鞘套會改變再入目標(biāo)的雷達(dá)目標(biāo)特性，從而對有源雷達(dá)的目標(biāo)探測帶來了一定欺騙和干擾作用。

等離子體鞘套及尾跡可能會對在其中傳播的電磁波產(chǎn)生衰減、折射、反射等效應(yīng)［3-4］。由于時域有限差分方法（Finite-Difference Time-Domain,FDTD）在色散介質(zhì)的仿真計算中存在的獨(dú)特優(yōu)越性，被廣泛用于等離子體的電磁模擬［5］。其中移位算子時域有限差分方法［6］因編程簡單、概念清晰且計算精度高而受到較多關(guān)注［7-11］。但這些研究中，普遍使用二維鞘套模型，少有人對等離子鞘套進(jìn)行三維建模仿真。

等離子體的電磁特性不利于有源雷達(dá)的目標(biāo)探測，但可用于無源雷達(dá)對空間目標(biāo)的探測。無源雷達(dá)，即自身不具備信號發(fā)射器，利用空間已有輻射源進(jìn)行目標(biāo)探測的雷達(dá)?；贕NSS 的無源雷達(dá)因具有成本低、信號源廣泛、覆蓋率高等特點，受到的關(guān)注日益增多。有研究發(fā)現(xiàn)，數(shù)類不同飛行器的飛行過程均會在短期內(nèi)降低經(jīng)過目標(biāo)飛行路徑上的GNSS 信號強(qiáng)度［12］。我國也有不少學(xué)者，在嘗試?yán)肎NSS 信號來進(jìn)行空間目標(biāo)探測［13-16］。但目前少有學(xué)者就高超聲速目標(biāo)尾跡對GNSS信號的衰減特性進(jìn)行研究。

由上可見，針對等離子體鞘套電磁特性的數(shù)值仿真研究中大多使用二維模型，對鞘套內(nèi)等離子體的三維分布重視不足。同時鞘套尾跡范圍較大，對GNSS 信號的衰減作用明顯，這些特征均十分適用于被動式目標(biāo)探測，但現(xiàn)有研究對鞘套尾跡的關(guān)注較少。本文借助計算流體動力學(xué)（Computational Fluid Dynamics,CFD），模擬了高超聲速再入體在臨近空間飛行時的三維等離子體流場，建立了三維等離子體鞘套電磁模型，并使用改進(jìn)的SO-FDTD算法，仿真計算了GNSS信號在目標(biāo)不同飛行條件下的等離子尾跡中的功率透射系數(shù)，繪制了透射系數(shù)分布圖。為基于GNSS 的被動式空間目標(biāo)探測提供了一定依據(jù)。

1 目標(biāo)繞流流場模型

圖1所示為文中建立的高超聲速目標(biāo)模型，其中模型總長0.9 m，直徑0.09 m，尾翼寬度0.045 m。在ESI_Group 軟件中，使用基于熱化學(xué)非平衡效應(yīng)的七組元空氣化學(xué)反應(yīng)模型［17］和Park 雙溫度模型［18］，以全N-S方程組對再入體熱化學(xué)非平衡流場進(jìn)行求解，仿真計算不同飛行條件下的目標(biāo)繞流流場。不同飛行條件下的背景大氣密度、壓強(qiáng)和溫度等相關(guān)參量，參考美國標(biāo)準(zhǔn)大氣模型（1976）進(jìn)行取值。

圖1 目標(biāo)模型示意圖

等離子角頻率以及等離子碰撞頻率是等離子體鞘套的兩個重要參數(shù)，能夠顯著影響等離子體鞘套的電磁特性。由文獻(xiàn)［4］，可基于仿真得出的流場參數(shù)計算等離子角頻率和等離子碰撞頻率，計算公式如式（1）、（2）所示：

式中：ωp為等離子角頻率；ven為等離子碰撞頻率；ne和nm分別為鞘套中的電子數(shù)密度以及中性粒子數(shù)密度；T為氣體溫度；me為電子質(zhì)量，me=9.1×10-31kg；e為電子電荷量，e=1.6×10-19C；ε0為真空介電常數(shù)，ε0=8.85×10-12F/m。

利用式（1）、（2）對高超聲速再入體繞流流場仿真結(jié)果進(jìn)行計算，可得等離子體鞘套的電磁參數(shù)分布。繞流流場隨目標(biāo)飛行高度、速度、攻角等條件的變化而存在區(qū)別，但流場結(jié)構(gòu)具有一定相似性，下面僅以海拔高度H=40 km、目標(biāo)飛行速度V=20 Ma、攻角為0°的繞流流場為例進(jìn)行說明。

圖2、圖3所示為再入模型繞流流場的參數(shù)分布。由圖2和圖3可知，等離子角頻率及等離子碰撞頻率在繞流流場頭部達(dá)到最大值，在流場身部及尾跡區(qū)域顯著下降。繞流流場在目標(biāo)身后不斷擴(kuò)散，尾跡外側(cè)的等離子角頻率低于尾跡中心，同時尾跡外側(cè)的等離子碰撞頻率高于尾跡中心。

圖2 Ma=20，H=40 km條件下的等離子角頻率分布圖

圖3 Ma=20，H=40 km條件下的等離子碰撞頻率分布圖

2 電磁波傳播的計算方法及模型

2.1 改進(jìn)的SO-FDTD方法

本文通過移位算子時域有限差分方法（SOFDTD），對電磁波與等離子體鞘套的相互作用進(jìn)行計算分析。在碰撞非磁化冷等離子介質(zhì)中，Maxwell方程組及其輔助方程為

式中：J為極化電流密度；E為電場強(qiáng)度；H為磁場強(qiáng)度；ωp和νen分別為等離子角頻率和等離子碰撞頻率；μ0為真空磁導(dǎo)率，μ0=4π×10-7N·A-2。

假定J與E的離散空間坐標(biāo)相同，E值更新于整數(shù)時間步，而J和H的值更新于半整數(shù)時間步，以x分量為例，可得式（3）、（4）兩個旋度方程的FDTD差分離散格式為

引入移位算子［6］：

利用移位算子對式（5）進(jìn)行化簡得

式中：p0=p1=2ε0ω2p；q0=νen+；q1=2νen；q2=νen-2/Δt。

如果直接使用式（9）進(jìn)行迭代計算，則程序內(nèi)存開銷較大，因為需要引入Jxn-1/2，Jxn-3/2，Enx和Enx-1四個輔助變量存儲以前時刻的場分量。這里采用文獻(xiàn)［8］中給出的內(nèi)存優(yōu)化算法，引入一個輔助變量Jx1，這樣式（9）可由式（10）、式（11）所代替：

可見，利用式（10）、式（11）進(jìn)行迭代計算，僅需使用1 個輔助變量Jx1即可，這種計算方式能夠節(jié)約內(nèi)存，特別是在三維模型中。但需注意的是，在計算式（10）之前，需用一個臨時變量保存Jn-1/2x，以防在計算式（11）時變量被覆蓋丟失。

2.2 電磁仿真模型的建立

SO-FDTD仿真程序中采用邊長為0.007 5 m的立方體網(wǎng)格對等離子體流場進(jìn)行建模，仿真程序通過Matlab 編寫。SO-FDTD 電磁仿真模型整體為三維矩形，其中再入體模型的軸向平行于電磁仿真模型中的x軸，且再入體模型的軸向位置位于電磁仿真模型y-z平面的中心。因本文主要關(guān)注流場尾跡導(dǎo)致的電磁波衰減，所以電磁仿真模型在x軸方向上的仿真范圍較大，為7 m，在y軸及z軸方向上仿真范圍均為1.7 m。電磁仿真模型內(nèi)共有49 916 846個網(wǎng)格，模型邊界使用了10層CPML 邊界。為了將CFD 計算得到的流場參數(shù)轉(zhuǎn)換到SOFDTD 電磁仿真模型中，本文使用距離反比法進(jìn)行了空間差值運(yùn)算。由于仿真實驗中的模型范圍較大，網(wǎng)格數(shù)較多，為加快程序運(yùn)行速度，在距離反比差值以及FDTD 迭代程序中均使用了GPU 并行加速技術(shù)。

在高度H=40 km，馬赫數(shù)Ma=20，攻角為0°的飛行條件下，頻率為1 575.42 MHz 的平面電磁波沿z軸正方向傳播時，電磁波在繞流流場中的歸一化電場幅度分布如圖4所示。從圖4中可以發(fā)現(xiàn)，在目標(biāo)中軸線所處的x-y平面上，目標(biāo)頭身部區(qū)域附近的電場幅度分布出現(xiàn)了條紋狀的“波紋”。這是由于高超聲速目標(biāo)頭身部區(qū)域激波的存在，因此目標(biāo)頭身部區(qū)域的介質(zhì)參數(shù)分布存在較大梯度，從而導(dǎo)致了入射電磁波的反射及折射現(xiàn)象。

圖4 平面波在繞流流場中的歸一化電場幅度分布示意圖

圖5為模型中目標(biāo)頭身部x=0.75 m 處的y-z平面上的電場分布示意圖。從圖5可見，反射及折射后改變傳播方向的電磁波又與原傳播路徑上的電磁波發(fā)生干涉，從而導(dǎo)致了圖4中出現(xiàn)的“波紋”。

圖5 目標(biāo)頭身部y-z平面切片上的電場幅度分布示意圖

3 等離子體流場尾跡中的GNSS信號衰減特性

在仿真程序中，設(shè)平面電磁波沿模型z軸正方向傳播。由于鞘套尾跡的范圍較大，對GNSS 信號的吸收衰減更容易被探測發(fā)現(xiàn)，所以本文重點關(guān)注鞘套尾跡對信號的影響。圖4所示結(jié)果的同等仿真條件下，模型中部x=3.50 m 處尾跡區(qū)域的y-z平面上的歸一化電場分布圖如圖6所示。

圖6 目標(biāo)尾跡中y-z平面切片上的電場幅度分布示意圖

由圖6可知，電磁波信號幅度在穿過等離子體鞘套后存在較明顯的衰減，電場幅度的周期性變化在模型中z＞1.40 m 的區(qū)域內(nèi)相對穩(wěn)定。所以在本文的后續(xù)實驗中，選擇在z=1.60 m 處的x-y平面上對透射信號進(jìn)行取樣，該平面可稱為取樣面。信號經(jīng)過等離子體流場后的衰減由歸一化功率透射系數(shù)Pt量化，Pt由式（12）確定：

式中，Ei為入射波電場強(qiáng)度，Et為透射波電場強(qiáng)度。

3.1 信號頻率變化的影響

設(shè)再入體飛行條件為：高度H=40 km、目標(biāo)飛行速度V=20 Ma，攻角為0°。入射波頻率f1設(shè)為北斗B2b 信號和伽利略E5b 信號所用的載波頻率1 207.14 MHz；入射波頻率f2設(shè)為北斗B1c 信號、GPS L1 信號和伽利略E1 信號所用的載波頻率1 575.42 MHz。仿真得到的信號透射系數(shù)在取樣面上的分布圖如圖7所示。

圖7 不同頻率GNSS信號的透射系數(shù)分布圖

由圖7可知，在本例給出的仿真條件下，GNSS信號經(jīng)過高超聲速再入體尾跡后的衰減較為明顯，尾跡中心區(qū)域的透射系數(shù)在0.2 以下。同時在目標(biāo)頭身部對應(yīng)的區(qū)域，透射系數(shù)的分布出現(xiàn)了3條相鄰的凹陷區(qū)域，由圖5可知這一現(xiàn)象同樣由電磁波的干涉所致。

由圖7（a）、（b）兩圖的對比可知，相同飛行條件下，頻率更低的GNSS 信號經(jīng)過鞘套尾跡后的透射系數(shù)更低，信號衰減更為嚴(yán)重。這是等離子體的高通濾波特性所帶來的結(jié)果，只有入射電磁波的角頻率大于等離子角頻率時，入射波才能在等離子體中透射。

3.2 飛行速度變化的影響

設(shè)入射GNSS信號的載波頻率為f=1 575.42 MHz，目標(biāo)飛行高度H=40 km，攻角為0°，目標(biāo)飛行速度V分別設(shè)為18、20、22 Ma，取樣面上的透射系數(shù)分布如圖8所示。

圖8 不同飛行速度下GNSS信號的透射系數(shù)分布圖

由圖8可知，在其他條件相同時，尾跡中的信號透射系數(shù)隨著飛行速度的增大而呈現(xiàn)整體下降的趨勢，但它們之間并非簡單的線性關(guān)系，透射系數(shù)的分布結(jié)構(gòu)會隨著飛行速度的增加而發(fā)生變化。在飛行速度為18 Ma 和20 Ma 時，取樣面尾跡中心的信號透射系數(shù)低于尾跡邊緣；而在飛行速度為22 Ma 時，情況相反，尾跡中部區(qū)域的信號透射系數(shù)不僅高于尾跡頭部及邊緣區(qū)域，甚至超過20 Ma飛行速度下相同區(qū)域的透射系數(shù)值。

上述現(xiàn)象，是不同運(yùn)動速度條件下的流場參數(shù)分布差異所導(dǎo)致的。在相對低速（18 Ma、20 Ma）條件下，尾跡邊界與中心的等離子碰撞頻率差異不大，所以等離子角頻率的分布決定了信號的衰減程度。由于尾跡中心的等離子角頻率高于尾跡邊緣，所以尾跡中心流場對透射信號的衰減作用更明顯。但是在相對高速（22 Ma）條件下，流場尾跡邊緣及頭部的等離子碰撞頻率加劇，由式（10）可知這些區(qū)域中等離子體對信號的衰減也隨之增加。同時尾跡邊界處碰撞頻率的增加也一定程度上加快了氣體分子復(fù)合速度，降低了尾跡中心區(qū)域的電子含量和等離子角頻率，從而導(dǎo)致了尾跡中后段信號透射系數(shù)的回升。

3.3 飛行高度變化的影響

設(shè)入射GNSS信號的載波頻率為f=1 575.42 MHz，飛行馬赫數(shù)Ma設(shè)為20，攻角為0°，再入體飛行高度H分別設(shè)為35、40、45 km，不同高度下的大氣參數(shù)如表1所示。

表1 不同高度下的大氣參數(shù)表

不同高度下，GNSS 信號透射系數(shù)仿真計算結(jié)果如圖9所示。

圖9 不同飛行高度下GNSS信號的透射系數(shù)分布圖

由圖9可知，與3.2節(jié)中的現(xiàn)象相似，在飛行速度一定的情況下，尾跡中的信號透射系數(shù)隨著飛行高度的降低而呈現(xiàn)整體下降的趨勢，因為隨著海拔高度的降低，鞘套中氣體分子經(jīng)高溫電離產(chǎn)生的等離子體濃度也隨之升高。但透射系數(shù)的分布結(jié)構(gòu)會隨著飛行高度的降低而發(fā)生變化。在飛行高度H為35 km 時，取樣面尾跡中心區(qū)域的信號透射系數(shù)高于尾跡邊緣及頭部區(qū)域；而在飛行高度為40 km 和45 km 時，尾跡中心區(qū)域的信號透射系數(shù)低于尾跡邊緣及頭部。這與3.2 節(jié)中出現(xiàn)的類似現(xiàn)象有著相同的原因，即更低的飛行高度條件下，尾跡頭部及邊界處的等離子碰撞頻率加劇，增強(qiáng)了這些區(qū)域?qū)﹄姶挪ǖ奈眨布涌炝宋槽E中氣體分子的復(fù)合速度，從而導(dǎo)致了信號透射系數(shù)分布結(jié)構(gòu)的變化。

4 結(jié)束語

本文通過熱化學(xué)非平衡模型模擬了導(dǎo)彈型再入體在臨近空間飛行的繞流流場，并利用改進(jìn)的SO-FDTD 方法對GNSS 信號在流場尾跡中的衰減特性作了計算和分析，得出了以下結(jié)論：

1）臨近空間高超聲速目標(biāo)繞流流場尾跡能夠給GNSS 信號帶來明顯的衰減。衰減特性因再入高度、速度、信號頻率和取樣點位置的不同而存在區(qū)別。

2）目標(biāo)飛行條件一定時，頻率越高的GNSS信號在繞流流場中的透射系數(shù)越高。

3）在其他飛行條件相同的情況下，高超聲速目標(biāo)運(yùn)動速度增加或高度降低時，GNSS 信號透射系數(shù)呈現(xiàn)整體降低的趨勢。但飛行條件變量與透射系數(shù)之間的關(guān)系并非簡單的線性關(guān)系，透射系數(shù)分布結(jié)構(gòu)會隨著尾跡流場參數(shù)的變化而發(fā)生變化。

4）一般情況下，高超聲速目標(biāo)尾跡中部區(qū)域?qū)π盘柕乃p效應(yīng)強(qiáng)于尾跡邊界及頭部區(qū)域。但當(dāng)其他飛行條件不變，目標(biāo)飛行速度超過或飛行高度低于某臨界值時，尾跡中部區(qū)域?qū)π盘柕乃p效應(yīng)會弱于尾跡邊界及頭部區(qū)域。

在未來的研究中，有條件的學(xué)者可利用風(fēng)洞等試驗設(shè)備對高超聲速目標(biāo)繞流流場進(jìn)行參數(shù)驗證。另外在后續(xù)電磁學(xué)仿真中，可增大仿真范圍，并將目標(biāo)外形、來流攻角、信號入射角度、信號的反射和折射納入考慮，為基于GNSS 的高超聲速目標(biāo)無源探測提供更為堅實的依據(jù)。