一種基于壓縮感知的三維導(dǎo)體目標(biāo)電磁散射問題的快速求解方法

丁亞輝 孫玉發(fā) 朱金玉

(安徽大學(xué)電子信息工程學(xué)院,合肥 230601)

1 引 言

矩量法[1?3]是求解導(dǎo)體目標(biāo)電磁散射問題的常用方法,具有計(jì)算精度高、所用格林函數(shù)自動滿足輻射條件、無須額外設(shè)置邊界條件等優(yōu)點(diǎn).但是矩量法在求解電大尺寸問題時(shí),會因?yàn)榛瘮?shù)個(gè)數(shù)過多、矩陣方程求解的復(fù)雜度過大導(dǎo)致計(jì)算時(shí)間難以接受.為此,有學(xué)者提出將壓縮感知[4?6](compressed sensing,CS)引入矩量法.一種方法是通過引入稀疏變換矩陣實(shí)現(xiàn)對未知電流的稀疏表示,從而在壓縮感知理論框架下構(gòu)造欠定方程,并優(yōu)化求解[7,8].這種方法已應(yīng)用到旋轉(zhuǎn)對稱目標(biāo)[9]和小波矩量法[10]中.另一種方法是首先構(gòu)建富含空間信息的激勵源,然后應(yīng)用矩量法計(jì)算得到表面感應(yīng)電流的觀測值,進(jìn)而可通過恢復(fù)算法獲得寬角度目標(biāo)表面感應(yīng)電流[11?13].本文只討論第一種方法,即構(gòu)造欠定方程的方法.此方法只能用于二維或二維半的物體上,對三維問題效果不佳,原因是三維目標(biāo)通過軟件離散化建模,場源非有序分布,導(dǎo)致目標(biāo)表面電流在常用的稀疏基上不稀疏.而一般獲取稀疏基的字典訓(xùn)練算法又因?yàn)闆]有訓(xùn)練集無法應(yīng)用在此方法上.

本文對上述方法加以改進(jìn),引入特征基函數(shù)[14?16]作為稀疏基,壓縮感知恢復(fù)算法采用廣義正交匹配追蹤算法(GOMP)[17],并計(jì)算了三維導(dǎo)體目標(biāo)的雷達(dá)散射截面(RCS),數(shù)值結(jié)果證明了該方法的準(zhǔn)確性與高效性.

2 基本原理

2.1 壓縮感知構(gòu)造欠定方程

為了提高矩量法的計(jì)算效率,引入信號處理領(lǐng)域中壓縮感知理論來降低內(nèi)存需求,減少電流求解時(shí)間.

在求解目標(biāo)表面電流時(shí),將阻抗矩陣視為測量矩陣,待求解的表面電流視為要重構(gòu)的信號,激勵電壓視為測量值.如果阻抗矩陣滿足約束等距特性[18],便可以從測量值恢復(fù)出原始信號.而壓縮感知中恢復(fù)算法只需用到低維的測量值,所以在應(yīng)用中降低了電流求解時(shí)間.

根據(jù)矩量法,可得到矩陣方程

式中Z,J,V分別是阻抗矩陣、電流向量和激勵向量,N是未知量的數(shù)目.通過對阻抗矩陣隨機(jī)抽取M行構(gòu)造欠定方程

通過應(yīng)用稀疏變換矩陣Ψ使信號J在變換域上是稀疏的,即

A是信號J在Ψ上的變換向量,其稀疏度為K.則(2)式可表示為

其中,再用恢復(fù)算法恢復(fù)出電流.本文采用GOMP方法作為恢復(fù)算法.與正交匹配追蹤(OMP)相比,GOMP在選擇原子時(shí)每次選擇最大的S個(gè)(本文取S=K/4),而OMP只選擇一個(gè),所以GOMP與OMP精度相同,但速度更快.

CS方法需要事先知道目標(biāo)表面電流的稀疏基,對于二維導(dǎo)體目標(biāo),由于目標(biāo)表面電流的變化緩慢,所以可以使用離散余弦變化基.而三維導(dǎo)體目標(biāo)表面劃分為Rao-Wilton-Glisson(RWG)基函數(shù),其目標(biāo)表面電流在一般的正交基下不稀疏,使得CS方法的應(yīng)用受到了限制.

2.2 特征基函數(shù)作為稀疏基

特征基函數(shù)法(CBFM)是一種基于區(qū)域分解的方法.特征基函數(shù)法首先將目標(biāo)劃分為若干個(gè)子域,在每個(gè)子域上構(gòu)造主要特征基函數(shù),并通過次要特征基函數(shù)來反映子域之間的相互作用,以此來實(shí)現(xiàn)矩陣方程的降階,更適合后續(xù)直接法求解.本文需要用到主要特征基函數(shù)和次要特征基函數(shù)作為稀疏基,得到稀疏基后就可以準(zhǔn)確地恢復(fù)出目標(biāo)表面電流.

特征基函數(shù)法首先將目標(biāo)劃分成m個(gè)子域,如圖1所示.為了保證人工劃分子域邊界電流的連續(xù)性,需要對每個(gè)子域進(jìn)行擴(kuò)展.圖中虛線表示的是擴(kuò)展子域邊界,實(shí)線表示的是原始子域邊界,子域擴(kuò)展寬度取為λ/10—λ/8.

圖1 目標(biāo)子域劃分示例Fig.1.Example of target subdomain partition.

當(dāng)求解主要特征基函數(shù)(PCBFs)時(shí),不考慮子域與子域之間的互耦作用,將每一個(gè)擴(kuò)展子域看作獨(dú)立子域來處理,分別計(jì)算其在對應(yīng)入射波激勵下的響應(yīng),于是有

式中Vi是擴(kuò)展子域的對應(yīng)入射波激勵;是擴(kuò)展子域i的自阻抗矩陣,維數(shù)是Ni×Ni,Ni是擴(kuò)展子域i包含的RWG基函數(shù)個(gè)數(shù).擴(kuò)展子域PCBFs去除其擴(kuò)展部分就可以得到原始子域的PCBFs.

當(dāng)求解擴(kuò)展子域的次要特征基函數(shù)(SCBFs)時(shí),需要分別考慮擴(kuò)展子域i與其他m?1個(gè)原始子域之間的互耦作用,即

當(dāng)擴(kuò)展子域i與原始子域j之間不重合時(shí),阻抗矩陣Zij是擴(kuò)展子域i與原始子域j之間的互阻抗,維數(shù)是Ni×Nj;當(dāng)擴(kuò)展子域i與原始子域j之間有重合時(shí),阻抗矩陣Zij是擴(kuò)展子域i與去掉重合部分的原始子域j之間的互阻抗,維數(shù)是Ni×(Nj?Nij),Nij是重合部分的RWG基函數(shù)個(gè)數(shù).擴(kuò)展子域SCBFs去除擴(kuò)展部分即可得其原始子域SCBFs.

在求得CBFs后,目標(biāo)表面電流可以表示為

式中,是第i個(gè)子域的特征基函數(shù),是此特征基函數(shù)的加權(quán)系數(shù).用矩陣可表示為

簡寫為

即

Jr就是壓縮感知中所需要的稀疏基,a是目標(biāo)表面電流在稀疏基上的稀疏表示.在求得稀疏基后,隨機(jī)抽取M行阻抗矩陣ZM×N和其對應(yīng)的激勵VM×1,就可以用恢復(fù)算法準(zhǔn)確地恢復(fù)出目標(biāo)表面電流了.

需要指出的是該方法在劃分子域時(shí),要綜合考慮特征基函數(shù)和恢復(fù)算法的復(fù)雜程度.如果目標(biāo)劃分的子域數(shù)量m過多,雖然每個(gè)子域所含未知數(shù)很少,便于求解PCBFs以及SCBFs,但是壓縮感知稀疏度K卻會增大,給恢復(fù)算法帶來困難;如果目標(biāo)劃分子域的數(shù)量m過少,則每個(gè)子域所含未知數(shù)就會增多,CBFs生成效率會降低.在分析過程中應(yīng)結(jié)合實(shí)際需求平衡m與計(jì)算效率之間的關(guān)系.在本文的算例中,壓縮感知中的N為未知量數(shù)目,M取為N/2,K則根據(jù)文獻(xiàn)[17]選擇恢復(fù)概率大于90%的取值.m的取值要根據(jù)目標(biāo)的大小、形狀等綜合考慮,取在附近.

3 數(shù)值算例與結(jié)果

為了驗(yàn)證本文方法的有效性,分別對導(dǎo)體球、組合導(dǎo)體目標(biāo)和導(dǎo)彈模型的RCS進(jìn)行了計(jì)算.所有算例均在Intel(R)Core(TM)i7-3820 CPU@3.60 GHz,56 GB RAM的PC機(jī)上完成,算例均采用雙精度浮點(diǎn)運(yùn)算.

算例1計(jì)算了導(dǎo)體球的雙站RCS,導(dǎo)體球的半徑為1.5 m,入射波頻率為300 MHz,用三角單元剖分,共10970個(gè)三角形,未知量數(shù)目為16455,子域劃分為27個(gè),相鄰子域擴(kuò)展大小為1.2倍公共邊平均長度,擴(kuò)展后未知量數(shù)目為23253.圖2給出了三種方法得到的雙站RCS,由圖可見壓縮感知方法與中位數(shù)倍數(shù)(MOM)計(jì)算結(jié)果符合良好.三種算法的電流計(jì)算時(shí)間和誤差列于表1,可見CS比CBFM計(jì)算時(shí)間減少,且精度更高.本文用相對均方根誤差函數(shù)表示恢復(fù)誤差,即

圖2 導(dǎo)體球雙站RCSFig.2.Bistatic RCS of conducting sphere.

表1 導(dǎo)體球電流計(jì)算時(shí)間Table 1.Comparison of computation time.

算例2討論了三種方法在不同未知量數(shù)目時(shí)的雙站RCS的計(jì)算時(shí)間,本算例的目標(biāo)為半徑1 m的導(dǎo)體球,入射波頻率和剖分的未知量數(shù)目如表2所列,圖3給出了三種方法的計(jì)算時(shí)間與未知量數(shù)目的關(guān)系.由圖3可見,隨著未知量數(shù)目的增加,壓縮感知方法所用時(shí)間增長得更緩慢.

圖3 導(dǎo)體球電流計(jì)算時(shí)間Fig.3.Computation time of current of conducting sphere.

算例3計(jì)算了九個(gè)組合導(dǎo)體目標(biāo)的雙站RCS,其中導(dǎo)體球的半徑為0.5 m,導(dǎo)體圓錐半徑為0.5 m,高為1 m,導(dǎo)體圓柱的半徑為0.5 m,高為1 m,正方體邊長為1 m,各導(dǎo)體底面中心相距1.5 m.入射波頻率為300 MHz,用三角單元剖分,共16840個(gè)三角形,未知量數(shù)目為25260,子域劃分為50個(gè),相鄰子域擴(kuò)展大小為1.2倍公共邊平均長度,擴(kuò)展后未知量數(shù)目為40956.圖4給出了三種方法得到的雙站RCS,三種算法的計(jì)算時(shí)間和恢復(fù)誤差如表3所列,可見CS比CBFM計(jì)算時(shí)間減少,且精度明顯提高.

表2 不同入射波頻率下的未知量數(shù)目Table 2.Number of unknown variable of different frequency.

圖4 組合導(dǎo)體目標(biāo)雙站RCSFig.4.Bistatic RCS of composite conductor target.

表3 組合導(dǎo)體目標(biāo)電流計(jì)算時(shí)間比較Table 3.Comparison of computation time.

算例4 計(jì)算了導(dǎo)彈模型的雙站RCS,其中導(dǎo)彈模型的長為1 m,寬為0.64 m,高為0.22 m.入射波頻率為1 GHz,用三角單元剖分,共10502個(gè)三角形,未知量數(shù)目為15753,子域劃分為28個(gè),相鄰子域擴(kuò)展大小為1.2倍公共邊平均長度,擴(kuò)展后未知量數(shù)目為24078.圖5給出了三種方法得到的雙站RCS,三種算法的計(jì)算時(shí)間和恢復(fù)誤差如表4所列.

圖5 導(dǎo)彈模型雙站RCSFig.5.Bistatic RCS of simple missile model.

表4 導(dǎo)彈模型電流計(jì)算時(shí)間比較Table 4.Comparison of computation time.

4 結(jié) 論

本文改進(jìn)了壓縮感知和矩量法結(jié)合構(gòu)造欠定方程的方法,引入特征基函數(shù)作為稀疏基,恢復(fù)算法采用廣義正交匹配追蹤算法.應(yīng)用改進(jìn)方法計(jì)算三維導(dǎo)體目標(biāo)的RCS,數(shù)值計(jì)算結(jié)果表明,本文方法在保證精度的同時(shí),計(jì)算效率也得以提高.

本文雖然提高了目標(biāo)RCS的計(jì)算效率,但是在計(jì)算特征基函數(shù)時(shí),仍然需要填充稠密的阻抗矩陣.所以,快速獲得特征基函數(shù)(即稀疏基)可以進(jìn)一步提高計(jì)算效率.壓縮感知理論的新進(jìn)展——結(jié)構(gòu)化壓縮感知[19,20],可以將與數(shù)據(jù)采集硬件及復(fù)雜信號模型相匹配的先驗(yàn)信息引入傳統(tǒng)壓縮感知,從而實(shí)現(xiàn)對更廣泛類型的信號準(zhǔn)確有效的重建.我們今后會進(jìn)一步研究,將結(jié)構(gòu)化壓縮感知應(yīng)用到本文工作中.

[1]Gibson W C 2014J.Electromagn.Waves Appl.1 181

[2]Andriulli F P,Cools K,Bagci H,Olyslager F,Buffa A,Christiansen S,Michelssen E 2012IEEE Trans.Antennas Propag.56 2398

[3]Chen Y,Zuo S,Zhang Y,Zhao X,Zhang H 2017IEEE Trans.Antennas Propag.65 3782

[4]Candè E J,Wakin M B 2008IEEE Signal Process.Mag.25 21

[5]Ji S,Xue Y,Carin L 2008IEEE Trans.Signal Process.65 3782

[6]Ender J H G 2010IEEE Trans.Signal Process.65 3782

[7]Wang Z,Wang B Z 2014Acta Phys.Sci.63 120202(in Chinese)[王哲,王秉中 2014物理學(xué)報(bào) 63 120202]

[8]Chai S R,Guo L X,Li J,Li K 2015Asia-Pacific Microwave ConferenceNanjing,China Dec.6–9,2015 p1

[9]Kong M,Chen M S,Wu B,Wu X 2017IEEE Antennas Wirel.Propag.Lett.1 99

[10]Wang Z,Wang B Z,Wen Y Q,Wang R 2015IEEE International Symposium on Antennas and Propagation&Usnc/ursi National Radio Science MeetingVancouver,BC,July 19–24,2015 p1488

[11]Chao X Y,Chen M S,Wu X L,Shen J 2013Chin.J.Electron.41 2361(in Chinese)[曹欣遠(yuǎn),陳明生,吳先良,沈晶2013電子學(xué)報(bào)41 2361]

[12]Du H M,Chen M S,Wu X L 2012International Conference on Microwave and Millimeter Wave TechnologyShenzhen,China May 5–8,2012 p1

[13]Chai S R,Guo L X 2015Acta Phys.Sin.64 060301(in Chinese)[柴水榮,郭立新 2015物理學(xué)報(bào) 64 060301]

[14]Prakash V V S,Mittra R 2003Microw.Opt.Technol.Lett.36 95

[15]Sun Y F,Chan C H,Mittra R,Tsang L 2003Antennas and Propagation Society International SymposiumColumbus,OH,USA June 22–27,2003 p1068

[16]Wang Z G 2014Ph.D.Dissertation(Hefei:Anhui University)(in Chinese)[王仲根 2014博士學(xué)位論文(合肥:安徽大學(xué))]

[17]Jian W,Kwon S,Shim B 2012IEEE Trans.Signal Process.60 6202

[18]Tropp J A,Gilbert A C 2007IEEE Trans.Inf.Theory53 4665

[19]Baraniuk R G,Cevher V,Duarte M F,Hegde C 2010IEEE Trans.Inf.Theory56 1982

[20]Duarte M F,Eldar Y C 2011IEEE Trans.Signal Process.59 4053