基于離散伴隨方程的三維雷達散射截面幾何敏感度計算

周琳，黃江濤，高正紅,*

1. 西北工業大學 航空學院，西安 710072 2. 中國空氣動力研究與發展中心，綿陽 621000

雷達散射截面(Radar Cross Section, RCS)反映了物體在給定方向上對入射雷達波散射的強弱，是衡量飛機隱身性能的重要指標,考慮隱身的飛行器設計常以減小雷達散射截面作為隱身設計的重要目標。現有飛行器隱身設計中多采用幾何光學法(GO)、物理光學法(PO)、幾何繞射理論(GTD)、物理繞射理論(PDT)等高頻近似算法評估散射體的雷達散射截面[1-2]。高頻算法根據高頻場的局部性原理，僅根據入射場獨立地近似確定表面感應電流[3]，計算速度快，所需內存小。近年來，隨著反隱身技術的發展，尤其是全數字化有源相控技術在低空警戒和監視雷達中的應用，對飛行器隱身性能的要求越來越高[4]。高頻方法雖然可以預估飛行器目標的鏡面散射、邊緣繞射等散射特性，但由于算法本身的局限，忽略了一些關鍵部件間的耦合關系，在弱散射源計算中的精度不足，使其在低RCS構型的評估和優化過程中可信度較低，無法滿足低RCS構型精細化設計的需求[5-6]。基于嚴格理論模型的數值解法如矩量法(MOM)及基于矩量法的快速解法(如多層快速多極子算法(MLFMA))[7-8]從電磁場積分方程出發，可以精確求解三維復雜外形目標的電磁散射。隨著高性能計算技術和低可探測飛行器的發展，基于精確建模的積分方程數值解法逐漸成為飛行器隱身設計中重要的電磁分析手段。

飛行器隱身性能與其外形密切相關，設計中常需處理隱身指標與氣動指標之間的矛盾[9]。現有的氣動隱身一體化設計多采用粒子群算法、遺傳算法、神經網絡算法等智能搜索算法[1, 6, 9-10]。智能搜索算法具有收斂到全局最優的能力，但優化效率較低，優化所需的雷達散射截面評估次數隨設計問題規模的增加迅速增加，與高精度電磁模擬方法結合時計算成本昂貴。基于梯度的優化算法效率較高[11]，但經典基于有限差分的RCS梯度計算代價大，限制了梯度算法在氣動/隱身優化中的應用，目前對基于梯度的氣動/隱身一體化優化研究較少。基于伴隨思想的梯度求解方法可以通過一次原方程求解和一次伴隨方程求解得到目標關于所有設計變量的導數，目前已在氣動外形優化[12]、氣動/結構耦合優化[13]、聲爆優化[14]等領域得到廣泛應用。基于伴隨方法的優化在隱身領域研究較少，Bondeson等[15]推導了二維積分形式的麥克斯韋連續伴隨方程，并對二維簡單外形(圓形)的RCS進行了優化；Wang和Anderson[16]推導了二維時域麥克斯韋方程的伴隨方程，通過給定翼型的目標電流分布，對翼型的幾何進行反設計；本文作者[17]在三維矩量法的基礎上推導了麥克斯韋方程的離散伴隨方程，驗證了伴隨方法的精度，但仍面臨矩量法求解計算成本高、工程應用困難的問題。

針對三維問題中高精度RCS評估梯度計算成本高的問題，本文引入麥克斯韋積分方程的伴隨方程，以三維矩量法為例詳細推導了麥克斯韋積分方程的離散伴隨方程和雷達散射截面的變分表達式。針對矩量法在三維復雜問題中內存、計算量大的問題，采用多層快速多極子算法(MLFMA)進行求解。選取雙椎體模型、導彈模型對求解器的計算精度、伴隨方法的可靠性進行驗證。計算結果表明本文采用的求解器有較高精度，基于矩量法、多層快速多極子算法伴隨梯度與差分梯度吻合良好，實現了通過兩次線性方程組求解獲得雷達散射截面關于所有設計變量的導數。證明了伴隨方法的可靠性和高效性，為基于梯度的高精度氣動/隱身一體化優化奠定基礎。

1 積分方程數值求解方法

本文采用矩量法和多層快速多極子算法求解麥克斯韋積分方程及其離散伴隨方程。矩量法是文獻[7]提出的一種用于嚴格計算電磁問題的數值方法。矩量法將積分形式的麥克斯韋方程離散后轉化為矩陣方程，通過求解線性方程組來得到目標表面感應電流，進而計算散射場。快速多極子算法(FMM)[8]和多層快速多極子算法[18]是建立在矩量法和線性方程組迭代求解算法基礎上的快速算法。

1.1 基于RWG基函數的矩量法

RWG(Rao-Wilton-Glisson)基函數[19]是當前使用較為廣泛的一種矩量法基函數，該基函數定義在具有公共邊的相鄰三角形上，可模擬任意形狀物體的表面電、磁流分布。對于RWG基函數矩量法，三角形剖分的網格邊長一般在[λ/12，λ/8]之間較為合適[20](λ為入射波長)。根據麥克斯韋方程和導體表面S上切向電場連續條件可得電場積分方程為

(1)

(2)

(3)

(4)

圖1 RWG基函數示意圖[21]Fig.1 Sketch of RWG basis function[21]

(5)

式中：In為第n個基函數的系數；N為三角形公共邊總數。采用RWG基函數作為檢驗函數(伽略金法)，檢測電場積分方程可得

(6)

整理寫成矩陣形式：

ZN×NIN×1=VN×1

(7)

式中：Z為矩量法阻抗矩陣；I和V分別為電流系數向量和激勵向量。阻抗元素和激勵項的表達式為[7, 22]

(8)

(9)

1.2 多層快速多極子算法

矩量法離散得到的阻抗矩陣是復數稠密矩陣，如果用直接法求解矩陣方程，則算法存儲量為O(e2)，計算量為O(e3)(相對于矩陣求逆運算量級)；如果用迭代法求解，存儲量為O(e2)，每步迭代的計算量為O(e2)(相對于矩陣矢量乘法運算量級)，其中e為未知量的個數。當頻率增加時，未知量個數迅速增加，對內存和計算速度提出了很高要求，限制了矩量法在高頻散射問題中的應用。為提高矩量法求解問題的規模，在矩量法和線性方程組迭代求解算法的基礎上發展了快速多極子算法和多層快速多極子算法，將計算量和存儲量進一步降低到O(e1.5)和O(elge)量級，使基于精確建模的積分方程數值解法求解電大尺寸目標的電磁散射問題成為可能[23]。

快速多極子算法把基函數分為G組，將組之間的作用分為近相互作用和遠相互作用，即將阻抗矩陣表示為

Z=Znear+Zfar

(10)

對于近相互作用，仍采用矩量法計算；對于遠相互作用，則采用快速多極子算法。快速多極子算法利用加法定理和平面波展開定理將格林函數展開為多極子形式，進而將矩陣與向量的乘積運算轉化為聚合-轉移-配置3個過程。標量格林函數G(r′,r)=e-jkR/R的多極子展開形式為

(11)

(12)

(13)

式中：

(14)

(15)

其中：ji和tj分別為電流J(r)的基函數和測試函數;右端第1項表示鄰近組的貢獻；右端第2項為各遠距離組的貢獻，通過式(14)、式(12)和式(15) 表示的聚合、轉移和配置3個步驟得到。

多層快速多極子算法是快速多極子算法的多層擴展。多層快速多極子算法基于樹形結構，通過在多個層級上分組，組間嵌套，逐層遞推來實現FMM[23]，進一步降低了內存需求，提高了矢量矩陣相乘的計算效率。由于FMM和MLFMA本質上加速的是線性方程求解過程中矩陣和向量的乘積運算，在計算過程中并不顯式存儲遠相互作用矩陣。

2 基于伴隨方程的梯度計算

隱身優化設計中常需要雷達散射截面關于設計變量的梯度信息。傳統的有限差分法計算梯度所需雷達散射截面的評估次數與設計變量個數成正比，當設計變量較多時，有限差分法的計算代價較高，難以在實際問題中應用。基于伴隨方程的梯度計算可以通過一次雷達散射截面求解、一次伴隨方程求解獲得雷達散射截面關于所有設計變量的梯度，有效提高梯度計算效率。本節在1.1節介紹的矩量法基礎上對離散伴隨方程和雷達散射截面的變分形式進行推導[17]，并對伴隨方法所需的計算量和存儲量進行簡要分析。

2.1 伴隨方程推導

典型的雷達散射截面優化設計問題形式為

(16)

需要滿足的約束為

R(I(D),D)=V-ZI=0

(17)

式中：σ為雷達散射截面；I為感應電流；D為設計變量。雷達散射截面是幾何外形和感應電流的函數，在矩量法求解中阻抗矩陣、激勵和解得的感應電流均只由幾何形狀決定，即

(18)

引入拉格朗日乘子Λ可以構造目標函數：

L=σ+ΛTR

(19)

對式(19)進行求導，有

(20)

從式(20)可以看出，若找到合適的Λ使式(20) 右端第1項為0，可完全消除感應電流對設計變量導數dI/dD的計算，大幅度降低計算量，即

(21)

注意到：

(22)

將式(22)代入式(21)即可得到伴隨方程：

(23)

伴隨方程的形式與正計算(RCS評估)的形式一致，可以直接采用正計算采用的數值方法(矩量法、多層快速多極子算法)求解。將式(23)代入式(20)，基于伴隨方法的目標梯度可以寫為

(24)

(25)

2.2 雷達散射截面變分推導

伴隨方程式的右端項激勵?σ/?I為雷達散射截面關于感應電流的導數，以1.1節中介紹的基于RWG基函數的三維矩量法為例對伴隨方程的激勵項進行推導。

通過矩量法或多層快速多極子算法解得表面感應電流分布后，則可以計算雷達散射截面標量(m2)[24]，即

(26)

根據復數的運算性質將模的平方寫成自身與其共軛相乘的形式：

(27)

(28)

式中:上劃線表示共軛。由鏈式求導法則：

(29)

復數求導需分別對復數的實部和虛部求導，即

(30)

則式(29)可以整理為

(31)

式(28)中：g離散后可以寫成感應電流和的形式，當采用1.1節介紹的RWG離散感應電流時，具體表達式為

(32)

(33)

將式(33)代入式(31)，有

(34a)

(34b)

將式(34)進一步代入式(30)即可得到伴隨方程右端項的表達式為

(35)

2.3 伴隨方法計算量和存儲量分析

基于伴隨方法的雷達散射截面導數求解主要由RCS正計算、伴隨方程求解、導數計算3部分構成。

在RCS正計算中，矩量法需要一次阻抗矩陣填充和一次線性方程組求解；多層快速多極子算法需要多次計算ZI，計算次數與所需迭代次數成正比。當采用迭代法求解時，2種算法的計算量和存儲量分別為O(e2)和O(elge)。在伴隨方程求解中，若正計算的阻抗矩陣為對稱陣，則矩量法可以直接使用正計算中的阻抗矩陣，不需要重新填充伴隨方程的阻抗矩陣，可通過一次線性方程組求解完成伴隨變量計算；多層快速多極子算法不顯式存儲阻抗矩陣，因此多層快速多極子算法在伴隨計算階段與正計算所需的計算量基本一致。

3 梯度驗證

梯度驗證采用基于電場積分方程的矩量法和快速多極子算法，線性方程組求解采用基于雙共軛梯度算法[25]的迭代法求解。雙共軛梯度算法較之普通的共軛梯度方法有較快收斂性，特別是在矩陣為對稱陣的情況下，可以有效降低方程組求解的計算量[23]，計算過程中取收斂閾值ε=10-5。本節采用計算電磁學經典算例及實際工程中復雜外形對伴隨梯度的可靠性和精度進行探討。

3.1 雙椎體模型

雙椎體模型長7.5 in(1 in=25.4 mm)，由2個 椎角分別為22.62°和46.4°的半椎拼接而成[26-27]，是EMCC(ElectroMagnetic Code Consortium)[27]提供用于校驗計算電磁學代碼的標準算例之一。雙椎體模型兩端存在尖點，高頻方法難以準確計算，需采用高精度數值求解方法。計算采用的雷達波頻率f=9 GHz，單站散射水平極化和垂直極化，網格平均尺寸約為λ/10，未知量總數e=11 094。采用矩量法、多層快速多極子算法計算，并用商用軟件和實驗值(EXP)[26-27]校核計算結果。幾種算法的計算結果如圖2和圖3所示。本文采用的矩量法和多層快速多極子算法均與實驗值及商用軟件FEKO矩量法的計算結果吻合良好，具有較高精度。

采用基于HMQ全局徑向基函數的域元法[28-29]對雙椎體模型進行參數化，設計變量分布如圖4所示。采用基于矩量法、多層快速多極子算法的伴隨方法和基于矩量法的有限差分法求解各設計變量在Z方向的梯度，認為基于矩量法的有限差分(FDD)結果為梯度的精確解。有限差分計算中擾動量Δx=10-3m。設計變量A和B(見圖4)在各入射角度的導數計算結果如圖5和圖6所示。

圖2 雙椎體水平極化計算結果Fig.2 Numerical results of double-ogive horizontal polarization

基于伴隨的矩量法和快速多極子算法計算得到的梯度與有限差分得到的梯度在趨勢上吻合較好。其中，基于伴隨的矩量法與有限差分的計算結果基本一致，具有較高精度；多層快速多極子算法得到的梯度與有限差分在梯度較小的角域存在一定誤差，在梯度峰值附近吻合良好。入射角為0°時的感應電流模值分布和伴隨電流模值分布如圖7所示。值得注意的是，在流場伴隨求解中，流場伴隨變量傳播方向與流場守恒變量相反[12]，而在電磁場求解中，電磁場伴隨變量分布與感應電流分布較為一致。在求解過程中，共計算180個入射角度，其中矩量法計算使用20核并行計算，最大內存占用1.85 G，總計算耗時42 min。多層快速多極子算法計算使用8核并行計算，最大內存使用18.3 MB，計算用時3.7 min。

圖3 雙椎體垂直極化計算結果Fig.3 Numerical results of double-ogive vertical polarization

圖4 雙椎體設計變量分布Fig.4 Design variables distribution of double-ogive

圖5 雙椎體設計變量A導數(水平極化)Fig.5 Gradients of double-ogive design variable A(horizontal polarization)

圖6 雙椎體設計變量B導數(水平極化)Fig.6 Gradients of double-ogive design variable B(horizontal polarization)

圖7 雙椎體感應電流及伴隨變量云圖Fig.7 Contours of surface current and adjoint variable of double-ogive

3.2 導彈模型

導彈模型是宏觀上的電大尺寸和細節上的電小尺寸共存的復雜模型，具有較高的實際工程意義，對電磁散射計算和梯度計算提出了較高要求。計算采用的導彈模型全長4 m，彈體直徑0.28 m。計算采用雷達波段頻率f=1 GHz，單站散射水平極化。網格平均尺寸約為λ/10，未知量e=26 157。采用矩量法、多層快速多極子算法計算，并與商用軟件的矩量法進行對比。幾種算法的計算結果對比如圖8所示，矩量法和多層快速多極子算法的計算結果一致性較好，與商用軟件的計算結果區別較小，具有較高的可信度。

采用基于HMQ全局徑向基函數的域元法對導彈模型進行參數化，設計變量分布如圖9所示。采用基于矩量法、多層快速多極子算法的伴隨方法求解各設計變量在z方向的梯度并與基于矩量法的有限差分計算結果對比，計算中差分擾動量Δx=10-2m。其中設計變量C和D(如圖9所示)在各入射角度的梯度計算結果如圖10和圖11 所示。基于伴隨方法的矩量法和多層快速多極子算法得到的梯度與有限差分方法得到的梯度在趨勢上較為一致。其中，基于伴隨的矩量法與有限差分計算得到的梯度吻合較好，誤差較小；基于多層快速多極子算法的伴隨求解在梯度較小及梯度起伏較為劇烈的區域存在一定誤差。由圖10 和圖11可以看到，雷達散射截面的導數隨入射角度變化在大小和正負號上均存在明顯起伏。圖12為入射角為0°時表面感應電流和伴隨變量的分布云圖。

圖8 導彈單站RCS結果(水平極化)Fig.8 RCS results of missile (horizontal polarization)

圖9 導彈模型表面網格Fig.9 Surface mesh of missile model

圖10 導彈設計變量C導數(水平極化)Fig.10 Gradients of missile design variable C(horizontal polarization)

圖11 導彈設計變量D導數(水平極化)Fig.11 Gradients of missile design variable D(horizontal polarization)

圖12 導彈感應電流及伴隨變量云圖Fig.12 Contours of surface current and adjoint variable of missile

在導數求解過程中，矩量法計算中內存使用10.4 G，計算使用20核并行計算，計算180個入射角度，總計算耗時608 min。多層快速多極子算法內存使用25.4 MB，計算使用8核并行計算，計算用時57 min。

4 結 論

本文提出了一種基于麥克斯韋積分方程離散伴隨方程的RCS梯度計算方法，實現了RCS梯度的高效、高精度求解，為基于梯度的電大尺寸目標高精度氣動/隱身一體化優化提供了基礎。

1) 麥克斯韋積分方程的離散伴隨方程形式與原積分方程一致，可直接采用原積分方程的數值解法如矩量法及多層快速多極子算法求解，程序實現容易。

2) 伴隨方程的計算量與設計變量個數基本無關，可以實現通過2次線性方程組求解得到雷達散射截面關于所有設計變量的梯度信息。伴隨求解的計算量與原方程求解基本一致，存儲量在原方程求解的基礎上增加不明顯。

3) 基于伴隨方程的梯度求解方法具有較高精度，多層快速多極子算法求取的梯度在精度上雖然略低于矩量法的計算結果，但在計算時間和內存需求上均明顯優于矩量法，更適合于電大尺寸復雜問題的評估和優化。