戰(zhàn)場電磁環(huán)境構(gòu)建中的拋物方程應(yīng)用研究

來嘉哲 方廣有 劉 帥

①(中國科學院電子學研究所 北京 100190)

②(中國科學院研究生院 北京 100049)

③(解放軍裝備學院 北京 101416)

1 引言

在信息化戰(zhàn)爭中，復(fù)雜的戰(zhàn)場電磁環(huán)境已經(jīng)成為影響各類電子裝備作戰(zhàn)效能發(fā)揮的主要要素，戰(zhàn)場電磁環(huán)境的相關(guān)問題研究已成為戰(zhàn)爭仿真的熱點。研究戰(zhàn)場電磁環(huán)境的構(gòu)建方法，模擬出符合實際戰(zhàn)場中復(fù)雜的電磁環(huán)境態(tài)勢，對于推動電子裝備的快速發(fā)展和作戰(zhàn)仿真的應(yīng)用具有重要的意義。

戰(zhàn)場電磁環(huán)境的構(gòu)建可以采用物理或數(shù)學的方法，按照預(yù)先設(shè)置的場景生成電子戰(zhàn)系統(tǒng)所面臨的戰(zhàn)場電磁環(huán)境，模擬戰(zhàn)場上的電子戰(zhàn)行為。戰(zhàn)場電磁環(huán)境的構(gòu)建按具體的實現(xiàn)方式分為以下4種：外場模擬、微波暗室仿真、注入式仿真和全數(shù)字仿真。前3種方式屬于實物半實物仿真的范疇，利用實際的電子設(shè)備形成射頻信號來模擬戰(zhàn)場電磁信號，逼真性好，可用于電子裝備的試驗，但是受試驗環(huán)境的影響和制約，難于體現(xiàn)出實際戰(zhàn)場電磁環(huán)境要素的影響。而全數(shù)字仿真方法可以任意設(shè)置輻射源和環(huán)境要素，手段靈活，利用計算機可快速構(gòu)建任意設(shè)置的戰(zhàn)場電磁環(huán)境，便于系統(tǒng)的分析和設(shè)計。

戰(zhàn)場電磁環(huán)境的數(shù)字化構(gòu)建中，最主要的是能夠計算出各類輻射源在實際傳播要素影響下的傳播過程，在可實用的電波傳播模型中，能夠快速準確地預(yù)測不規(guī)則地形與復(fù)雜大氣條件下電波傳播特性的方法是拋物方程(Parabolic Equation,PE)[1]。與時域有限差分法、射線跟蹤法等方法相比，拋物方程方法計算效率更高更精確，特別適合計算遠距離電波場，得到了人們廣泛關(guān)注。

國外較早開展了基于PE模型的電磁環(huán)境研究，如美國空間和海軍作戰(zhàn)系統(tǒng)中心開發(fā)的 AREPS[2](Advanced Refractive Effects Prediction System)，主要是基于拋物方程方法，可預(yù)測復(fù)雜環(huán)境下電波傳播特性。由于最常用于求解 PE模型的算法是一種步長迭進的方法，難于實現(xiàn)算法并行化，對于戰(zhàn)場電磁環(huán)境的數(shù)字化快速構(gòu)建造成了困難。

本文提出了一種基于并行計算機架構(gòu)的數(shù)字化戰(zhàn)場電磁環(huán)境構(gòu)建方法，在輸入各個輻射源參數(shù)、地形數(shù)據(jù)和大氣參數(shù)的情況下，利用拋物方程模型計算電波傳播的過程，得出整個區(qū)域的場分布，從而模擬出復(fù)雜的戰(zhàn)場電磁環(huán)境。在計算的過程中，設(shè)計了并行算法來充分利用計算集群的運算能力，加速模型的實際求解過程。計算機仿真結(jié)果驗證了該方法具有較好的性能，能滿足大區(qū)域戰(zhàn)場電磁環(huán)境的構(gòu)建需求。文中第2節(jié)對拋物方程模型的基本原理和數(shù)值算法進行了介紹，通過與雙射線模型的比較，驗證了該模型的正確性；第3節(jié)給出了并行算法的基本設(shè)計思路和實現(xiàn)流程；第4節(jié)給出了系統(tǒng)的運行流程和仿真計算結(jié)果，證明了并行算法的有效性；第5節(jié)給出了結(jié)論。

2 拋物方程模型

研究在直角坐標系(x,y,z)中電波的傳播問題，假設(shè)電波只是在2維平面中(x,z)傳播，x為傳播方向，電磁場的時諧因子為e?iωt,i為復(fù)數(shù)單位，ω為角頻率，場量ψ在傳播過程中滿足以下2維標量方程[3]

式(3)中的兩個因式相互獨立，分別代表了電波的前向和后向傳播過程，在只考慮前向傳播的情況下，式(4)已能很好的描述電波的傳播過程：

偽微分算子Q中包含了根式和二階微分形式，在進行求解時無法得出簡單的解析表達式，因此需要對算子作相應(yīng)的近似處理。最常用的方法是采用Taylor近似法，利用Taylor展開的一階近似表達式：

式(6)即為標準拋物方程(Standard Parabolic Equation,SPE)，與式(2)所表示的橢圓方程相比，兩式的主要差別在于標準拋物方程忽略了高階偏微分項，即基于這一事實，忽略了項。SPE的誤差項主要是Taylor展開中忽略的二階無窮小量O(z2)。以傳播仰角為α的平面波為例，其諧函數(shù)表達式通常為u(x,z)=exp[ikx(cosθ?1)+ikzs inθ)],可得出。因此忽略的二次項，即整個方程的誤差正比于 s in4α。在傳播仰角為10°時誤差項正比于10?3,35°時誤差項正比于 1 0?1，可見SPE適合于傳播仰角較小的場合，以目前的非接觸作戰(zhàn)模式為例，輻射源和目標的距離都在上百公里，能夠滿足SPE的傳播仰角需求。

SPE常用的求解方法是分步傅里葉法(Split-Step Fourier Transform,SSFT)[4]。SSFT是一種步進迭代的算法，其基本思路是在當前計算點將高度域變換成頻率域，利用頻率域的快速算法完成計算后變換回高度域，再向前逐步推進，從而計算出整個區(qū)域的場分布。該算法求解速度快，穩(wěn)定性好，本文采用該方法完成戰(zhàn)場電磁環(huán)境的構(gòu)建。

SSFT算法定義傅里葉變換對如下：

對該微分方程求解，可以得出

上式的兩個指數(shù)項分別代表了折射項和繞射項，作傅里葉逆變換后可得出

SSFT算法的步進過程如圖1所示。

圖1 SSFT算法示意圖

實際的戰(zhàn)場上受電磁環(huán)境影響最大的裝備類型是雷達和通信裝備，所以戰(zhàn)場電磁環(huán)境數(shù)字化構(gòu)建過程中考慮的頻率范圍是30 MHz～12.5 GHz。大氣環(huán)境對電波傳播的影響主要體現(xiàn)在式(10)中折射指數(shù)n的變化上，實際環(huán)境中n的值會隨高度和傳播距離發(fā)生變換，通常可以利用通用的折射率模型計算得出，也可依據(jù)實際測量值輸入。戰(zhàn)場環(huán)境的地理覆蓋范圍較大，其地理數(shù)據(jù)可通過 DEM 數(shù)字地圖文件來獲取，并利用數(shù)據(jù)內(nèi)插的方法構(gòu)建出電波傳播的地形剖面。實際的地表起伏不定，可以采用移位變換地形模型來處理地表起伏[5]，并結(jié)合阻抗邊界條件，采用離散混合傅里葉變換算法(Discrete Mixed Fourier Transform,DMFT)來解決邊界條件[6]。

PE方程的推導(dǎo)和 SSFT算法中都采取了一定的近似方法，但其電波傳播預(yù)測的影響程度很小。下文以雙射線模型作為基準參考模型，驗證其正確性。雙射線模型的傳播因子理論值為

式中 ΔR為直達波與反射波的程差，k為自由空間波數(shù)，Γ為地面的歸一化反射系數(shù)，取輻射源天線高度20 m，電磁波頻率為2 GHz，地面介質(zhì)導(dǎo)電率σ=30 mS/m,εr=5，接收點距發(fā)射天線距離為10,000 m，利用SSFT算法的計算結(jié)果和雙射線模型的理論值對比如圖2所示。

圖2 拋物方程與雙射線模型結(jié)果比較

接收點的高度范圍是從地表至300 m，圖2顯示了2種方法計算值的符合程度較高，說明拋物方程模型和SSFT的精度能夠滿足戰(zhàn)場電磁環(huán)境的構(gòu)建需求，所需要解決的問題就是在頻率高、傳播距離遠的情況下如何快速完成傳播過程的計算工作。在單臺工作站計算能力有限的情況下，可以利用高性能計算集群，構(gòu)建適合并行計算的傳播算法，提高算法中最為耗時的離散傅里葉變換(DFT)的計算速度，加速SSFT的求解過程。

3 并行算法設(shè)計

為解決SSFT求解過程中 DFT運算的耗時問題，本文采用了矩陣轉(zhuǎn)置的并行算法，該算法的基本思想是：將1維序列分解成較小規(guī)模的子數(shù)組，構(gòu)成 2維矩陣，繼而對 2維矩陣行、列分別計算DFT[7]。由于矩陣的各行之間、各列之間DFT計算是獨立的，便可以將2維矩陣按行或列分發(fā)給不同的進程，實現(xiàn)數(shù)據(jù)并行計算。該算法可擴展性強，將通信集中于一步內(nèi)完成。

設(shè)N=N1×N2，記

并行計算DFT有如下5步：

(1)每個進程負責矩陣u(n1,n2)的N1P行，并對每行做N2點DFT；

(2)對結(jié)果的每個元素乘以相應(yīng)的旋轉(zhuǎn)因子En1n2=exp(?(2πj/N)n1n2)；

(3)將所得中間結(jié)果矩陣進行轉(zhuǎn)置，使其按列重新分配；

(4)每個進程負責矩陣的N2P列，并對每列做N1點DFT；

(5)最后再次進行轉(zhuǎn)置操作，使矩陣元素在各進程的分布回到原矩陣的順序。

在求解拋物方程的SSFT算法中，需要大量的DFT計算，本文基于矩陣轉(zhuǎn)置的并行計算其核心是：先不做DFT并行算法的第(5)步，而是將繞射項指數(shù)轉(zhuǎn)置，直接與第(4)步所得結(jié)果相乘，這樣避免了第(5)步轉(zhuǎn)置操作所帶來的全局通信開銷，這樣處理提高了并行計算效率。完成繞射項指數(shù)計算后，按照(1)～(4)步相反的順序做IDFT變換，只是其中使用了旋轉(zhuǎn)因子的共軛項，最后計算折射項[8]。

仿真實現(xiàn)中并行算法流程如圖3所示。

圖3 矩陣轉(zhuǎn)置SSFT算法并行方案流程圖

本方案中計算并行性主要體現(xiàn)在如下3方面：

(1)由于折射項、繞射項的計算不依賴于輸入序列，可以在迭代開始前將其計算出來存儲到各節(jié)點的本地存儲器上；

(2)各個進程分別計算出所有的旋轉(zhuǎn)因子，但計算時只取用該進程所對應(yīng)的部分，這樣做同時簡化了程序?qū)崿F(xiàn)難度；

(3)行、列變換在各進程內(nèi)獨立進行。

4 系統(tǒng)仿真結(jié)果

系統(tǒng)采用的集群系統(tǒng)是銀河高性能計算集群，計算單元包含8組×10片刀片式計算機，系統(tǒng)峰值性能為11.2512 TFlops。系統(tǒng)實現(xiàn)過程中使用MPI作為并行計算的消息傳遞接口，它對并行算法的表達能力強，性能高，可移植性好，已經(jīng)成為并行編程事實上的標準。系統(tǒng)的運行流程如圖4所示。

圖4 系統(tǒng)運行流程圖

在完成大氣環(huán)境、地形和輻射源等參數(shù)的輸入之后，系統(tǒng)進行MPI的初始化，將任務(wù)分解到各個運算進程中去，各個進程完成初步運算，并利用MPI的通信機制進行數(shù)據(jù)交換，完成并行算法中的矩陣轉(zhuǎn)置，并逐步向前推進，完成整個區(qū)域的戰(zhàn)場電磁環(huán)境構(gòu)建。

以高度域上的離散點數(shù)為 212個點，推進步數(shù)4000步為例，采用基于矩陣轉(zhuǎn)置的SSFT整體優(yōu)化并行算法對SSFT算法進行求解，其中矩陣的行數(shù)和列數(shù)均取64。在并行機上的測試結(jié)果如表1所示。

表1 并行算法的加速比和效率(212點)

若取高度域上離散點數(shù)為214點，其余的條件不變，在并行機上的測試結(jié)果如表2所示。

表2 并行算法的加速比和效率(214點)

圖 5為在不同的輸入條件下系統(tǒng)的加速比曲線，相比于單機上實現(xiàn)的串行算法，利用并行機和并行算法，速度可以提升到4倍以上。同時仍然可以看出并行算法的加速比相對于理想的線性加速比還是存在著一定的差距，主要的原因在于：在使用矩陣轉(zhuǎn)置的并行方法時，每一步迭代過程都需要進行2次矩陣行和列的轉(zhuǎn)置，用于矩陣中各個元素的DFT和IDFT實現(xiàn)，導(dǎo)致了算法的通信開銷。但是相比于單純的使用FFTW的MPI版本等現(xiàn)成軟件的粗粒度并行方法，其通信開銷已經(jīng)減少了近一半。

圖5 加速比曲線

從圖5還可看出，隨著高度域離散點數(shù)的增加，系統(tǒng)加速比有所增加，這說明并行算法中計算耗費的時間相對于通信時間的比例在增加。綜上所述，隨著用于戰(zhàn)場電磁環(huán)境構(gòu)建任務(wù)的計算節(jié)點的增加，運算速度在不斷的增加，而且隨著采樣點數(shù)的增加，計算強度越來越大，系統(tǒng)的加速比也就越大，因此系統(tǒng)很適合在大區(qū)域戰(zhàn)場電磁環(huán)境的數(shù)字化構(gòu)建。

5 結(jié)論

本文針對利用 PE方程構(gòu)建戰(zhàn)場電磁環(huán)境過程中，對應(yīng)的SSFT算法不易實現(xiàn)并行化的特點，提出了一種基于矩陣轉(zhuǎn)置的并行實現(xiàn)方法。該方法將1維向量構(gòu)造為2維矩陣，并按行或列分發(fā)給并行系統(tǒng)中的各個進程，在各進程中完成繞射項和折射項的計算，降低了并行計算的通信開銷，提高了SSFT算法的并行效率。實驗結(jié)果表明，該方法提高了系統(tǒng)的并行程度和計算速度，具有一定的實際應(yīng)用價值。