基于伴隨方程和自動微分的雷達散射截面表面靈敏度計算

周琳，陳憲，黃江濤，*，鄧俊，高正紅

1.中國空氣動力研究與發展中心，綿陽 621000

2.西北工業大學 航空學院，西安 710072

雷達隱身技術是提高目標生存能力的技術，低可探測特征是未來軍用飛機應具備的重要特征之一［1-2］。隨著反隱身技術的發展，飛行器的氣動/隱身設計受到越來越多重視。基于伴隨方程的雷達散射截面（RCS）梯度求解可以通過一次控制方程求解、一次伴隨方程求解實現目標靈敏度的快速計算，是氣動/隱身精細化設計的基礎。

近年來，基于電磁場積分方程伴隨方程的靈敏度計算方法及相應的梯度優化受到廣泛關注［3-7］。基于伴隨方程的電磁場優化早期主要應用在天線設計領域，Georgieva 等［8］在2002 年將伴隨方程引入電磁場頻域計算，研究了基于矩量法的伴隨方程求解，并基于伴隨方法對Yagi-Uda 天線進行了優化；Toivanen 等［9］將自動微分法引入基于矩量法的電磁場伴隨方程的設計變量靈敏度求解中；Kataja 等［10-11］基于矩量法（MOM）［12］、多層快速多極子算法（MLFMA）［13-14］推導了電場積分方程的連續伴隨表達式。

基于伴隨方法的RCS 梯度計算在飛行器隱身優化問題中研究較晚，公開文獻較少。Bondeson 等［15］推導了二維積分形式麥克斯韋方程的連續伴隨方程，并對簡單外形的RCS 進行了優化；周琳等［3-5］在MOM、MLFMA 的基礎上推導了混合場積分方程的離散伴隨方程，并對典型飛翼布局進行了氣動隱身優化。Li 等［6］在MOM 的基礎上提出一種基于伴隨方程和自動微分的高效設計變量梯度求解方法，并對機翼進行了氣動隱身優化。

雖然近年來基于伴隨方程的氣動隱身優化取得了較大進展，但尚無公開文獻探討RCS 的表面靈敏度求解方法。表面靈敏度可以反映表面網格位置的改變對目標函數的影響，一方面對優化設計有重要指導作用；另一方面設計變量的梯度可以根據表面靈敏度直接求出，從而克服經典基于伴隨方程的設計變量RCS 梯度求解需反復填充阻抗矩陣的問題［3-4，8］，在經典伴隨方程梯度求解的基礎上進一步提升設計變量梯度的求解效率。飛行器氣動優化領域伴隨方法已較為成熟，已采用連續伴隨方程［16］、自動微分方法［17］、網格伴隨方程等方法［18-19］對目標的表面靈敏度的計算方法進行了深入研究，相對于氣動問題，電磁場求解問題的阻抗矩陣是稠密矩陣，對于表面網格量龐大的復雜問題，計算和存儲開銷較大，且電磁場問題中的阻抗積分項與幾何信息的依賴關系復雜，運算多為復數運算，難以手動推導得到。以上2 點都使電磁場問題的表面靈敏度求解更具挑戰。

本文在矩量法的基礎上提出基于自動微分的雷達散射截面表面靈敏度計算方法，以及適用于表面靈敏度存儲的稀疏矩陣存儲技術。該方法將基于自動微分法的阻抗矩陣偏導數求解嵌入阻抗矩陣并行填充循環中。對于任意數量設計變量和入射角度，將需進行一次阻抗矩陣、阻抗矩陣偏導數求解，從而使RCS 梯度計算時間與設計變量數無關；當入射角度改變時，計算所有設計變量梯度約為8 次矩陣乘，顯著提升基于伴隨的梯度求解效率。

1 基于伴隨方程的RCS 梯度計算

典型雷達散射截面σ的減縮問題可表示為

式中：I為表面感應電流的基函數系數；D為設計變量向量；V-ZI=0 為離散形式的麥克斯韋方程組。

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

式（2）對設計變量D求導，有

令式（3）右端第1 項為0 即可得到伴隨方程，又由

伴隨方程可以表示為

求式（5）得到伴隨變量Λ后，設計變量的RCS 梯度可以表示為

設計變量D的梯度也可以根據表面靈敏度間接求解。將式（6）中的D用r替換得到基于伴隨方程的表面靈敏度計算方法

設計變量D與表面網格r的關系可以表示為

式中：Np為表面網格點數量；為RCS關于表面網格點rj=［x，y，z］在3 個坐標方向的靈敏度；為表面網格點rj與設計變量的關系，由網格變形算法決定。由于設計變量位置、入射角度的改變不影響?Z/?r的結果，因此梯度計算過程中只需要一次?Z/?r求解。而且基于表面靈敏度的設計變量梯度求解方法不需要反復填充阻抗矩陣Z，矩陣-矢量乘積Z的計算次數少，從而提升梯度求解計算效率。

2 RCS 表面靈敏度存儲方法

采用表面靈敏度計算設計變量梯度的關鍵在于表面靈敏度的計算、存儲及表面靈敏度與感應電流的矩陣-矢量乘。本節首先針對表面靈敏度存儲量大的問題提出了表面靈敏度矩陣的稀疏存儲技術。

式中：Ne為基函數個數。由于σ和V的求解表達式較為簡潔，不涉及復雜的積分計算，因此?σ/?r、?V/?r可以通過自動微分直接求解，通過向量-向量相乘計算ΛT?V/?r，計算量較小。式（9）的難點在于?Z/?r的求解和的計算。理論上?Z/?r是維度為Ne×Ne×Np的大規模三維矩陣，無法直接存儲，本節通過分析?Z?r的特點，提出一種點索引稀疏矩陣存儲技術，將的計算量和存儲量降到優化可承受的范圍。

以電磁場積分方程求解最常使用的RWG（Rao-Wilton-Glisson）基函數［20］為例進行推導，RWG 基函數定義在一對相鄰三角形的公共邊上（見圖1），第n條邊對應的基函數fn(r)為

圖1 基函數RWG 示意圖［21］Fig.1 RWG basis function diagram［21］

由RWG 基函數的特點可知，每一阻抗矩陣元素Zmn至多只與8 個表面網格點的位置有關（圖2）。因此，雖然?Z/?r理論上是維度為Ne×Ne×Np的大型三維矩陣，但?Z/?r中非零元素數量只有Ne×Ne×8，因此可以建立以網格點為索引的稀疏矩陣存儲技術，從而顯著降低?Z/?r的存儲需求。

圖2 Zmn存儲示意圖Fig.2 Zmn storage diagram

圖3 雙錐體外形及設計變量位置Fig.3 Design variables of double-ogive test case

在程序處理中，開辟整型索引數組Zdiff_index（Ne，Ne，8），并按圖2 中Index 數組中點的記錄順序存儲與阻抗元素Zmn有關的8 個網格點的全局編號。開辟復數數組Zdiff（Ne，Ne，8）存儲Zmn關于每個相關網格點的導數?Zmn/?r。若與基函數m、n相關的4 個三角形間存在重合頂點，則與基函數Zmn相關的網格點數量將小于8，此時在索引數組Zdiff_index 和Zdiff 數組的相應位置補0。?V/?r的稀疏存儲方法與此非常相似，不再贅述。

3 RCS 表面靈敏度計算方法

以電場積分方程為例詳細分析表面靈敏度的計算技術。為將?Z/?x的計算嵌入阻抗矩陣Z的填充循環［22］，從而利用Z計算的中間結果，實現?Z/?r的高效求解，首先將Zmn的求解拆分成乘積項M1、M2和M2的組合，并將乘積項進一步拆分成積分項G 和H的組合。其中乘積項的表達式為

積分項G和H的表達式為

當場三角形與源三角形重合時積分項需采用奇異積分計算方法，不重合時可以采用高斯積分法。則的線性組合。

在阻抗矩陣Zmn填充循環的基礎上采用自動微分法求解?Z/?r，采用Tapenade 自動微分工具［23］。自動微分方法通過分析程序源碼，將復雜計算分解為簡單運算的組合，根據鏈式求導法則得到復雜目標的梯度，自動微分法不產生截斷誤差，避免了差分法依賴步長的問題。?Z/?r的具體程序求解方法見算法1。其中最多與一條公共邊相連的2 個三角形有關（6 個網格點），當考慮積分項在r=［x，y，z］3 個方向的導數時，分別為維度為1×3、3×3 的矩陣。

?M1/?r、?M2/?r、?M3/?r最多與一條公共邊相連的2 個三角形有關（6 個網格點），根據鏈式求導法則，結合?G/?r、?H/?r的求解結果，可以得到

采用式（20）～式（22）結合算法1 可以得到?Z/?r在r=［x，y，z］3 個方向上的導數。由于實際優化中通常重點關注目標函數在z方向的梯度，因此程序實現時可以僅對z方向的梯度進行計算和存儲，使向量-向量乘和矩陣-向量乘為標量乘，從而進一步降低計算量和存儲量。同時，采用自動微分工具計算積分項和乘積項的偏導數時，需重復計算積分項和乘積項求解的中間量，可以通過修改自動微分得到的源碼，降低重復計算量，從而進一步降低計算量。

完成?Z/?r的計算和存儲后，需根據式（9）求解目標函數對網格位置的梯度dσ/dr，其核心是基于稀疏矩陣存儲格式計算對于每個網格點ri，可以寫為

從算法1 可以看到，提出的表面靈敏度計算方法中矩陣元素填充方法與矩量法基本一致，因此可以直接采用矩量法的并行策略［22］，較易實現并行求解，從而進一步提高計算效率。激勵向量的自動微分梯度與阻抗矩陣的梯度計算非常類似，這里不再贅述。

4 表面靈敏度計算量分析

本節重點關注基于經典伴隨方法和基于表面靈敏度方法求解設計變量梯度的計算量。工程問題中的隱身優化目標通常為給定角域的RCS 均值，記優化問題考慮的設計變量數為N，考慮的入射角度數為Na。2 種方法的計算量對比見表 1。

表1 設計變量梯度計算量對比Table 1 Gradient calculation amount comparison of design variables

梯度計算過程中需要擾動設計變量，生成新外形表面網格并計算新外形的阻抗矩陣Z(D+ΔD)，并進行矩陣-向量乘、向量-向量乘得到ΛTZ(D+ΔD)。由于每個入射角度的感應電流不同，因此計算N個設計變量在Na個入射角度的梯度需要進行N次表面網格變形計算、N次矩陣元素填充、N×Na次矩陣-矢量相乘運算、N×Na次矢量-矢量相乘運算。

對于基于表面靈敏度的梯度計算方法，表面靈敏度計算過程時需采用自動微分法計算?Z/?r，計算量大于阻抗矩陣填充，通過修改自動微分代碼，復用阻抗矩陣填充的中間量可以提升?Z/?r的計算效率。由于?Z/?r與設計變量數量、入射角度無關，在梯度計算過程中僅需一次?Z/?r求解，因此雖然?Z/?r的計算量較大，但對梯度計算效率影響不大。?Z/?r所需的存儲量約為阻抗矩陣存儲量的8 倍。得到?Z/?r后計算該過程需要8Na次矩陣-向量相乘和8Na次向量-向量相乘運算。得到表面靈敏度后，根據式（8）計算設計變量梯度，計算時需依次擾動設計變量，由網格變形算法（如自由曲面造型法）進行表面網格變形計算得到dr/dD，與網格靈敏度dσ/dr相乘得到設計變量的梯度dσ/dD。在基于伴隨方程的氣動優化中，通常引入網格伴隨方程［18-19］實現空間網格的dr/dD高效求解。由于基于麥克斯韋積分方程的RCS 求解僅需要對物體的表面進行網格劃分，計算所需的網格量遠小于氣動求解，采用有限差分法計算dr/dD計算量較小。

5 算例驗證

本節采用雙錐體、球體、F22 外形和飛翼模型對提出的表面靈敏度計算方法進行驗證，并探討典型飛行器的表面敏感性特征。下文算例均以x-y平面為工作平面，計算RCS 在z方向的表面靈敏度。

1）雙錐體模型

雙椎體模型長0.190 5 m，由2 個錐角分別為22.62°、46.4°的半錐拼接而成（見圖 3），是EMCC（Electromagnetic Code Consortium）［24-25］提供用于校驗計算電磁學代碼的標準算例之一。采用雙錐體算例對采用的電磁場求解器和提出的表面靈敏度計算方法進行驗證。雙錐體算例計算采用的雷達波頻率f=9 GHz，單站散射垂直極化（VV），定義沿x軸正向入射時入射角φ=0°，沿x軸負向入射時φ=180°。計算采用的網格平均尺寸約為λ/10，未知量總數N=11 094。

圖4 為采用矩量法求解器的計算結果與實驗結果［25］的對比，可見計算結果與實驗結果吻合較好，證明采用的求解器具有較高精度，滿足RCS評估和靈敏度計算的需求。

圖4 雙錐體外形RCS（f=9 GHz，VV 極化）Fig 4 RCS of double-ogive test case（f=9 GHz，VV）

采用域元法進行表面網格變形，采用表面靈敏度計算得到的設計變量梯度與有限差分計算結果對比見圖5。可以看到，2 種梯度計算方法得到的結果具有良好的一致性，證明提出的表面靈敏度計算方法能夠通過表面靈敏度的計算結果獲取任意設計變量的梯度。

圖5 表面靈敏度計算梯度與有限差分結果對比Fig.5 Comparisons of design variable gradient calculation by surface sensitivity and finite difference

圖6 雙錐體外形表面靈敏度（f=9 GHz，VV 極化）Fig.6 Surface sensitivity of double-ogive（f=9 GHz，VV）

圖7 表面網格點梯度對比Fig.7 Gradients of surface mesh

圖8 球體算例的表面電流和表面靈敏度（f=1 GHz，VV）Fig.8 Surface current and surface sensitivity（f=1 GHz，VV）

圖9 表面靈敏度計算中間量分布（f=1 GHz，VV）Fig.9 Intermediate quantity results（f=1 GHz，VV）

圖10 類F22 戰斗機表面靈敏度（沿鼻錐方向入射）Fig.10 Surface sensitivity of F22（incident towards nose）

圖11 典型飛翼布局表面靈敏度（沿鼻錐方向入射）Fig.11 Surface sensitivity of flying wing（incident towards nose）

采用本文方法計算φ=0°，90°的表面靈敏度見圖 6。與表面電流的分布特點相似，單一入射頻點和入射角度的表面靈敏度具有正負交替的條紋，且正負條紋間隔約為λ4。

對比采用經典伴隨方法和表面靈敏度方法求解設計變量梯度的計算效率和內存需求。計算采用Intel i7-9700 處理器8 核心并行計算。在計算效率方面，雙錐體外形的阻抗矩陣Z填充用時11 s；僅計算z方向靈敏度時填充?Z/?r和Z矩陣用時35 s；除去求解控制方程和伴隨方程的耗時，采用經典伴隨方法計算36 個設計變量在180 個入射角度的梯度用時413 s，采用基于表面靈敏度的梯度計算方法用時19 s。在內存開銷方面，經典矩量方法內存占用1.78 G，表面靈敏度方法內存占用16.72 G。可以看到，雖然偏導數矩陣?Z/?r的填充速度較慢，但由于基于表面靈敏度方法求解設計變量梯度時不需要填充擾動外形的阻抗矩陣，因此當設計變量較多時效率具有優勢。

2）球體模型

球體模型半徑r=0.2 m，計算入射頻率f=1 GHz，單站散射垂直極化（VV），表面網格尺寸Δx=0.05 m；表面網格的總數N=197。

電磁波沿x軸負向入射時，采用經典伴隨法（Adjoint_diff）和本文表面靈敏度計算方法（Adjoint_surface）求解模型在z方向的表面靈敏度dσ/dz，計算結果對比見圖 7。可以看出2 種方法計算得到的梯度完全重合，證明本文表面靈敏度方法具有較高精度。表面感應電流和表面靈敏度計算結果見圖 8，表面靈敏度空間分布特征和表面電流存在一定差異，對于表面感應電流，入射波直接照射的區域明顯高于其他區域，而表面靈敏度的幅值與電流強度不直接相關，感應電流強的區域不一定具有更大的表面靈敏度，因此在優化過程中難以直接根據電流的分布特點判斷區域的RCS 減縮潛力和方向，需進行表面靈敏度求解。

為進一步探討表面靈敏度的特性，從式（9）可以看到，表面靈敏度dσ/dz的計算結果由?σ/?z、線性疊加得到，計算電磁波沿x軸負向入射、垂直極化時每一部分的空間分布特征，計算結果如圖 9。可以看到，3 部分計算結果均具有明顯的相位特征，空間分布與dσ/dz基本一致，與感應電流有差異顯著。

3）F22 戰斗機外形

分析典型戰斗機的RCS 表面靈敏度特征，采用F22 戰斗機外形，計算模型機身總長約19 m，半展長約7 m。考慮入射頻率f=200 MHz，計算網格未知量數N=27 595。

沿鼻錐方向入射時水平極化（HH）和垂直極化（VV）的表面靈敏度計算結果見圖 10。可以看出，與雙錐體和球體算例相似，2 種極化方式入射時表面靈敏度均具有正負相間的特征，但2 種極化的表面靈敏度空間分布存在較大差異。對于垂直極化，戰斗機機身的表面靈敏度大于機翼的表面靈敏度，且厚度較大的區域（如座艙）靈敏度更高，平尾邊緣的靈敏度較大。對于水平極化，靈敏度較大的區域主要為機身和機翼的邊緣，平尾的表面靈敏度較小，水平極化的靈敏在幅值上整體略小于垂直極化的表面靈敏度，且與幾何的厚度沒有明顯關聯。表面靈敏度的計算結果表明，靈敏度的空間分布與極化方式關聯明顯，極化方式的RCS 減縮需重點關注的區域有所區別。

4）典型飛翼外形

計算采用的飛翼布局根弦長約12 m，半展長約10 m，采用NACA65016 對稱翼型生成。考慮入射頻率f=200 MHz，計算網格的未知量總數N=28 086，分別計算入射波為水平極化（HH）、垂直極化（VV）時的表面靈敏度。

沿鼻錐方向入射時的表面靈敏度計算結果見圖 11，其中水平極化的靈敏度明顯高于垂直極化的靈敏度，2 種極化方式的靈敏度幅值的差異主要受當前狀態的RCS 影響，當沿鼻錐入射時，兩種極化的RCS 分別為RCSHH=-4.9 dB、RCSVV=-21.0 dB，RCS 與靈敏度正相關。垂直極化靈敏度較大的區域主要為中央體上下表面，厚度越大的區域靈敏度越大，該現象與F22戰斗機的計算結果一致；而水平極化靈敏度較大的區域為機翼的前后緣，中央體前緣的靈敏度顯著高于其他區域。

圖12 分別為y=6 m 站位VV、HH 極化方式入射時表面靈敏度的計算結果，由于布局采用對稱翼型，2 種極化方式入射時上下表面的靈敏度符號相反、幅值幾乎一致，且都表現出了強烈的正負震蕩特性，正負波峰的間隔約為入射波長的1/4。

圖12 y=6 m 表面靈敏度分布Fig.12 Surface sensitivity at y=6 m

圖13 為入射波水平極化（HH）時前向威脅扇區±30°角域的平均感應電流和表面靈敏度。可以看到，機翼前緣的平均感應電流較強，從前緣向后緣電流幅度逐漸降低。角域平均靈敏度仍表現出正負相間的特征，但由于不同角度入射時表面靈敏度存在相位疊加，平均表面靈敏度的明暗條紋不如單一角度入射時顯著。圖 14 為Slice 1、Slice 2 剖面（見圖 11）的上下表面設計變量在前向威脅扇區的平均梯度，其中網格變形采用自由曲面造型技術（FFD），設計變量位置如圖 11。雖然表面靈敏度呈現出顯著正負交替的特征，設計變量梯度變化較為和緩。上下表面設計變量梯度不完全對稱是由于設計變量上下表面在z方向的位置不完全對稱造成。

圖13 角度平均感應電流和平均表面靈敏度（f=200 MHz，HH）Fig.13 Angle-averaged current and sensitivity（f=200 MHz，HH）

圖14 Slice 1 和Slice 2 上下表面設計變量平均梯度Fig.14 Angle-averaged design variable gradients at Slice 1 and Slice 2

6 結論

提出了一種基于伴隨方程和自動微分的雷達散射截面表面靈敏度計算方法，并結合基函數特點提出了適用于表面靈敏度的稀疏矩陣存儲和矩陣-矢量相乘方法，降低了表面靈敏度的計算、存儲開銷。該方法可以通過一次阻抗矩陣偏導數求解得到所有網格點的表面靈敏度，當改變入射角度時，求解所有設計變量梯度的計算量約為8 次矩陣-矢量乘，極大提高了基于伴隨方程的梯度計算效率。

計算典型入射條件下雙錐體、球體、F22 戰斗機、典型飛翼布局的表面靈敏度，證明了提出的表面靈敏度計算方法的準確性，該方法能夠通過表面靈敏度快速完成設計變量梯度計算。表面靈敏度的空間分布特征受極化方式影響顯著，具有明顯的相位特征，可以為隱身優化設計提供指導。

本文僅研究了基于矩量法的表面靈敏度計算方法，但是矩量法的計算量和存儲量較大，限制了其在電大尺寸問題中的應用，后續將進一步研究將該方法擴展到基于多層快速多極子算法的RCS 表面靈敏度求解中。