袁倩 周培陽 何姿? 陳學文 丁大志
1) (南京理工大學電子工程與光電技術(shù)學院,南京 210094)
2) (華中科技大學物理學院,武漢光電國家研究中心,武漢 430074)
半空間環(huán)境如地/空、海/空等是研究復雜環(huán)境下目標電磁散射特性的基礎(chǔ)[1].工業(yè)界和學術(shù)領(lǐng)域通過大量的研究,如在地、海背景中的目標探測與識別[2]、雷達目標的電磁隱身設(shè)計[3,4]、復雜環(huán)境下的通信系統(tǒng)研究及天線設(shè)計[5]、地球物理探測與遙感遙測[6,7],對復雜環(huán)境中目標的電磁特性進行了理論研究和數(shù)值分析,為后期提供可靠高效的分析和改進方法[8].
在科學計算的快速推動下,平面分層結(jié)構(gòu)越來越被人們所關(guān)注,因此提出了微分方法和積分方法來處理平面分層結(jié)構(gòu)問題.以時域有限差分法(finite difference time domain method,FDTD)、有限元法(finite element method,FEM)[9-12]為代表的微分類方法用來處理開放問題的電磁仿真.離散網(wǎng)格必須覆蓋研究目標及其周圍的空域,并且必須在距離目標一定距離處用適當?shù)倪吔鐥l件截斷,因為其求解區(qū)域包含所有目標,導致其計算效率低下.以矩量法(method of moment,MoM)[13]為代表的積分類方法,目標的網(wǎng)格離散時,其未知量相比微分類方法大幅減小.電場積分方程(electric field integral equation,EFIE)因具有良好的普適性可以來處理任意目標,但因其算子條件數(shù)不好,迭代收斂慢;盡管磁場積分方程(magnetic field integral equation,MFIE)條件數(shù)有所改善,但依然無法處理開放體等非常規(guī)目標,且兩者的內(nèi)諧振問題均無法消除.為了解決上述問題,通常選擇混合場積分方程(combined field integral equation,CFIE)[14],在CFIE 的基礎(chǔ)上,針對半空間金屬目標依然有計算效率低的問題,1909 年,德國物理學家Sommerfeld 使用赫茲勢方法系統(tǒng)研究了有損介質(zhì)上的水平極化和垂直極化電偶極子[15],提出了索末菲積分.索末菲積分高振蕩,慢收斂,所以計算極其復雜且耗時極長,早期學者主要以漸近近似方法研究索末菲積分,但該方法參數(shù)逼近較大,所以應(yīng)用范圍有限[16].數(shù)值計算方法因計算機水平的發(fā)展而成為主流方法.早期的研究主要針對簡單的輻射和散射問題而提出公式方法,故系統(tǒng)性理論欠缺.70 年代,King 等[17]就已經(jīng)分析了有關(guān)半空間天線的輻射特性;后來Xu 等[18]人又對半空間二維柱體的散射問題進行了系統(tǒng)研究,對半空間任意形狀物體散射問題和建模的關(guān)鍵問題推導出與之相對應(yīng)的并矢格林函數(shù)[19],該問題從Sommerfeld做出開創(chuàng)性工作以來,相關(guān)研究一直得到國內(nèi)外電磁學學者的深入研究.學者把半空間垂直偶極子源推導到任意極化和任意位置,還把半空間背景擴展成多層介質(zhì),然后建立了相關(guān)的格林函數(shù)表達式.到了90 年代初,Michalski 等[20]在前人的基礎(chǔ)上使用傳輸線模型推導出適用于任意分層介質(zhì)中導體目標的混合位積分方程,還有相關(guān)格林函數(shù)的譜域形式,因為此方法兼容了自由空間中的混合位方法,所以該方法得到了廣泛的使用.之后Carin 等[21]將此方法應(yīng)用到埋地目標的電磁建模中,最終得到了很好的效果.隨著半空間目標剖分未知量的增加,Geng 等[22]將多層快速多極子方法引入到分析地面上方或下方目標的電磁散射中.在國內(nèi),關(guān)于半空間及多層媒質(zhì)背景下目標電磁散射建模也取得了較大的進展.1991 年,國內(nèi)南京理工大學的Yang、Chow 和Fang 教授使用Prony 方法和指數(shù)方程,推導出了新的離散復鏡像方法(discrete complex image method,DCIM)[23].之后又經(jīng)過眾多學者對離散復鏡像采樣路徑進行優(yōu)化,使DCIM 得到廣泛應(yīng)用,其中包括Aksun[24]提出優(yōu)化路徑的二階DCIM、Yuan 等[25]提出的直接離散復鏡像法、武漢大學朱國強教授[26]為了避免準靜態(tài)項和表面波的提取改進了DCIM 方法,一系列關(guān)于DCIM 的工作促進了該方法的發(fā)展和應(yīng)用.還有眾多學者對索末菲積分的研究做了很多顯著工作,諸如電子科技大學潘錦教授等[27]提出的截斷技術(shù)、胡俊教授[28]基于快速漢克爾變化提出的加速方法、東南大學洪偉教授等[29]為了增加遠場精度提出的全模式方法、北京理工大學盛新慶教授、吳比翼博士等[30]使用基因遺傳優(yōu)化算法來提高DCIM 的計算效率等.基于對格林函數(shù)特性(柱對稱性、平移不變性和相對緩慢的空間變化特性)的認識,研究學者意識到可以利用插值方法來計算半空間格林函數(shù).文獻[31]就曾將插值方法用于透地天線的仿真,Ling 等[32]討論了用于微帶天線仿真的格林函數(shù)插值方法,Chen等[33]也將插值方法用于混合位積分方程(mixed potential integral equation,MPIE)對微帶耦合介質(zhì)諧振腔的分析.
本文使用CFIE 處理半空間上方金屬目標電磁散射問題,并給出一些有效的方法來克服數(shù)值計算上存在的困難.首先給出分析半空間大尺度金屬目標的CFIE,分析了半空間金屬目標散射特性,但是還存在計算效率過低的問題.為了加快求解速度,本文使用了半空間格林函數(shù)的列表與空間插值方法對阻抗矩陣生成過程進行優(yōu)化,從而提高阻抗矩陣的生成速度.該方法在保證精度的同時,能顯著降低求解問題的時間.
如圖1 所示,半空間環(huán)境下任意形狀的金屬目標電磁散射問題,金屬目標在有耗半空間上方.上下半空間的介質(zhì)參數(shù)分別是 (ε1,μ1,σ1) 和 (εh,μh,σh),金屬目標上某點處指向金屬體外部的單位法向量用n1表示,初始入射場用(Einc,Hinc)表示,半空間分界面產(chǎn)生的反射場用(Eref,Href)表示.利用惠更斯等效原理建立S面外部空間區(qū)域中場源的等效問題,如圖2 所示,金屬體外等效為外表面S上的等效電流J1在半空間環(huán)境中產(chǎn)生的場,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]為
結(jié)合外表面S上的邊界條件,可以建立外部等效問題的電場積分方程EFIE 和磁場積分方程MFIE為[34]
CFIE 能有效避免內(nèi)諧振現(xiàn)象,并且具有矩陣性態(tài)好的優(yōu)勢,所以本文引入了CFIE,CFIE 的選取在實際中有TENE,NENH,NETH 及TENH四種類型[35].采用的是TENH 型的CFIE,因為其形成的矩陣性態(tài)好,并能真正消除內(nèi)諧振現(xiàn)象.CFIE可以表示為
離散相鄰三角形公共邊上的RWG 基函數(shù),得到如下矩陣方程:
其中In為未知電流向量,m,n=1,2,···,N對應(yīng)于公共邊編號,Zmn和Vm分別為阻抗矩陣和激勵向量,其元素分別為
?G(kρ)其中,代表譜域格林函數(shù),n=0,1,2 時Sn表示為0 階、1 階和2 階索末菲積分(Sommerfeld integrals,SI)的形式,SI 由于其高振蕩和慢衰減特性,只有采用有效的數(shù)值積分方法才能解決直接數(shù)值積分耗時巨大的問題.常用的數(shù)值積分方法有最陡下降路徑法(steepest descent path,SDP)和離散復鏡像法(DCIM);其他方法包括快速漢克爾變換技術(shù)和窗口函數(shù)方法.本文采用DCIM 對索末菲積分進行快速計算.DCIM 很好地避免了索末菲積分的直接數(shù)值計算,利用一組復鏡像來近似譜域格林函數(shù),之后通過廣義函數(shù)束法(generalized pencil of function,GPOF)或者矩陣束方法(matrix pencil,MP) 得到復鏡像的各參數(shù),最后通過索末菲恒等式得到SI 的解析表達式[36]:
本文半空間環(huán)境下金屬目標采用DCIM 計算格林函數(shù)的索末菲積分,但對于大尺度半空間金屬目標,逐點式計算格林函數(shù)(在計算面積分的過程中通過場源坐標逐個計算索末菲積分)的方法使得CFIE 的阻抗矩陣填充時間依舊很長,且計算量巨大,顯然僅使用DCIM 并不能解決索末菲積分計算次數(shù)多,時間消耗長的問題.考慮到由于半空間環(huán)境下格林函數(shù)的對稱性和平移不變性(場源相對位置相同,索末菲積分也會相同)的性質(zhì),本文給出了一種格林函數(shù)高效計算方法—四維空間插值.其基本思路是:首先在一些離散的空間點上進行采樣,用DCIM 方法計算出這些采樣點處的索末菲積分值并建立四維插值表,在阻抗矩陣填充過程中,找到插值表中該點場源相對位置附近的插值節(jié)點處的值,進行線性的Lagrange 插值.可以看出,插值方法計算格林函數(shù)至少具有以下特點:
1) 降低調(diào)用索末菲積分的次數(shù),達到調(diào)用次數(shù)的最小化;
2) 空間插值有效利用半空間格林函數(shù)的柱對稱性和平移不變性,以避免重復多次計算相同位置索末菲積分;
3) 可以通過精確的積分方法與的插值方法相結(jié)合來保證格林函數(shù)計算精度.
對于任何半空間環(huán)境下金屬目標,插值表采樣范圍需要目標的空間位置來決定,接著在該范圍離散的坐標點上對索末菲積分進行采樣,建立插值表.對于大尺度半空間金屬目標問題,要建立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′) .
在阻抗矩陣生成過程中,利用Lagrange 插值公式,可以得到插值表中任意場源點組合對應(yīng)的索末菲積分值,進一步計算出格林函數(shù).對于和K?e所對應(yīng)的索末菲積分值Sn,要通過二維插值得到;對于所對應(yīng)的索末菲積分值Sn,則要進行二維和四維插值.Lagrange 插值方法由于其只需要區(qū)間端點的信息并使用線性插值公式,不需要對索末菲積分進行逐次計算,可以獲得最快的計算速度;另外,由于是對場、源點的相對位置進行插值,實際上充分利用了半空間格林函數(shù)的柱對稱性和平移性,從而進一步降低計算量.
假設(shè)場點和源點分別位于立方體l 和k中(這兩個立方體之間至少相隔一個立方體長度),所有立方體內(nèi)取一些插值點,然后采用插值近似的方法計算場點r和源點r′之間積分核的相互作用,則可以表示為如下形式:
每一維上的插值基函數(shù)[32]表示為
以(x—x′)方向上的插值基函數(shù)為例,在(x—x′)方向上有四個插值基函數(shù):
同樣在(y—y′),(z—z′)和(z+z′)方向上各有4 個插值基函數(shù),可以得到此時在插值點(n— 1,j— 1,s— 1,d— 1)處四維張量積為
取出插值表中插值點處響應(yīng)的函數(shù)值,利用Lagrange 插值公式,進一步計算出格林函數(shù)值:
為了驗證本文方法,本節(jié)數(shù)值算例運行的平臺為Intel(R) Core(TM) i5-7200U CPU @ 2.50 GHz,內(nèi)存為12 GB.
算例1首先,如圖6 所示,考慮半空間上方一個金屬立方體.上半空間是空氣,下半空間的相對介電常數(shù)為εr=6.38-j0.663,金屬立方體的邊長為1 m,金屬立方體到下平面距離0.2 m,入射平面波為600 MHz,入射角θ=60°,?=0°,散射角θ=60°,?=-180°—180° .分別采用了本文提出的插值加速方法、傳統(tǒng)的CFIE 方法(未插值)、以及商業(yè)軟件FEKO 計算了該模型的雷達散射截面.
圖6 金屬立方體模型示意圖Fig.6.Schematic diagram for the metal cube model.
圖7 仿真結(jié)果表明,本文提出的加速方法與傳統(tǒng)CFIE 方法以及FEKO 仿真結(jié)果均吻合,證明了本文方法的正確性.
圖7 金屬立方體雙站RCS 仿真結(jié)果對比圖Fig.7.RCS comparison among the proposed method,the traditional CFIE method and the FEKO for the PEC block model.
算例2如圖8 所示,下圖為半空間上方一個金屬slicy.上半空間是空氣,下半空間的相對介電常數(shù)為εr=5.0-j0.2,金屬slicy 的長為2 m,寬為2 m,高為1.5 m,金屬slicy 到下平面距離0.2 m,入射平面波為300 MHz,入射角θ=60°,?=0°,散射角θ=60°,?=-180°—180° .分別采用了本文提出的插值加速方法與傳統(tǒng)的CFIE 方法(未插值)計算了該模型的雷達散射截面.
圖8 金屬slicy 模型示意圖Fig.8.Schematic diagram for the metal slicy model.
圖9 仿真結(jié)果表明,本文提出的加速方法與傳統(tǒng)CFIE 方法(未插值)計算結(jié)果相吻合.
圖9 金屬slicy 模型雙站RCS 仿真結(jié)果對比圖Fig.9.RCS comparison between the proposed method and the traditional CFIE method for the PEC slicy model.
算例3如圖10 所示,下圖為半空間上方一個金屬船.上半空間是空氣,下半空間的相對介電常數(shù)為εr=5.0-j0.2,金屬船的長為7 m,寬為1 m,高為1.7 m,金屬船到下平面距離0.2 m,入射平面波為300 MHz,入射角θ=60°,?=0°,散射角θ=60°,?=-180°—180°.分別采用了本文提出的插值加速方法與傳統(tǒng)的CFIE 方法(未插值)計算了該模型的雷達散射截面.
圖11 仿真結(jié)果表明,本文提出的加速方法與傳統(tǒng)CFIE 方法(未插值)計算結(jié)果吻合.并且,由于激勵、背景和結(jié)構(gòu)的對稱性,所得的RCS 也具有較好的對稱性.
圖11 金屬船雙站RCS 仿真結(jié)果對比圖Fig.11.RCS comparison between the proposed method and the traditional CFIE method for the PEC ship model.
表1 所示為本文方法與FEKO、傳統(tǒng)CFIE 的計算資源比較.從表1 可以看出,在計算半空間上方金屬船時,本文方法在傳統(tǒng)CFIE 基礎(chǔ)上提高了18684/1476≈12.6 倍,本文提出的插值加速方法,保證了計算精度的基礎(chǔ)上,比FEKO(CFIE)和傳統(tǒng)CFIE 方法明顯節(jié)省了計算時間.與傳統(tǒng)CFIE 方法相比,插值表存儲需要一定的計算內(nèi)存;另外,FEKO(CFIE)也運用了插值方法和本文插值方法計算內(nèi)存基本一致.因此,可以通過精確的積分方法與插值方法相結(jié)合來保證格林函數(shù)的計算精度和效率.隨著金屬目標電尺寸的增大,與傳統(tǒng)CFIE 方法相比,插值加速方法在格林函數(shù)計算時間上的減少將越發(fā)顯著,對于提高計算效率具有重要意義.
表1 本文方法與FEKO、傳統(tǒng)CFIE 的計算資源比較Table 1.Computational comparison among the proposed method,FEKO and the traditional CFIE.
本文討論了半空間環(huán)境下大尺度金屬目標電磁散射快速分析方法,在CFIE 的基礎(chǔ)上,利用半空間格林函數(shù)的性質(zhì),在阻抗矩陣生成過程中,通過二維和四維空間插值方法,避免了索末菲積分的大量重復計算.本文提出的插值加速方法可節(jié)省積分的次數(shù),同時,插值加速方法可以建立通用的插值表,從而達到與自由空間問題接近的計算效率.