基于時間反轉(zhuǎn)算子分解的選擇性聚焦方法

付永慶，劉 偉

（哈爾濱工程大學 信息與通信工程學院，哈爾濱150001）

時間反轉(zhuǎn)理論由光學中的相位共軛產(chǎn)生，F(xiàn)ink等［1］于1989年驗證了時間反轉(zhuǎn)的空間聚焦特性并首次給出了時間反轉(zhuǎn)鏡（Time reversal mirror，TRM）的概念。隨后，國內(nèi)外專家在水聲領域做了大量的相關研究及實驗并取得了一定的研究成果［2－3］。到目前為止，時間反轉(zhuǎn)理論已經(jīng)在超聲碎石［4］、無損探傷［5］、水下目標探測［6］、無線通信［7－8］和 成像［9－10］等 領 域 得 到 廣 泛 探 討 與 研 究。時間反轉(zhuǎn)技術最主要的優(yōu)點是不需要介質(zhì)性質(zhì)和陣列分布等先驗知識就可以實現(xiàn)自適應聚焦［11－12］，因此越來越受到人們的重視。對于信號反射強度不同的多個目標，為了實現(xiàn)在目標上的選擇性聚焦，可通過有限次數(shù)的時間反轉(zhuǎn)迭代逐漸削弱反射系數(shù)較小的目標處的相對聚焦信號能量，最終實現(xiàn)在強反射目標上的選擇性聚焦［13］，但通過迭代時間反轉(zhuǎn)使信號聚焦于弱反射目標處是不可能實現(xiàn)的。

為了實現(xiàn)在信號反射強度不同的多目標處的選擇性聚焦，本文提出了一種基于電磁平面波波動方程的時間反轉(zhuǎn)算子分解算法，仿真實驗顯示，根據(jù)時間反轉(zhuǎn)算子分解的特征值和特征向量可以有效地實現(xiàn)在反射系數(shù)不同的目標處的選擇性聚焦。

1 波動方程的時間反轉(zhuǎn)不變性

時間反轉(zhuǎn)理論以波動方程的時間反轉(zhuǎn)不變性為基礎，把天線陣列接收到的目標信號進行時間反轉(zhuǎn)后重新發(fā)射到空間中，則時間反轉(zhuǎn)后信號經(jīng)過反向傳輸并實現(xiàn)在目標處的重聚焦。

假設電磁場中的標量位函數(shù)為φ（r，t），則在無源場區(qū)域滿足的標量波動方程為：

式中：μ 和ε 分別為介質(zhì)的介電常數(shù)和磁導率；r為電磁波傳播的位移矢量。

考慮平面電磁波，則式（1）所示標量波動方程的通解為：

式中：f1和f2為任意函數(shù)為電磁波傳播速度。

可以看出，式（3）和式（4）分別表示傳播方向相反的兩類波形，且它們都是波動方程的解。

對式（3）和式（4）做時間反轉(zhuǎn)處理，可得：

2 時間反轉(zhuǎn)算子的分解算法

2.1 時間反轉(zhuǎn)算子

考慮一個由N 個天線陣元組成的TRM 陣列，有N×N 個交叉陣元脈沖響應，如圖1所示。令klm（）t 表示當一個時域delta函數(shù)作用于陣元m 時，陣元l接收到的信號。根據(jù)第1節(jié)對電磁波波動方程的描述可將平面電磁波在空間中的傳輸函數(shù)寫成：

式中：τ為電磁波的傳播時延。

圖1 陣元間脈沖響應Fig.1 Interelement impulse response

從陣元m 到陣元l之間的交叉脈沖響應可以表示為：

式中：τm、τl分別為目標到第m 個和第l個陣元的傳播時延；＊表示卷積。

令em（）t（1≤m≤N），是天線陣列的發(fā)射信號，信號經(jīng)過目標的反射后被天線陣列重新接收，則接收信號rl（）t（1≤l≤N），可以表示為：

式（9）寫成頻域形式為：

用矩陣形式表示：

式中：E（）ω 和R（）ω 分別表示N 個天線陣元的發(fā)射和接收信號向量；K（）ω 為傳輸矩陣，維數(shù)為陣元個數(shù)。

根據(jù)互易定理，從陣元m 到陣元l 之間的交叉脈沖響應與從陣元l到陣元m 之間的交叉脈沖響應相同，因此傳輸矩陣K（）ω 是對稱陣。

設E0（）ω 為天線陣列第一次發(fā)射的信號向量，則經(jīng)過目標反射后天線陣列的接收信號可用傳輸矩陣表示為：

時域內(nèi)的時間反轉(zhuǎn)操作等同于頻域內(nèi)的相位共軛，則天線陣列第二次發(fā)射的信號向量為第一次接收信號向量的相位共軛：

經(jīng)過兩次時間反轉(zhuǎn)處理后，天線陣列的發(fā)射信號向量可以表示為：

以此類推，經(jīng)過2n 次和2n＋1次時間反轉(zhuǎn)處理后，天線陣列的發(fā)射信號向量可以分別表示為：

式中：K＊（）ω K（）ω 稱為時間反轉(zhuǎn)算子。

因為K（）ω 具有對稱性，所以有：

式中：T 表示對矩陣的轉(zhuǎn)置運算。

式（17）說明K＊（）ω K（）ω 是埃爾米特矩陣并可以進行對角化，其特征向量是正交的，特征值為實數(shù)。將時間反轉(zhuǎn)算子進行特征值分解可得到與目標數(shù)目相同的非零主特征值，其對應的主特征向量與目標存在著一一對應的關系。令λ為時間反轉(zhuǎn)算子的特征值，v（）ω 為相應的特征向量，則由式（18）可知λ為非負值。

2.2 時間反轉(zhuǎn)算子的分解算法

假設目標域中共有d 個目標，各目標對信號的反射率分別為C1，C2，…，Cd。第i（1≤i≤d）個目標與接收陣列第n（1≤n≤N）個陣元之間的傳輸函數(shù)為hin（）t，其傅里葉變換為Hin（）ω，則第i個目標與接收天線陣之間的傳輸向量為：

如果點目標理想可分辨，即在其中一個目標上的時間反轉(zhuǎn)聚焦不會給其他目標提供能量，則Hi（）ω 正交，且傳輸矩陣可以寫為如下形式［13］：

從式（20）可以看出傳輸矩陣由3個矩陣組成：描述從陣列到目標的前向傳輸矩陣，目標反射系數(shù)組成的反射矩陣和描述從目標到陣列的反向傳輸矩陣。

則S（）ω 的第n 個元素可以表示為：

由Hi（）ω 的正交性可得：

寫成矩陣形式有：

進而可得到：

可見，時間反轉(zhuǎn)算子對應某目標的特征向量是該目標到天線陣列之間的傳輸向量的共軛，那么該特征向量與目標區(qū)域中的觀測點對應的傳輸函數(shù)的內(nèi)積在其對應的目標處將趨于最大，可實現(xiàn)在該目標處的信號能量聚焦。假設目標域中共有M 個目標觀測點，于是可得目標域中第j（1≤j≤M）個觀測點處時間反轉(zhuǎn)算子分解的選擇性聚焦目標函數(shù)公式為：

式中：Hjn（）ω 為天線陣列第n 個陣元到目標域中第j個觀測點的傳輸函數(shù)；vi（）n 為時間反轉(zhuǎn)算子分解后第i個非零特征值所對應特征向量中的第n 個元素。

根據(jù)以上分析，可將基于時間反轉(zhuǎn)算子分解的選擇性聚焦算法分為以下幾個步驟：

（1）利用天線陣列的第一個陣元發(fā)射中心頻率為ω0的目標探測信號，經(jīng)目標反射后被天線陣列接收。對接收信號進行傅里葉變換可得第一個陣元與所有陣元之間的交叉脈沖響應。

（2）用同樣的信號對天線陣列的其他N－1個陣元進行激勵，重復步驟（1）的操作可獲得傳輸矩陣K（）ω 。

（3）計算時間反轉(zhuǎn)算子K＊（）ω K（）ω，并利用式（27）將其進行特征值分解獲得主特征值與主特征向量，其中主特征值的個數(shù)即為目標個數(shù)。

（4）確定目標搜索域及各觀測點。

（5）根據(jù)式（26）并利用目標所對應的主特征向量計算目標域中各觀測點的目標函數(shù)值，可取得在該目標處的選擇性聚焦。

3 仿真分析

為了驗證基于時間反轉(zhuǎn)算子分解的選擇性聚焦特性，建立如圖2所示的均勻圓形天線陣列模型，其中N 個陣元均勻分布在以r 為半徑的圓周上，且陣列的幾何中心與坐標原點重合。仿真實驗中取N ＝32，r為20m，天線陣列發(fā)射中心頻率為30 MHz的目標探測波。考慮以坐標原點為中心，半徑為20km 的圓周區(qū)域為待偵測目標區(qū)域，且在整個圓周域上平均分布著360個觀測點，偵測步長為1°，目標位于待偵測區(qū)域內(nèi)，且方位角分別為22°，45°，70°。

圖2 均勻圓形TRMFig.2 Uniform circular TRM

為實現(xiàn)時間反轉(zhuǎn)算子分解的選擇性聚焦，首先通過時間反轉(zhuǎn)陣列的發(fā)射和接收信號獲得傳輸矩陣，并計算得到時間反轉(zhuǎn)算子，然后對時間反轉(zhuǎn)算子進行特征值分解，其中主特征值的個數(shù)為目標個數(shù)，而主特征值對應的主特征向量則包含了對應目標的方位信息。利用主特征向量并根據(jù)式（26）計算各觀測點的目標函數(shù)值，其中函數(shù)值最大點所對應的觀測點的方位即為該主特征值所對應目標的方位，因此利用每個主特征值進行目標域中各觀測點的目標函數(shù)值計算，可分別實現(xiàn)在不同目標處的選擇性聚焦。

首先取C1∶C2∶C3＝1∶0.5∶0.3，通過仿真計算獲得陣元間傳輸矩陣K（）ω，則時間反轉(zhuǎn)算子可以寫成K＊（）ω K（）ω ，對其進行特征值分解后，歸一化特征值分布如圖3（a）所示。從圖中可以看出：由于目標域中只有3個目標，所以時間反轉(zhuǎn)算子分解后只有3個非零特征值。令3個大小依次排列的特征值對應的特征向量分別記為V1、V2和V3，并將這3個特征向量作為權向量分別計算目標域中的目標函數(shù)，可得如圖3（b）所示目標域中的歸一化目標函數(shù)信號強度分布情況。從圖中可以看出，利用弱反射目標所對應的特征向量可以實現(xiàn)在弱目標處的有效聚焦，而且利用每個主特征向量可以分別實現(xiàn)在相應目標處的準確聚焦，聚焦方位分別為22°、45°、70°，這說明每個主特征向量分別為天線陣列提供了相應目標的方位信息，可以有效地聚焦在相應的目標處，很好地實現(xiàn)了選擇性聚焦。

為分析基于時間反轉(zhuǎn)算子分解的選擇性聚焦算法在不同信噪比下的聚焦性能，與圖3實驗相同，取C1∶C2∶C3＝1∶0.5∶0.3，3個主特征向量V1、V2、V3分別對應目標1、2、3，并在仿真時加入信噪比不同的高斯白噪聲。利用不同主特征向量實現(xiàn)其對應目標處的選擇性聚焦，并執(zhí)行100次蒙特卡洛實驗，圖4為不同信噪比情況下基于時間反轉(zhuǎn)算子分解的選擇性聚焦算法的均方根誤差。從圖中可以看出：隨著信噪比的增加，該算法選擇性聚焦的均方根誤差逐漸變小，而且在信噪比為0dB的情況下，不同目標的選擇性聚焦誤差都在1°以內(nèi)，顯示了該算法良好的聚焦性能。

圖3 C1∶C2∶C3＝1∶0.5∶0.3時，時間反轉(zhuǎn)算子分解的選擇性聚焦Fig.3 Selective focusing of decomposition of time reversal operator in C1∶C2∶C3＝1∶0.5∶0.3

圖4 不同信噪比下選擇性聚焦的均方根誤差Fig.4 RMSE of selective focusing in different SNR

利用與圖3同樣的方法可以獲得C1∶C2∶C3＝0.5∶1∶0.3和C1∶C2∶C3＝0.3∶0.5∶1條件下的仿真結果，如圖5和圖6所示。從圖中可以看出：改變目標的反射系數(shù)，時間反轉(zhuǎn)算子分解算法依然可以很好地實現(xiàn)在各不同目標處的選擇性聚焦，只是特征向量V1、V2、V3所對應的目標位置發(fā)生了改變。綜合圖3、圖5和圖6可以看出：非零特征值的大小與對應目標的信號反射率有關，目標反射率越大，對應的特征值也越大，因此可以根據(jù)特征值的大小選擇相應的特征向量實現(xiàn)在不同目標處的聚焦。

圖5 C1∶C2∶C3＝0.5∶1∶0.3時，時間反轉(zhuǎn)算子分解的選擇性聚焦Fig.5 Selective focusing of decomposition of time reversal operator in C1∶C2∶C3＝0.5∶1∶0.3

圖6 C1∶C2∶C3＝0.3∶0.5∶1時，時間反轉(zhuǎn)算子分解的選擇性聚焦Fig.6 Selective focusing of decomposition of time reversal operator in C1∶C2∶C3＝0.3∶0.5∶1

4 結束語

提出了一種基于平面電磁波波動方程時間反轉(zhuǎn)不變性的時間反轉(zhuǎn)算子分解選擇性聚焦方法，該方法通過對時間反轉(zhuǎn)算子的分解可以獲得與目標數(shù)目一致的主特征值和與目標方位一一對應的主特征向量。根據(jù)每個目標所對應的特征向量計算目標域中所有觀測點處的目標函數(shù)值可以實現(xiàn)在各目標處的選擇性聚焦。仿真實驗結果顯示該方法具有有效性，利用弱反射目標所對應的特征向量可實現(xiàn)在弱目標處的聚焦，利用所有的主特征向量可以實現(xiàn)在不同目標上的選擇性聚焦，另外，主特征值的大小和各目標所對應的特征向量次序都與目標的反射系數(shù)有關。該方法可以用于對多個目標中的弱目標進行檢測或者是對多目標進行選擇性檢測的情況，具有很好的應用前景。

［1］Fink M，Prada C，Wu F，et al.Self focusing in inhomogeneous media with time reversal acoustic mirrors［C］∥IEEE Ultrasonics Symposium Proceedings，Montreal，Que：IEEE，1989：681－686.

［2］Song H C，Kuperman W A，Hodgkiss W S，et al.Iterative time reversal in the ocean［J］.Journal of the Acoustical Society of America，1999，105（6）：3176－3184.

［3］Sheng X L，Bao X Z，Hui J Y，et al.Underwater Passive Location Technology Using Dummy Timereversal Mirror［M］.Lecture Notes in Computer Science，Berlin：Springer，2008：605－612.

［4］Thomas J L，F(xiàn)ink M.Self focusing on extended objects with time reversal mirror，application to lithotripsy［C］∥IEEE Ultrasonics Symposium Proceedings，Cannes，F(xiàn)rance：IEEE，1994：1809－1814.

［5］趙乃至，閻石，齊霽.基于時間反轉(zhuǎn)法的管道周向裂縫損傷檢測［J］.沈陽建筑大學學報，2013，29（1）：44－49.Zhao Nai－zhi，Yan Shi，Qi Ji.The pipeline circumferential cracks damage detection based on time reversal method［J］.Journal of Shenyang Jianzhu University，2013，29（1）：44－49.

［6］李壯，喬鋼，王健培，等.基于虛擬時間反轉(zhuǎn)鏡的短基線定位研究［J］.應用聲學，2012，31（4）：256－261.Li Zhuang，Qiao Gang，Wang Jian－pei，et al.Short baseline positioning based on virtual time reversal mirror［J］.Applied Acoustics，2012，31（4）：256－261.

［7］Khan W N，Zubair M，Wyne S，et al.Performance evaluation of time－reversal on measured 60 GHz wireless channels［J］.Wireless Personal Communications，2013，71（1）：701－707.

［8］Chen Y，Yang Y H，Han F，et al.Time－reversal wideband communications［J］.IEEE Signal Processing Letters，2013，20（12）：1219－1222.

［9］Choi H，Ogawa Y T，Ohgane T.Time－reversal MUSIC imaging with time－domain gating technique［J］.IEICE Transactions on Communications，2012，E95－B（7）：2377－2385.

［10］Liu X F，Wang B Z，Li J L.Transmitting－mode time reversal imaging using MUSIC algorithm for surveillance in wireless sensor network［J］.IEEE Transactions on Antennas and Propagation，2012，60（1）：220－230.

［11］Fu Y Q，Jiang Y L，Liu Z Y.Near－field source localization method and application using the time reversal mirror technique［J］.Journal of Electronics（China），2011，28（4）：531－538.

［12］夏云龍，付永慶.時間反轉(zhuǎn)算子特征值分解算法［J］.哈爾濱工程大學學報，2008，29（12）：1340－1343.Xia Yun－long，F(xiàn)u Yong－qing.Eigenvalue decomposition algorithm for time reversal operators［J］.Jousnal of Harbin Engineering University，2008，29（12）：1340－1343.

［13］Prada C，Mannecille A，Spoliansky D，et al.Decomposition of the time reversal operator：Detection and selective focusing on two scatters［J］.Journal of the Acoustical Society of America，1996，99（4）：2067－2076.

［14］Qin X，Reyes P.Conditional quantile analysis for crash count data［J］.Journal of Transportation Engineering，2011，137（9）：601－607.

［15］王軍雷，孫小端，賀玉龍，等.交通事故宏觀計量經(jīng)濟學模型［J］.交通運輸工程學報，2012，12（2）：70－75.Wang Jun－lei，Sun Xiao－duan，He Yu－long，et al.Macroscopic econometrics model of traffic accident［J］.Journal of Traffic and Transportation Engineering，2012，12（2）：70－75.

［16］徐婷，孫小端，王偉力，等.基于Panel Data的高速公路事故預測模型［J］.北京工業(yè)大學學報，2010，36（4）：495－499.Xu Ting，Sun Xiao－duan，Wang wei－li，et al.Highway accidents statistical analysis with panel data model［J］.Journal of Beijing University of Technology，2010，36（4）：495－499.

［17］Kweon Y J，Kockelman K M.Spatially disaggregate panel models of crash and injury counts：the effect of speed limits and design［R］.Transportation Research Board 83rd Annual Meeting，2004.

［18］Openshaw S.The modifiable areal unit problem［J］.Concepts and Techniques in Modern Geography，1984.

［19］Hausman J，Hall B，Griliches Z.Econometric models for count data with an application to the patents－R＆D relationship［J］.Econometrica，1984，52（4）：909－938.

［20］Oh J，Lyon C，Washington S P，et al.Validation of the FHWA crash models for rural intersections：lessons learned［J］.Transportation Research Record，2003，1840：41－49.