基于矩陣降維的MoM-PO電磁散射混合算法

李建周, 張海軒, 陳書慧

(西北工業(yè)大學(xué)電子信息學(xué)院, 陜西 西安 710129)

0 引 言

眾所周知,矩量法(method of moments, MoM)是一種基于離散積分方程的數(shù)值方法,廣泛被用于計算具有較小電尺寸的任意三維目標(biāo),通過分解子區(qū)域和選用基函數(shù),構(gòu)造阻抗矩陣進(jìn)行電流的求解計算。但是,當(dāng)未知量增加時,內(nèi)存和時間消耗巨大,計算效率低,無法用于計算電大尺寸目標(biāo)[1-2]。為了減少計算量,加快方程求解,已有文獻(xiàn)提出了諸如高階MoM[3-5]、自適應(yīng)算法[6-10]、快速多極子法[11-14]和預(yù)處理技術(shù)[15-17]等一些改進(jìn)算法,但是也依然存在著不能有效降低計算復(fù)雜度從而完全有效地處理電大尺寸目標(biāo)等問題[18]。

對電大尺寸目標(biāo),高頻算法[18-21]是一種常用的計算方法,但無法有效處理精細(xì)結(jié)構(gòu)目標(biāo)。因此,MoM常常和物理光學(xué)法(physical optic, PO)結(jié)合形成混合算法[22-29]。混合算法將目標(biāo)分成矩量區(qū)和物理光學(xué)區(qū):矩量區(qū)包含精細(xì)的復(fù)雜結(jié)構(gòu),物理光學(xué)區(qū)則包含大尺寸平滑結(jié)構(gòu),以期兼顧兩種算法的精度和效率。文獻(xiàn)[22]中,使用矩量區(qū)表面電流系數(shù)來表征物理光學(xué)區(qū)的電流,該方法不僅忽略了物理光學(xué)區(qū)的自耦合作用,還忽略了物理光學(xué)區(qū)對矩量區(qū)的耦合作用,會影響結(jié)果的精確性。文獻(xiàn)[25]提出的高階混合算法和文獻(xiàn)[27]提出的自適應(yīng)積分方法無法對任意的三維復(fù)雜目標(biāo)進(jìn)行精確的建模和模擬,導(dǎo)致這些方法的適用性不夠好[2]。在文獻(xiàn)[18]提出的迭代混合算法中,先不考慮耦合,單獨(dú)計算矩量區(qū)的電流后,在矩量區(qū)的影響下計算物理光學(xué)區(qū)電流,然后由物理光學(xué)區(qū)電流對激勵進(jìn)行修正后計算矩量區(qū)的新電流,這樣不斷迭代計算,直至兩區(qū)域電流解趨于穩(wěn)定。文獻(xiàn)[27]相比文獻(xiàn)[18]在矩量區(qū)計算中使用快速多極子算法進(jìn)一步降低計算成本,但是在面對電大尺寸矩量區(qū)時,兩區(qū)域之間的大量迭代仍會導(dǎo)致矩量區(qū)的矩陣計算成本很高[29]。

本文提出了一種新的矩陣降維混合算法,該方法立足于對矩陣的分塊和變換,依據(jù)MoM和RWG(Rao-Wilton-Glisson)基函數(shù)展開目標(biāo)表面電流后,對阻抗矩陣依照降維區(qū)域(物理光學(xué)區(qū))和原始區(qū)域(矩量區(qū))進(jìn)行分塊,同時未知量(即基函數(shù)的展開系數(shù))也被分為兩個部分。之后,依照兩區(qū)域之間電流的相互關(guān)系建立兩部分基函數(shù),展開系數(shù)之間的線性關(guān)系,構(gòu)造耦合轉(zhuǎn)移矩陣和激勵轉(zhuǎn)移矩陣,并以此通過對矩陣方程的變形,將對整個阻抗矩陣的元素填充過程轉(zhuǎn)化為較小矩陣塊之間的簡單運(yùn)算過程。同時,為了簡化和加速矩陣塊填充,還引入了等效偶極矩法[30-31],將基函數(shù)視為偶極子模型,避免了矩陣塊填充時的雙重積分計算。通過以上處理,阻抗矩陣的維數(shù)降低為矩量區(qū)自耦合矩陣塊的維數(shù),大大減少了后續(xù)求解計算的復(fù)雜度。在該算法中,耦合轉(zhuǎn)移矩陣和激勵轉(zhuǎn)移矩陣具有明確的物理含義,前者表征了矩量區(qū)對物理光學(xué)區(qū)的耦合作用,后者為入射波對物理光學(xué)區(qū)的激勵,在矩陣的變形過程中還考慮了物理光學(xué)區(qū)對矩量區(qū)的耦合作用,而僅忽略掉物理光學(xué)區(qū)的自耦合。同時,因為對整個目標(biāo)的表面電流進(jìn)行RWG基函數(shù)的展開,矩量區(qū)和物理光學(xué)區(qū)之間的表面電流連續(xù)性也可以得到保證,保證了計算結(jié)果的精確性。算例結(jié)果表明,與傳統(tǒng)MoM算法相比,該算法在保證計算精度的同時大大提高了計算效率。

1 算 法

1.1 矩陣方程降維方法

MoM算法使用RWG基函數(shù)對目標(biāo)表面電流進(jìn)行離散,最終將散射問題轉(zhuǎn)變?yōu)橐曰瘮?shù)展開系數(shù)為未知量的矩陣方程求解問題。求解的計算量主要有兩部分:首先是阻抗矩陣的填充,阻抗矩陣中各個元素的計算需要大量的二重積分運(yùn)算;其次是矩陣方程的求解,高階阻抗矩陣方程的求解運(yùn)算往往耗費(fèi)大量時間。

在矩陣降維混合算法中,首先對電大尺寸目標(biāo)進(jìn)行區(qū)域劃分,如圖1所示。模型在入射電場Ein與入射磁場Hin作用下,一部分(圖1中彈翼和彈頭部分)被稱為原始區(qū)域,含有M個未知量,原始區(qū)域上有電流JM;另一部分(圖1中彈體部分)被稱為降維區(qū)域,含有N個未知量。降維區(qū)域上有電流JN,且M+N=X為未知量總個數(shù)。由尾翼部分細(xì)節(jié)圖可以看出,面元對j和面元對l位于原始區(qū)中,而面元對k位于降維區(qū)中。

圖1 模型區(qū)域劃分及面元剖分Fig.1 Model area division and mesh dissection

按照區(qū)域劃分,可以將矩量法的阻抗矩陣分為4個矩陣塊:

(1)

式中:Zij為阻抗矩陣元素;Ii為電流系數(shù);Vi為入射波激勵。將式(1)簡寫為

(2)

式中:依照MoM算法以及RWG基函數(shù)的物理含義,可以得出Z1是原始區(qū)自耦合矩陣;Z2是降維區(qū)對原始區(qū)的耦合作用矩陣;Z3是原始區(qū)對降維區(qū)的耦合作用矩陣;Z4是降維區(qū)自耦合矩陣;IM、VM分別是原始區(qū)電流展開系數(shù)向量和激勵向量;IN、VN分別是降維區(qū)電流展開系數(shù)向量和激勵向量。為了能夠?qū)?M+N)×(M+N)維的阻抗矩陣降低為M×M維,僅考慮阻抗矩陣的前M行,并通過矩陣運(yùn)算變形為

Z1IM+Z2IN=VM

(3)

然后,構(gòu)造兩組未知量之間的線性關(guān)系,假設(shè)降維區(qū)的N個未知量可以由原始區(qū)的M個未知量線性表示,即

IN=TIM+P

(4)

式中:T為系數(shù)矩陣;P為常數(shù)向量。將式(4)代入式(3),可以簡化為只有M個未知量時的矩陣方程:

Z′IM=V′

(5)

式中:V′=VM-Z2P是新的M維激勵列向量;Z′=Z1+Z2T是新的M×M維阻抗矩陣。由式(2)～式(5)的矩陣變形可以看出,通過構(gòu)造矩陣T和矩陣P,可以實現(xiàn)對原阻抗矩陣的降維,可通過求解式(5)的低維矩陣方程得到原始區(qū)的M個未知量,再由式(4)求出剩余的N個降維區(qū)的未知量,從而得到整個目標(biāo)表面的電流系數(shù)。

相比矩量法中(M+N)×(M+N)的阻抗矩陣填充,在矩陣降維之后,保持不變的部分為M×M維的Z1和M×N維的Z2。而原始區(qū)對降維區(qū)的耦合矩陣Z3,變?yōu)榱送S數(shù)的矩陣T,激勵對降維區(qū)的作用VN則變?yōu)橥S數(shù)的矩陣P。將矩陣T稱為耦合轉(zhuǎn)移矩陣,矩陣P稱為激勵轉(zhuǎn)移矩陣,這兩個矩陣的構(gòu)造是實現(xiàn)本算法的關(guān)鍵。在求解未知量時,因M+N階的高階矩陣Z變?yōu)镸階的低階矩陣Z′,求解計算量將大大削減。在內(nèi)存占用方面,只需存儲阻抗矩陣中前M行的M×N個元素,以及一個N×M的耦合轉(zhuǎn)移矩陣T,相比傳統(tǒng)矩量法可以節(jié)省N×N個元素的存儲空間。

1.2 耦合轉(zhuǎn)移矩陣T和激勵轉(zhuǎn)移向量P的構(gòu)造

圖2 降維區(qū)第k個公共邊上的面元對Fig.2 Pair of triangular patches on the k th common edge in the dimension reduction region

(6)

(7)

另外對pk處的電流進(jìn)行RWG基函數(shù)展開,還可以得到

在美國緬因州，有一所神奇的高中。這所學(xué)校招收的大多數(shù)學(xué)生幾乎沒有希望被錄取進(jìn)入可靠的大學(xué)，更不用說競爭激烈的學(xué)校。但這里98%以上的畢業(yè)生進(jìn)入了四年制大學(xué)，其中包括多所名校。98%的大學(xué)錄取率一直保持了幾十年。本書是兩位作者作為父母、作為教育者幾十年經(jīng)驗和經(jīng)歷的總結(jié)，因而充滿著真知灼見，既有理想的光芒，又有落地的操作性。對家長和教育者來說，本書都是智慧的寶庫，也是實踐的手冊。

(8)

(9)

由式(7)～式(9)可得

(10)

在降維區(qū)內(nèi),從k=1,2,…,N循環(huán)各個公共邊從而得到IN。對比式(4),可以將矩陣T和矩陣P中的元素分別定義為

(11)

(12)

式中:k為波數(shù);m表示等效偶極子矩量,是有效偶極子電流和有效偶極子長度的乘積。綜上,由式(11)和式(12)即可構(gòu)建出矩陣T和矩陣P。

2 算 例

本節(jié)給出幾個數(shù)值算例,通過與FEKO軟件中的MoM算法和PO算法相比,評估本文方法的精度和效率。

2.1 算例1

如圖3所示,模型由圓柱和棱臺組成。圓柱高和直徑均為1 m,底面圓心位于坐標(biāo)原點處。棱臺高0.1 m,上下底面邊長分別為0.05 m和0.1 m,底面中心坐標(biāo)為(0.6,0,0.5)。設(shè)置頻率為2 GHz的θ極化入射波沿著φ=0°、θ=0°～180°方向入射。如圖3(c)所示,將棱臺部分設(shè)置為原始區(qū)域,共剖分出340個面元對,有510個未知量;圓柱部分設(shè)置為降維區(qū)域,剖分面元對38 042個,有57 063個未知量。分別用FEKO中的MoM算法和MoM-PO混合算法計算單站雷達(dá)散射截面積(radar cross section, RCS),并與本文的矩陣降維算法對比,結(jié)果如圖3(d)所示。可以看出,本文算法與FEKO中的矩量法以及混合算法的結(jié)果都吻合得非常好。其中,矩量法用時41 min,而本文的矩陣降維算法用時9 min,計算時間節(jié)省了78.04%,具有很好的計算速度優(yōu)勢。

圖3 圓柱與棱臺組合體及其RCSFig.3 Combination of cylinder and prism and its RCS

2.2 算例2

如圖4所示,模型由電大尺寸導(dǎo)體球、圓錐和正方體組成。其中,導(dǎo)體球半徑為1.3 m,球心坐標(biāo)為(0,0,0.7)。圓錐高為0.2 m,底面半徑為0.1 m,底面圓心坐標(biāo)為(0,0,2.02)。正方體的中心坐標(biāo)在(1.3,0.1,0.1)處,棱長為0.1 m。設(shè)置頻率為1.5 GHz的θ極化入射波沿著φ=0°、θ=0°～180°方向入射。如圖4(c)將圓錐和正方體部分設(shè)置為原始區(qū)域,有742個面元對、1 113個未知量,導(dǎo)體球部分設(shè)置為降維區(qū)域,有93 374個面元對、140 061個未知量。分別用MoM算法、PO和本文方法計算單站RCS,結(jié)果如圖4(d)所示。其中,MoM用時6 h,矩陣降維算法用時為94 min。可以看出,矩陣降維混合算法在計算精度方面優(yōu)于PO方法,與MoM結(jié)果更接近,且計算用時相比MoM大大縮減。

圖4 球、圓錐與正方體組合體及其RCSFig.4 Combination of sphere, cone, and cube and its RCS

2.3 算例3

導(dǎo)彈模型如圖5(a)所示,將導(dǎo)彈的弧形彈頭和8片彈翼設(shè)置為原始區(qū)域,剖分三角面元對34 916個、未知量51 809個;其余彈體部分設(shè)置為降維區(qū)域,剖分三角面元對49 106個、未知量73 296個。設(shè)置頻率為6 GHz的θ極化入射波沿著φ=0°、θ=0°～180°方向入射,使用FEKO中的MoM算法進(jìn)行仿真,將單站RCS結(jié)果與本文的矩陣降維算法進(jìn)行對比,結(jié)果如圖5(b)所示。可以看出,MoM與矩陣降維算法的結(jié)果基本一致,本文算法在面對包含精細(xì)結(jié)構(gòu)的電大尺寸目標(biāo)時,計算精度良好。

圖5 導(dǎo)彈模型及其RCSFig.5 Missile model and its RCS

3 結(jié) 論

本文在傳統(tǒng)MoM算法構(gòu)造阻抗矩陣方程的基礎(chǔ)上,提出了矩陣降維混合算法。通過對阻抗矩陣的分塊,依托于各區(qū)域間的表面電流關(guān)系,通過矩陣方程的變形合理地簡化阻抗矩陣,推導(dǎo)了新的阻抗矩陣方程表達(dá)式,能夠有效地降低阻抗矩陣維數(shù),提高阻抗矩陣填充以及方程求解的計算效率。同時,依托于明確的物理含義,保證計算的結(jié)果與傳統(tǒng)矩量法的精度可比擬。通過若干算例,驗證了矩陣降維混合算法的正確性。通過計算耗時可以看出,該算法與傳統(tǒng)的數(shù)值算法相比,具有計算效率高的優(yōu)點,可以應(yīng)用于計算包含精細(xì)結(jié)構(gòu)的電大尺寸目標(biāo)問題。