航空發動機排氣系統單站RCS 準確高效仿真方法與試驗驗證

陳新蕾，盧立昌，吉洪湖，顧長青，高帆，施小娟

1.南京航空航天大學 電子信息工程學院 雷達成像與微波光子學教育部重點實驗室，南京 211106

2.南京航空航天大學 能源與動力學院，南京 210016

航空發動機排氣系統是典型的腔體結構，是戰斗機后向雷達散射截面（Radar Cross Section， RCS）的主要散射源之一［1-4］。有效降低排氣系統的RCS 是提高作戰飛機隱身能力和作戰能力的一項重要措施。而在對發動機排氣腔體進行隱身設計的過程中，對其電磁散射特性進行準確高效的仿真計算是一個重要的環節。

發動機排氣系統內部具有復雜結構，散射物理過程具有多次反射的特點，是電大散射問題，其面臨著準確求解計算量大、計算時間長的問題。針對這些問題，本文研究排氣系統電磁散射問題的準確快速求解算法。

目前用于排氣系統電磁散射計算方法主要有彈跳射線法（Shooting and Bouncing Rays，SBR）［5］、迭代物理光學法（Iterative Physical Optics，IPO）［6-8］、矩 量 法（Method of Moments，MoM）［9-10］等。SBR 和IPO 屬于高頻漸近方法，計算速度、快占用內存少，但是其計算精度低且不好控制。MoM 是一種基于頻域電磁積分方程的準確數值方法，將積分方程離散成矩陣方程求解。對于內部具有復雜結構的發動機排氣腔體，與高頻漸近方法相比，MoM 可以給出更準確的結果［11］。但是傳統的MoM 迭代求解方法的計算和存儲復雜度很高，導致其難以計算電大排氣腔體。多層快速多極子算法（Multilevel Fast Multipole Algorithm，MLFMA）［12］是MoM 的一種快速求解算法，可以顯著降低MoM 的計算和存儲復雜度，但是該算法并不改變基函數的數目，在計算電大腔體時仍然花費較多的計算時間和內存。特征基函數法（Characteristic Basis Function Method，CBFM）［13-15］可以顯著降低基函數的數目，與MLFMA 相結合可以獲得更高的計算效率并節省更多的計算內存［16-17］。

本文利用插值分解（Interpolative Decomposition，ID）［18-19］進 一 步 改 進CBFM-MLFMA 的計算效率，并將改進后的方法應用于航空發動機排氣系統的電磁散射特性仿真中。插值分解可以將CBFM-MLFMA 中單站激勵矩陣壓縮分解成2 個小矩陣的乘積，因此CBFM-MLFMA 只需求解較少數的平面波激勵，而其他平面波激勵的解可以通過已求出解的插值得到，而不需要再進行耗時的迭代求解。 本文利用ID-CBFMMLFMA 對發動機排氣系統的單站RCS 進行仿真計算，顯著提高了計算效率，并與試驗測試結果進行對比，驗證了算法的精度。

1 排氣系統電磁散射問題

研究發動機排氣系統腔體電磁散射時，為了排除發動機排氣腔體外表面散射的影響，常將發動機排氣腔體安裝在一個低散射載體中。圖1 給出了一個典型的軸對稱排氣系統安裝在一個低散射載體（見圖2）中形成的電磁散射研究模型示意圖，主要包括腔體內壁面、中心錐、火焰穩定器、支板、低散射載體等組成。渦輪葉片近似處理成金屬平板。圖3 為研究模型實物照片。

圖1 軸對稱排氣系統與低散射載體示意圖Fig.1 Schematic diagram of an axisymmetric exhaust system with a low scattering carrier

圖2 低散射載體模型示意圖Fig.2 Schematic diagram of a low scattering carrier

圖3 軸對稱排氣系統的實物照片Fig.3 Photo of an axisymmetric exhaust system

為了建立求解發動機排氣系統電磁散射的積分方程，將圖1 所示的排氣系統抽象成圖4。其中，Sin表示腔體的內表面；Sout表示低散射載體的外表面，n表示腔體表面的單位外法向。Ei表示外來平面波的電場，Es表示腔體目標的散射電場。

圖4 發動機排氣系統電磁散射的示意圖Fig.4 Schematic diagram of electromagnetic scattering from an engine exhaust system

根據等效原理對圖4 中的腔體散射問題建立嚴格的電磁積分方程，得到電場積分方程（Electric Field Integral Equation，EFIE）和磁場積分方程（Magnetic Field Integral Equation，MFIE）如下：

式中：J(r)表示目標表面的感應電流；Hi表示外來平面波的磁場；Hs表示腔體目標的散射磁場。根據場源關系，Es(r)和Hs(r)可以表示為

其中：G(r，r′)表示自由空間的格林函數。可以發現，與SBR 和IPO 方法不同，方程（1）和（2）是嚴格的，沒有近似，因此可以用于求解不同頻率和不同入射波的電磁散射問題。

由于單獨使用方程（1）或者方程（2）都會存在偽內諧振問題，這里使用混合場積分方程（Combined Field Integral Equation，CFIE），它是方程（1）和（2）的線性組合。

求解方程（5）可以得到腔體表面的感應電流J(r)，進而可以求出腔體的散射場Es(r)及其雷達散射截面RCS （σ）

2 計算方法

針對準確計算發動機排氣系統電磁散射時間長的問題，本文利用一種并行的插值分解-特征基函數法-多層快速多極子算法（ID-CBFMMLFMA）對發動機排氣系統的單站RCS 進行快速準確仿真。

2. 1 特征基函數法-多層快速多極子算法

利用矩量法（MoM）［9-10］可以將電磁散射積分方程轉化成矩陣方程來求解。首先，利用RWG 基函數對感應電流進行展開

式 中：fn(r) 表 示 第n 個Rao-Wilton-Glisson（RWG）基函數［20］；In為待求的展開系數。然后，將式（7）代入CFIE，并利用Galerkin 檢測，得到最終需求解的矩陣方程

式中：Z 和V 分別稱為廣義阻抗矩陣和廣義電壓矩陣。

CBFM 是一種宏基函數法，通過在較大區域上定義宏基函數，稱之為特征基函數（CBF），可以實現未知量的縮減，從而降低阻抗矩陣所占內存，并減少阻抗矩陣方程的迭代求解時間。并且所構造的特征基函數具有激勵無關（Excitation Independent，EI）的 特 性［13］，其 主 要 特 點是：CBF 只需要構造一次，就可以應用于不同極化和方向的平面波照射的散射問題，而不需要重復構造CBF，這種特性非常適合于單站RCS的計算。

在CBFM 中，首先要對目標進行分塊，然后在每一塊上構造出相應的特征基函數。CBF 的生成分為2 步，第1 步生成一組通過平面波照射得到的初始CBF；第2 步使用奇異值分解去除初始CBF 中的線性相關部分。通常，保留的CBF數目要遠小于RWG 的數目［13-15］。

當所有塊上的CBF 都生成好之后，可以用CBF 壓縮傳統MoM 的阻抗矩陣，得到縮減后的阻抗矩陣方程

式中：ZR表示利用CBF 壓縮后的矩陣，稱之為縮減 矩陣；VR是CBFM 中的電壓 向量；α 為 待求的CBF 的線性表出系數。

式中：J 由各塊保留的CBF 聚集而成。

圖5 給出了用CBF 壓縮MoM 阻抗矩陣的示意圖，圖中等號左邊的藍色小方塊代表傳統MoM 中阻抗矩陣的子矩陣，2 個細長的塊對角矩陣表示CBF 矩陣，等號右邊的橙色方塊表示壓縮后的阻抗矩陣的子矩陣。可以看出CBFM 可以顯著降低內存需求，同時與之相關的矩陣向量積（Matrix Vector Product，MVP）的計算時間也會顯著減少。

圖5 CBFM 對MoM 阻抗矩陣的壓縮Fig.5 Compression of MoM impedance matrix by using CBFM

為了進一步提高計算效率和降低內存需求，使 用 多 層 快 速 多 極 子 算 法（MLFMA）［11，15］對CBFM 進行加速。要使用MLFMA，需要將目標進行多層八叉樹分組。基于組與組之間的相對位置關系，將CBFM 的縮減阻抗矩陣ZR分解為和兩部分

利用加法定理，可以將CBFM 中與遠場阻抗相關的MVP 轉化成“聚集-轉移-發散”的形式［16-17］，這樣可以大大提高MVP 的計算效率。在此基礎上，再利用多層分組、插值、反插值技術，MLFMA 可以將CBFM 的MVP 的計算和存儲復雜度降低至O(NCBFlg NCBF)，這里NCBF為CBF 基函數的總數目。

與傳統的MLFMA 相比，CBFM-MLFMA在MVP 計算時間和內存消耗方面也更具優勢。如果CBFM 的基函數壓縮率為，也就是說基函數數目變為原來的1 R。那么，在CBFM-MLFMA 中，聚集和發散函數耗費的內存和迭代時間會變為傳統MLFMA 的1 R，近場阻抗矩陣耗費的內存和迭代時間會變成傳統MLFMA 的1 R2。例如壓縮率R=10，那么聚集和發散函數耗費的內存和MVP 時間會變為原來的10%，耗費的內存和MVP 時間會變成原來的1%，效果很可觀。

2. 2 縮減激勵矩陣的低秩壓縮

CBFM 降低了基函數的數目，MLFMA 降低了迭代求解中MVP 的計算復雜度，但是當計算單站RCS 時，由于需要對多入射角進行迭代求解，而迭代求解時間正比于入射角個數，因此這導致計算時間仍然可能很長。當入射平面波采樣過密的時候，激勵矩陣往往存在低秩 特 性，可 以 利 用 低 秩 分 解 算 法［15，19，21］進 行 壓縮 分 解［18，22-23］。這 里 利 用 基 于 選 主 元 的QR 分解對CBFM-MLFMA 中的激勵矩陣進行壓縮分解。

式中：Q 是酉矩陣；R 是上三角矩陣。利用誤差門限ε 對R 進行截斷，當R 中的第p 行的對角線元素滿足以下條件

將R 中相應的行去掉。這里，rp表示R 對角線上的第p 個值。ε 通常取值為10-3～10-4。

假設R 矩陣最終只有k 個行被保留，那么VR可以近似表示為

式中：R11是上三角矩陣；；Q1是Q 的一個子矩陣；是VR的一個子矩陣。如圖6 所示，就是VR的插值分解［16-17］。

圖6 CBFM-MLFMA中單站激勵矩陣的插值分解示意圖Fig.6 Schematic diagram of interpolative decomposition of monostatic excitation matrix in CBFMMLFMA

因此，只需要利用CBFM-MLFMA 求解下面的矩陣方程

則最終需要的解可以表示為可以看出，利用插值分解只需要求解k 次矩陣方程，因此計算效率得到提高。

2. 3 并行加速計算

為了充分利用CPU 的計算資源，利用共享內存的OpenMP 技術將ID-CBFM-MLFMA 程序進行并行。將ID-CBFM-MLFMA 算法中的4 個主要步驟進行并行：① 特征基函數的生成；② 聚集和發散函數的生成；③ 近場縮減矩陣的生成；④ 迭代求解時的MVP。

在特征基函數的生成步驟，需要對每組進行CBF 構造。由于每組CBF 的構造互不相關，因此非常適合并行。按組對CBF 構造進行并行，每個線程不存在競爭，可以獲得很高的并行效率。與CBF 構造步驟相類似，聚集和發散函數的生成以及近場縮減矩陣的生成也是組與組之間互不相關的，非常適合并行，這里仍然是按組進行并行。在MVP 步驟，由于MLFMA 算法結構復雜，并行起來較為復雜。這里使用按組和按平面波2 種并行方式相結合的方式進行并行實現。

3 仿真結果與試驗驗證

試驗件的結構圖如圖1 所示，而實物圖如圖3所示。該試驗件為純金屬結構，放置在微波暗室中。RCS 測試系統的示意如圖7 所示，包含天線、泡沫轉臺等。測量時方位角取值從-180°～180°，采樣間隔為0.2°，共1 801 個角度。

圖7 RCS 測試系統示意圖Fig.7 Schematic diagram of an RCS test system

利用并行ID-CBFM-MLFMA 對軸對稱排氣系統的單站RCS 進行仿真計算，并與試驗測試結果進行比較以驗證算法精度。仿真模型的結構如圖1 所示，是帶有低散射載體的排氣系統。網格剖分尺寸為0.1 波長。本文計算的單站RCS 均為VV 極化。仿真時方位角取值從-180°～180°，采樣間隔為1°，共361 個角度。

3. 1 計算精度

首先，考慮軸對稱排氣系統在0°俯仰角的時候，在多個頻率下的單站RCS。圖8～圖10 分別給 出 了 ID-CBFM-MLFMA、傳 統 CBFMMLFMA 對軸對稱排氣腔體試驗件在6 GHz、10 GHz、35 GHz 頻率下，單站RCS 隨方位角的變化進行仿真計算。 可以看出ID-CBFMMLFMA 的仿真結果與傳統CBFM-MLFMA 的仿真結果幾乎完全重合，ID-CBFM-MLFMA 的仿真結果和試驗測試結果吻合良好。

圖8 6 GHz 頻率下VV 極化單站RCS 隨方位角的變化Fig.8 VV-polarized monostatic RCS versus azimuth angle at 6 GHz

圖9 10 GHz 頻率下VV 極化單站RCS 隨方位角的變化Fig.9 VV-polarized monostatic RCS versus azimuth angle at 10 GHz

圖10 35 GHz頻率下VV 極化單站RCS 隨方位角的變化Fig.10 VV-polarized monostatic RCS versus azimuth angle at 35 GHz

圖11 對比了35 GHz 頻率下載體以及排氣系統（含載體）的單站RCS。可以看出，在0°方位角附近角域，載體的單站RCS 遠小于整體目標的單站RCS。因此，在這些角域可以認為整體目標的單站RCS 主要是由排氣腔體貢獻的。

圖11 35 GHz頻率下VV極化載體與排氣系統的單站RCS對比Fig.11 VV-polarized monostatic RCS comparison of exhaust system and carrier at 35 GHz

3. 2 計算效率

接下來將ID-CBFM-MLFMA 與傳統的CBFM-MLFMA 進行比較來說明ID-CBFMMLFMA 的計算效率。所使用的計算機配置是Intel （R） Core （TM） i9-10900K CPU，主 頻3.7 GHz，10 核20 線程，總內存128 GB。并行線程數均設置為20。GMRES 收斂迭代收斂精度設置為0.001。使用了塊對角預處理器。計算的是軸對稱排氣系統10 GHz 頻率下VV 極化單站RCS。仿真時方位角取值從-180°～180°，采樣間隔為1°，共361 個角度。插值分解的截斷門限設置為0.001。

使用插值分解對激勵矩陣進行壓縮，得到有效秩為75。這意味著借助于插值分解，可以將迭代求解矩陣方程的次數減少到75 次，而傳統的CBFM-MLMFA 則需要迭代求解矩陣方程361 次，因此ID-CBFM-MLFMA 可以顯著降低迭代求解的時間，這對于計算發動機排氣腔體的寬角RCS 特性很有幫助。表1 給出了2 種方法計算時間的比較，可以看出，與CBFMMLFMA 相 比，ID-CBFM-MLFMA 在 計 算CBF、近場縮減矩陣、聚集轉移發散函數時耗費的時間相同，但是顯著降低了迭代求解時的時間。

表1 計算10 GHz頻率下VV 極化RCS的計算時間的比較Table 1 Comparison of calculation time for VVpolarized RCS at 10 GHz

4 結 論

本文將并行的ID-CBFM-MLFMA 算法應用于航空發動機排氣腔體的電磁散射特性的仿真計算，并與試驗測試結果進行對比驗證。

1）與試驗測試結果對比表明ID-CBFMMLFMA 具有良好的計算精度。和傳統CBFMMLFMA 方法的精度幾乎相同。

2）在計算航空發動機排氣系統模型多角度單站RCS 時，與傳統的CBFM-MLFMA 相比，ID-CBFM-MLFMA 具有更高的計算效率。在本文的算例中，迭代求解時間減少了78%。

3）理論上，并行ID-CBFM-MLFMA 是一種通用的算法，也可用于飛行器整機、發動機進氣道等各種目標電磁散射特性的快速準確仿真。