一種基于壓縮感知的三維導體目標電磁散射問題的快速求解方法

丁亞輝 孫玉發 朱金玉

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

1 引 言

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

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

2 基本原理

2.1 壓縮感知構造欠定方程

為了提高矩量法的計算效率,引入信號處理領域中壓縮感知理論來降低內存需求,減少電流求解時間.

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

根據矩量法,可得到矩陣方程

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

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

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

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

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

2.2 特征基函數作為稀疏基

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

特征基函數法首先將目標劃分成m個子域,如圖1所示.為了保證人工劃分子域邊界電流的連續性,需要對每個子域進行擴展.圖中虛線表示的是擴展子域邊界,實線表示的是原始子域邊界,子域擴展寬度取為λ/10—λ/8.

圖1 目標子域劃分示例Fig.1.Example of target subdomain partition.

當求解主要特征基函數(PCBFs)時,不考慮子域與子域之間的互耦作用,將每一個擴展子域看作獨立子域來處理,分別計算其在對應入射波激勵下的響應,于是有

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

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

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

在求得CBFs后,目標表面電流可以表示為

式中,是第i個子域的特征基函數,是此特征基函數的加權系數.用矩陣可表示為

簡寫為

即

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

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

3 數值算例與結果

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

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

圖2 導體球雙站RCSFig.2.Bistatic RCS of conducting sphere.

表1 導體球電流計算時間Table 1.Comparison of computation time.

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

圖3 導體球電流計算時間Fig.3.Computation time of current of conducting sphere.

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

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

圖4 組合導體目標雙站RCSFig.4.Bistatic RCS of composite conductor target.

表3 組合導體目標電流計算時間比較Table 3.Comparison of computation time.

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

圖5 導彈模型雙站RCSFig.5.Bistatic RCS of simple missile model.

表4 導彈模型電流計算時間比較Table 4.Comparison of computation time.

4 結 論

本文改進了壓縮感知和矩量法結合構造欠定方程的方法,引入特征基函數作為稀疏基,恢復算法采用廣義正交匹配追蹤算法.應用改進方法計算三維導體目標的RCS,數值計算結果表明,本文方法在保證精度的同時,計算效率也得以提高.

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

[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物理學報 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)[曹欣遠,陳明生,吳先良,沈晶2013電子學報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物理學報 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博士學位論文(合肥:安徽大學)]

[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