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

付永慶，劉 偉

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖1 陣元間脈沖響應(yīng)Fig.1 Interelement impulse response

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

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

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

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

用矩陣形式表示：

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

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

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

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

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

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

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

因?yàn)镵（）ω 具有對稱性，所以有：

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

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

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

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

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

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

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

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

寫成矩陣形式有：

進(jìn)而可得到：

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

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

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

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

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

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

（4）確定目標(biāo)搜索域及各觀測點(diǎn)。

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

3 仿真分析

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

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

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

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

為分析基于時間反轉(zhuǎn)算子分解的選擇性聚焦算法在不同信噪比下的聚焦性能，與圖3實(shí)驗(yàn)相同，取C1∶C2∶C3＝1∶0.5∶0.3，3個主特征向量V1、V2、V3分別對應(yīng)目標(biāo)1、2、3，并在仿真時加入信噪比不同的高斯白噪聲。利用不同主特征向量實(shí)現(xiàn)其對應(yīng)目標(biāo)處的選擇性聚焦，并執(zhí)行100次蒙特卡洛實(shí)驗(yàn)，圖4為不同信噪比情況下基于時間反轉(zhuǎn)算子分解的選擇性聚焦算法的均方根誤差。從圖中可以看出：隨著信噪比的增加，該算法選擇性聚焦的均方根誤差逐漸變小，而且在信噪比為0dB的情況下，不同目標(biāo)的選擇性聚焦誤差都在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條件下的仿真結(jié)果，如圖5和圖6所示。從圖中可以看出：改變目標(biāo)的反射系數(shù)，時間反轉(zhuǎn)算子分解算法依然可以很好地實(shí)現(xiàn)在各不同目標(biāo)處的選擇性聚焦，只是特征向量V1、V2、V3所對應(yīng)的目標(biāo)位置發(fā)生了改變。綜合圖3、圖5和圖6可以看出：非零特征值的大小與對應(yīng)目標(biāo)的信號反射率有關(guān)，目標(biāo)反射率越大，對應(yīng)的特征值也越大，因此可以根據(jù)特征值的大小選擇相應(yīng)的特征向量實(shí)現(xiàn)在不同目標(biāo)處的聚焦。

圖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 結(jié)束語

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

［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］.沈陽建筑大學(xué)學(xué)報(bào)，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］.應(yīng)用聲學(xué)，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］.哈爾濱工程大學(xué)學(xué)報(bào)，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ì)量經(jīng)濟(jì)學(xué)模型［J］.交通運(yùn)輸工程學(xué)報(bào)，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的高速公路事故預(yù)測模型［J］.北京工業(yè)大學(xué)學(xué)報(bào)，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.