基于混合場積分方程的半空間上方金屬目標電磁散射特性高效分析*

袁倩 周培陽 何姿? 陳學文 丁大志

1) (南京理工大學電子工程與光電技術學院,南京 210094)

2) (華中科技大學物理學院,武漢光電國家研究中心,武漢 430074)

1 引言

半空間環境如地/空、海/空等是研究復雜環境下目標電磁散射特性的基礎[1].工業界和學術領域通過大量的研究,如在地、海背景中的目標探測與識別[2]、雷達目標的電磁隱身設計[3,4]、復雜環境下的通信系統研究及天線設計[5]、地球物理探測與遙感遙測[6,7],對復雜環境中目標的電磁特性進行了理論研究和數值分析,為后期提供可靠高效的分析和改進方法[8].

在科學計算的快速推動下,平面分層結構越來越被人們所關注,因此提出了微分方法和積分方法來處理平面分層結構問題.以時域有限差分法(finite difference time domain method,FDTD)、有限元法(finite element method,FEM)[9-12]為代表的微分類方法用來處理開放問題的電磁仿真.離散網格必須覆蓋研究目標及其周圍的空域,并且必須在距離目標一定距離處用適當的邊界條件截斷,因為其求解區域包含所有目標,導致其計算效率低下.以矩量法(method of moment,MoM)[13]為代表的積分類方法,目標的網格離散時,其未知量相比微分類方法大幅減小.電場積分方程(electric field integral equation,EFIE)因具有良好的普適性可以來處理任意目標,但因其算子條件數不好,迭代收斂慢;盡管磁場積分方程(magnetic field integral equation,MFIE)條件數有所改善,但依然無法處理開放體等非常規目標,且兩者的內諧振問題均無法消除.為了解決上述問題,通常選擇混合場積分方程(combined field integral equation,CFIE)[14],在CFIE 的基礎上,針對半空間金屬目標依然有計算效率低的問題,1909 年,德國物理學家Sommerfeld 使用赫茲勢方法系統研究了有損介質上的水平極化和垂直極化電偶極子[15],提出了索末菲積分.索末菲積分高振蕩,慢收斂,所以計算極其復雜且耗時極長,早期學者主要以漸近近似方法研究索末菲積分,但該方法參數逼近較大,所以應用范圍有限[16].數值計算方法因計算機水平的發展而成為主流方法.早期的研究主要針對簡單的輻射和散射問題而提出公式方法,故系統性理論欠缺.70 年代,King 等[17]就已經分析了有關半空間天線的輻射特性;后來Xu 等[18]人又對半空間二維柱體的散射問題進行了系統研究,對半空間任意形狀物體散射問題和建模的關鍵問題推導出與之相對應的并矢格林函數[19],該問題從Sommerfeld做出開創性工作以來,相關研究一直得到國內外電磁學學者的深入研究.學者把半空間垂直偶極子源推導到任意極化和任意位置,還把半空間背景擴展成多層介質,然后建立了相關的格林函數表達式.到了90 年代初,Michalski 等[20]在前人的基礎上使用傳輸線模型推導出適用于任意分層介質中導體目標的混合位積分方程,還有相關格林函數的譜域形式,因為此方法兼容了自由空間中的混合位方法,所以該方法得到了廣泛的使用.之后Carin 等[21]將此方法應用到埋地目標的電磁建模中,最終得到了很好的效果.隨著半空間目標剖分未知量的增加,Geng 等[22]將多層快速多極子方法引入到分析地面上方或下方目標的電磁散射中.在國內,關于半空間及多層媒質背景下目標電磁散射建模也取得了較大的進展.1991 年,國內南京理工大學的Yang、Chow 和Fang 教授使用Prony 方法和指數方程,推導出了新的離散復鏡像方法(discrete complex image method,DCIM)[23].之后又經過眾多學者對離散復鏡像采樣路徑進行優化,使DCIM 得到廣泛應用,其中包括Aksun[24]提出優化路徑的二階DCIM、Yuan 等[25]提出的直接離散復鏡像法、武漢大學朱國強教授[26]為了避免準靜態項和表面波的提取改進了DCIM 方法,一系列關于DCIM 的工作促進了該方法的發展和應用.還有眾多學者對索末菲積分的研究做了很多顯著工作,諸如電子科技大學潘錦教授等[27]提出的截斷技術、胡俊教授[28]基于快速漢克爾變化提出的加速方法、東南大學洪偉教授等[29]為了增加遠場精度提出的全模式方法、北京理工大學盛新慶教授、吳比翼博士等[30]使用基因遺傳優化算法來提高DCIM 的計算效率等.基于對格林函數特性(柱對稱性、平移不變性和相對緩慢的空間變化特性)的認識,研究學者意識到可以利用插值方法來計算半空間格林函數.文獻[31]就曾將插值方法用于透地天線的仿真,Ling 等[32]討論了用于微帶天線仿真的格林函數插值方法,Chen等[33]也將插值方法用于混合位積分方程(mixed potential integral equation,MPIE)對微帶耦合介質諧振腔的分析.

本文使用CFIE 處理半空間上方金屬目標電磁散射問題,并給出一些有效的方法來克服數值計算上存在的困難.首先給出分析半空間大尺度金屬目標的CFIE,分析了半空間金屬目標散射特性,但是還存在計算效率過低的問題.為了加快求解速度,本文使用了半空間格林函數的列表與空間插值方法對阻抗矩陣生成過程進行優化,從而提高阻抗矩陣的生成速度.該方法在保證精度的同時,能顯著降低求解問題的時間.

2 半空間CFIE

如圖1 所示,半空間環境下任意形狀的金屬目標電磁散射問題,金屬目標在有耗半空間上方.上下半空間的介質參數分別是 (ε1,μ1,σ1) 和 (εh,μh,σh),金屬目標上某點處指向金屬體外部的單位法向量用n1表示,初始入射場用(Einc,Hinc)表示,半空間分界面產生的反射場用(Eref,Href)表示.利用惠更斯等效原理建立S面外部空間區域中場源的等效問題,如圖2 所示,金屬體外等效為外表面S上的等效電流J1在半空間環境中產生的場,J1=n1×H1,其中(E1,H1)表示外表面S上的總場.

圖1 半空間上方金屬目標示意圖Fig.1.Schematic diagram of the metal target above the half space.

圖2 外部等效模型示意圖Fig.2.Schematic diagram of external equivalent model.

外表面S上的總場(E1,H1)[34]為

結合外表面S上的邊界條件,可以建立外部等效問題的電場積分方程EFIE 和磁場積分方程MFIE為[34]

CFIE 能有效避免內諧振現象,并且具有矩陣性態好的優勢,所以本文引入了CFIE,CFIE 的選取在實際中有TENE,NENH,NETH 及TENH四種類型[35].采用的是TENH 型的CFIE,因為其形成的矩陣性態好,并能真正消除內諧振現象.CFIE可以表示為

離散相鄰三角形公共邊上的RWG 基函數,得到如下矩陣方程:

其中In為未知電流向量,m,n=1,2,···,N對應于公共邊編號,Zmn和Vm分別為阻抗矩陣和激勵向量,其元素分別為

?G(kρ)其中,代表譜域格林函數,n=0,1,2 時Sn表示為0 階、1 階和2 階索末菲積分(Sommerfeld integrals,SI)的形式,SI 由于其高振蕩和慢衰減特性,只有采用有效的數值積分方法才能解決直接數值積分耗時巨大的問題.常用的數值積分方法有最陡下降路徑法(steepest descent path,SDP)和離散復鏡像法(DCIM);其他方法包括快速漢克爾變換技術和窗口函數方法.本文采用DCIM 對索末菲積分進行快速計算.DCIM 很好地避免了索末菲積分的直接數值計算,利用一組復鏡像來近似譜域格林函數,之后通過廣義函數束法(generalized pencil of function,GPOF)或者矩陣束方法(matrix pencil,MP) 得到復鏡像的各參數,最后通過索末菲恒等式得到SI 的解析表達式[36]:

3 半空間格林函數的空間插值算法

本文半空間環境下金屬目標采用DCIM 計算格林函數的索末菲積分,但對于大尺度半空間金屬目標,逐點式計算格林函數(在計算面積分的過程中通過場源坐標逐個計算索末菲積分)的方法使得CFIE 的阻抗矩陣填充時間依舊很長,且計算量巨大,顯然僅使用DCIM 并不能解決索末菲積分計算次數多,時間消耗長的問題.考慮到由于半空間環境下格林函數的對稱性和平移不變性(場源相對位置相同,索末菲積分也會相同)的性質,本文給出了一種格林函數高效計算方法—四維空間插值.其基本思路是:首先在一些離散的空間點上進行采樣,用DCIM 方法計算出這些采樣點處的索末菲積分值并建立四維插值表,在阻抗矩陣填充過程中,找到插值表中該點場源相對位置附近的插值節點處的值,進行線性的Lagrange 插值.可以看出,插值方法計算格林函數至少具有以下特點:

1) 降低調用索末菲積分的次數,達到調用次數的最小化;

2) 空間插值有效利用半空間格林函數的柱對稱性和平移不變性,以避免重復多次計算相同位置索末菲積分;

3) 可以通過精確的積分方法與的插值方法相結合來保證格林函數計算精度.

3.1 插值表的建立

對于任何半空間環境下金屬目標,插值表采樣范圍需要目標的空間位置來決定,接著在該范圍離散的坐標點上對索末菲積分進行采樣,建立插值表.對于大尺度半空間金屬目標問題,要建立1個二維插值表[4]用于1個二維插值表和1 個四維插值表用于一般用等間隔采樣方式建立插值表,為了保證插值精度,要求對索末菲積分進行采樣的間隔要足夠小,本文使用0.09 個波長的采樣間隔.

如圖3 所示,表示雷達散射截面(radar cross section,RCS)計算誤差和使用空間插值方法的運行時間隨插值間隔的變化.從圖3 可以看出,當間隔為0.05 個波長時,誤差最小,但運行時間長;隨著插值間隔的增大,誤差變大,時間減少.因此本文選取0.09 個波長的插值間隔,在保證精度的同時也縮短了運行時間.

圖3 RCS 計算誤差和運行時間隨插值間隔變化的折線圖Fig.3.Line chart of RCS calculation error and running time for different interpolation intervals.

圖4 二維插值表示意圖,圖中插值點在二維插值表(ρ,z)中取值Fig.4.Schematic diagram of two-dimensional interpolation table,the interpolation points in the figure are taken in the two-dimensional interpolation table (ρ,z) .

圖5 四維插值表示意圖,圖中插值點在四維插值表(x-x′,y-y′,z-z′,z+z′)中取值Fig.5.Schematic diagram of the 4D interpolation table,the interpolation points in the figure are taken in the four-dimensional interpolation table (x-x′,y-y′,z-z′,z+z′) .

3.2 格林函數插值計算

在阻抗矩陣生成過程中,利用Lagrange 插值公式,可以得到插值表中任意場源點組合對應的索末菲積分值,進一步計算出格林函數.對于和K?e所對應的索末菲積分值Sn,要通過二維插值得到;對于所對應的索末菲積分值Sn,則要進行二維和四維插值.Lagrange 插值方法由于其只需要區間端點的信息并使用線性插值公式,不需要對索末菲積分進行逐次計算,可以獲得最快的計算速度;另外,由于是對場、源點的相對位置進行插值,實際上充分利用了半空間格林函數的柱對稱性和平移性,從而進一步降低計算量.

假設場點和源點分別位于立方體l 和k中(這兩個立方體之間至少相隔一個立方體長度),所有立方體內取一些插值點,然后采用插值近似的方法計算場點r和源點r′之間積分核的相互作用,則可以表示為如下形式:

每一維上的插值基函數[32]表示為

以(x—x′)方向上的插值基函數為例,在(x—x′)方向上有四個插值基函數:

同樣在(y—y′),(z—z′)和(z+z′)方向上各有4 個插值基函數,可以得到此時在插值點(n— 1,j— 1,s— 1,d— 1)處四維張量積為

取出插值表中插值點處響應的函數值,利用Lagrange 插值公式,進一步計算出格林函數值:

4 數值仿真和分析

為了驗證本文方法,本節數值算例運行的平臺為Intel(R) Core(TM) i5-7200U CPU @ 2.50 GHz,內存為12 GB.

算例1首先,如圖6 所示,考慮半空間上方一個金屬立方體.上半空間是空氣,下半空間的相對介電常數為εr=6.38-j0.663,金屬立方體的邊長為1 m,金屬立方體到下平面距離0.2 m,入射平面波為600 MHz,入射角θ=60°,?=0°,散射角θ=60°,?=-180°—180° .分別采用了本文提出的插值加速方法、傳統的CFIE 方法(未插值)、以及商業軟件FEKO 計算了該模型的雷達散射截面.

圖6 金屬立方體模型示意圖Fig.6.Schematic diagram for the metal cube model.

圖7 仿真結果表明,本文提出的加速方法與傳統CFIE 方法以及FEKO 仿真結果均吻合,證明了本文方法的正確性.

圖7 金屬立方體雙站RCS 仿真結果對比圖Fig.7.RCS comparison among the proposed method,the traditional CFIE method and the FEKO for the PEC block model.

算例2如圖8 所示,下圖為半空間上方一個金屬slicy.上半空間是空氣,下半空間的相對介電常數為εr=5.0-j0.2,金屬slicy 的長為2 m,寬為2 m,高為1.5 m,金屬slicy 到下平面距離0.2 m,入射平面波為300 MHz,入射角θ=60°,?=0°,散射角θ=60°,?=-180°—180° .分別采用了本文提出的插值加速方法與傳統的CFIE 方法(未插值)計算了該模型的雷達散射截面.

圖8 金屬slicy 模型示意圖Fig.8.Schematic diagram for the metal slicy model.

圖9 仿真結果表明,本文提出的加速方法與傳統CFIE 方法(未插值)計算結果相吻合.

圖9 金屬slicy 模型雙站RCS 仿真結果對比圖Fig.9.RCS comparison between the proposed method and the traditional CFIE method for the PEC slicy model.

算例3如圖10 所示,下圖為半空間上方一個金屬船.上半空間是空氣,下半空間的相對介電常數為εr=5.0-j0.2,金屬船的長為7 m,寬為1 m,高為1.7 m,金屬船到下平面距離0.2 m,入射平面波為300 MHz,入射角θ=60°,?=0°,散射角θ=60°,?=-180°—180°.分別采用了本文提出的插值加速方法與傳統的CFIE 方法(未插值)計算了該模型的雷達散射截面.

圖11 仿真結果表明,本文提出的加速方法與傳統CFIE 方法(未插值)計算結果吻合.并且,由于激勵、背景和結構的對稱性,所得的RCS 也具有較好的對稱性.

圖11 金屬船雙站RCS 仿真結果對比圖Fig.11.RCS comparison between the proposed method and the traditional CFIE method for the PEC ship model.

表1 所示為本文方法與FEKO、傳統CFIE 的計算資源比較.從表1 可以看出,在計算半空間上方金屬船時,本文方法在傳統CFIE 基礎上提高了18684/1476≈12.6 倍,本文提出的插值加速方法,保證了計算精度的基礎上,比FEKO(CFIE)和傳統CFIE 方法明顯節省了計算時間.與傳統CFIE 方法相比,插值表存儲需要一定的計算內存;另外,FEKO(CFIE)也運用了插值方法和本文插值方法計算內存基本一致.因此,可以通過精確的積分方法與插值方法相結合來保證格林函數的計算精度和效率.隨著金屬目標電尺寸的增大,與傳統CFIE 方法相比,插值加速方法在格林函數計算時間上的減少將越發顯著,對于提高計算效率具有重要意義.

表1 本文方法與FEKO、傳統CFIE 的計算資源比較Table 1.Computational comparison among the proposed method,FEKO and the traditional CFIE.

5 結論

本文討論了半空間環境下大尺度金屬目標電磁散射快速分析方法,在CFIE 的基礎上,利用半空間格林函數的性質,在阻抗矩陣生成過程中,通過二維和四維空間插值方法,避免了索末菲積分的大量重復計算.本文提出的插值加速方法可節省積分的次數,同時,插值加速方法可以建立通用的插值表,從而達到與自由空間問題接近的計算效率.