陳新蕾，盧立昌，吉洪湖，顧長(zhǎng)青，高帆，施小娟

1.南京航空航天大學(xué) 電子信息工程學(xué)院 雷達(dá)成像與微波光子學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，南京 211106

2.南京航空航天大學(xué) 能源與動(dòng)力學(xué)院，南京 210016

航空發(fā)動(dòng)機(jī)排氣系統(tǒng)是典型的腔體結(jié)構(gòu)，是戰(zhàn)斗機(jī)后向雷達(dá)散射截面（Radar Cross Section， RCS）的主要散射源之一［1-4］。有效降低排氣系統(tǒng)的RCS 是提高作戰(zhàn)飛機(jī)隱身能力和作戰(zhàn)能力的一項(xiàng)重要措施。而在對(duì)發(fā)動(dòng)機(jī)排氣腔體進(jìn)行隱身設(shè)計(jì)的過(guò)程中，對(duì)其電磁散射特性進(jìn)行準(zhǔn)確高效的仿真計(jì)算是一個(gè)重要的環(huán)節(jié)。

發(fā)動(dòng)機(jī)排氣系統(tǒng)內(nèi)部具有復(fù)雜結(jié)構(gòu)，散射物理過(guò)程具有多次反射的特點(diǎn)，是電大散射問(wèn)題，其面臨著準(zhǔn)確求解計(jì)算量大、計(jì)算時(shí)間長(zhǎng)的問(wèn)題。針對(duì)這些問(wèn)題，本文研究排氣系統(tǒng)電磁散射問(wèn)題的準(zhǔn)確快速求解算法。

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

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

1 排氣系統(tǒng)電磁散射問(wèn)題

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

圖1 軸對(duì)稱排氣系統(tǒng)與低散射載體示意圖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 軸對(duì)稱排氣系統(tǒng)的實(shí)物照片F(xiàn)ig.3 Photo of an axisymmetric exhaust system

為了建立求解發(fā)動(dòng)機(jī)排氣系統(tǒng)電磁散射的積分方程，將圖1 所示的排氣系統(tǒng)抽象成圖4。其中，Sin表示腔體的內(nèi)表面；Sout表示低散射載體的外表面，n表示腔體表面的單位外法向。Ei表示外來(lái)平面波的電場(chǎng)，Es表示腔體目標(biāo)的散射電場(chǎng)。

圖4 發(fā)動(dòng)機(jī)排氣系統(tǒng)電磁散射的示意圖Fig.4 Schematic diagram of electromagnetic scattering from an engine exhaust system

根據(jù)等效原理對(duì)圖4 中的腔體散射問(wèn)題建立嚴(yán)格的電磁積分方程，得到電場(chǎng)積分方程（Electric Field Integral Equation，EFIE）和磁場(chǎng)積分方程（Magnetic Field Integral Equation，MFIE）如下：

式中：J(r)表示目標(biāo)表面的感應(yīng)電流；Hi表示外來(lái)平面波的磁場(chǎng)；Hs表示腔體目標(biāo)的散射磁場(chǎng)。根據(jù)場(chǎng)源關(guān)系，Es(r)和Hs(r)可以表示為

其中：G(r，r′)表示自由空間的格林函數(shù)?？梢园l(fā)現(xiàn)，與SBR 和IPO 方法不同，方程（1）和（2）是嚴(yán)格的，沒(méi)有近似，因此可以用于求解不同頻率和不同入射波的電磁散射問(wèn)題。

由于單獨(dú)使用方程（1）或者方程（2）都會(huì)存在偽內(nèi)諧振問(wèn)題，這里使用混合場(chǎng)積分方程（Combined Field Integral Equation，CFIE），它是方程（1）和（2）的線性組合。

求解方程（5）可以得到腔體表面的感應(yīng)電流J(r)，進(jìn)而可以求出腔體的散射場(chǎng)Es(r)及其雷達(dá)散射截面RCS （σ）

2 計(jì)算方法

針對(duì)準(zhǔn)確計(jì)算發(fā)動(dòng)機(jī)排氣系統(tǒng)電磁散射時(shí)間長(zhǎng)的問(wèn)題，本文利用一種并行的插值分解-特征基函數(shù)法-多層快速多極子算法（ID-CBFMMLFMA）對(duì)發(fā)動(dòng)機(jī)排氣系統(tǒng)的單站RCS 進(jìn)行快速準(zhǔn)確仿真。

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

利用矩量法（MoM）［9-10］可以將電磁散射積分方程轉(zhuǎn)化成矩陣方程來(lái)求解。首先，利用RWG 基函數(shù)對(duì)感應(yīng)電流進(jìn)行展開

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

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

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

在CBFM 中，首先要對(duì)目標(biāo)進(jìn)行分塊，然后在每一塊上構(gòu)造出相應(yīng)的特征基函數(shù)。CBF 的生成分為2 步，第1 步生成一組通過(guò)平面波照射得到的初始CBF；第2 步使用奇異值分解去除初始CBF 中的線性相關(guān)部分。通常，保留的CBF數(shù)目要遠(yuǎn)小于RWG 的數(shù)目［13-15］。

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

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

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

圖5 給出了用CBF 壓縮MoM 阻抗矩陣的示意圖，圖中等號(hào)左邊的藍(lán)色小方塊代表傳統(tǒng)MoM 中阻抗矩陣的子矩陣，2 個(gè)細(xì)長(zhǎng)的塊對(duì)角矩陣表示CBF 矩陣，等號(hào)右邊的橙色方塊表示壓縮后的阻抗矩陣的子矩陣?？梢钥闯鯟BFM 可以顯著降低內(nèi)存需求，同時(shí)與之相關(guān)的矩陣向量積（Matrix Vector Product，MVP）的計(jì)算時(shí)間也會(huì)顯著減少。

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

為了進(jìn)一步提高計(jì)算效率和降低內(nèi)存需求，使 用 多 層 快 速 多 極 子 算 法（MLFMA）［11，15］對(duì)CBFM 進(jìn)行加速。要使用MLFMA，需要將目標(biāo)進(jìn)行多層八叉樹分組?；诮M與組之間的相對(duì)位置關(guān)系，將CBFM 的縮減阻抗矩陣ZR分解為和兩部分

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

與傳統(tǒng)的MLFMA 相比，CBFM-MLFMA在MVP 計(jì)算時(shí)間和內(nèi)存消耗方面也更具優(yōu)勢(shì)。如果CBFM 的基函數(shù)壓縮率為，也就是說(shuō)基函數(shù)數(shù)目變?yōu)樵瓉?lái)的1 R。那么，在CBFM-MLFMA 中，聚集和發(fā)散函數(shù)耗費(fèi)的內(nèi)存和迭代時(shí)間會(huì)變?yōu)閭鹘y(tǒng)MLFMA 的1 R，近場(chǎng)阻抗矩陣耗費(fèi)的內(nèi)存和迭代時(shí)間會(huì)變成傳統(tǒng)MLFMA 的1 R2。例如壓縮率R=10，那么聚集和發(fā)散函數(shù)耗費(fèi)的內(nèi)存和MVP 時(shí)間會(huì)變?yōu)樵瓉?lái)的10%，耗費(fèi)的內(nèi)存和MVP 時(shí)間會(huì)變成原來(lái)的1%，效果很可觀。

2. 2 縮減激勵(lì)矩陣的低秩壓縮

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

式中：Q 是酉矩陣；R 是上三角矩陣。利用誤差門限ε 對(duì)R 進(jìn)行截?cái)啵?dāng)R 中的第p 行的對(duì)角線元素滿足以下條件

將R 中相應(yīng)的行去掉。這里，rp表示R 對(duì)角線上的第p 個(gè)值。ε 通常取值為10-3～10-4。

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

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

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

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

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

2. 3 并行加速計(jì)算

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

在特征基函數(shù)的生成步驟，需要對(duì)每組進(jìn)行CBF 構(gòu)造。由于每組CBF 的構(gòu)造互不相關(guān)，因此非常適合并行。按組對(duì)CBF 構(gòu)造進(jìn)行并行，每個(gè)線程不存在競(jìng)爭(zhēng)，可以獲得很高的并行效率。與CBF 構(gòu)造步驟相類似，聚集和發(fā)散函數(shù)的生成以及近場(chǎng)縮減矩陣的生成也是組與組之間互不相關(guān)的，非常適合并行，這里仍然是按組進(jìn)行并行。在MVP 步驟，由于MLFMA 算法結(jié)構(gòu)復(fù)雜，并行起來(lái)較為復(fù)雜。這里使用按組和按平面波2 種并行方式相結(jié)合的方式進(jìn)行并行實(shí)現(xiàn)。

3 仿真結(jié)果與試驗(yàn)驗(yàn)證

試驗(yàn)件的結(jié)構(gòu)圖如圖1 所示，而實(shí)物圖如圖3所示。該試驗(yàn)件為純金屬結(jié)構(gòu)，放置在微波暗室中。RCS 測(cè)試系統(tǒng)的示意如圖7 所示，包含天線、泡沫轉(zhuǎn)臺(tái)等。測(cè)量時(shí)方位角取值從-180°～180°，采樣間隔為0.2°，共1 801 個(gè)角度。

圖7 RCS 測(cè)試系統(tǒng)示意圖Fig.7 Schematic diagram of an RCS test system

利用并行ID-CBFM-MLFMA 對(duì)軸對(duì)稱排氣系統(tǒng)的單站RCS 進(jìn)行仿真計(jì)算，并與試驗(yàn)測(cè)試結(jié)果進(jìn)行比較以驗(yàn)證算法精度。仿真模型的結(jié)構(gòu)如圖1 所示，是帶有低散射載體的排氣系統(tǒng)。網(wǎng)格剖分尺寸為0.1 波長(zhǎng)。本文計(jì)算的單站RCS 均為VV 極化。仿真時(shí)方位角取值從-180°～180°，采樣間隔為1°，共361 個(gè)角度。

3. 1 計(jì)算精度

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

圖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 對(duì)比了35 GHz 頻率下載體以及排氣系統(tǒng)（含載體）的單站RCS?？梢钥闯觯?°方位角附近角域，載體的單站RCS 遠(yuǎn)小于整體目標(biāo)的單站RCS。因此，在這些角域可以認(rèn)為整體目標(biāo)的單站RCS 主要是由排氣腔體貢獻(xiàn)的。

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

3. 2 計(jì)算效率

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

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

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

4 結(jié) 論

本文將并行的ID-CBFM-MLFMA 算法應(yīng)用于航空發(fā)動(dòng)機(jī)排氣腔體的電磁散射特性的仿真計(jì)算，并與試驗(yàn)測(cè)試結(jié)果進(jìn)行對(duì)比驗(yàn)證。

1）與試驗(yàn)測(cè)試結(jié)果對(duì)比表明ID-CBFMMLFMA 具有良好的計(jì)算精度。和傳統(tǒng)CBFMMLFMA 方法的精度幾乎相同。

2）在計(jì)算航空發(fā)動(dòng)機(jī)排氣系統(tǒng)模型多角度單站RCS 時(shí)，與傳統(tǒng)的CBFM-MLFMA 相比，ID-CBFM-MLFMA 具有更高的計(jì)算效率。在本文的算例中，迭代求解時(shí)間減少了78%。

3）理論上，并行ID-CBFM-MLFMA 是一種通用的算法，也可用于飛行器整機(jī)、發(fā)動(dòng)機(jī)進(jìn)氣道等各種目標(biāo)電磁散射特性的快速準(zhǔn)確仿真。