高效的MIMO雷達運動目標三維成像方法

王偉，胡子英，岳佳男

（哈爾濱工程大學自動化學院，黑龍江 哈爾濱 150001）

1 引言

多輸入多輸出（MIMO, multiple-input multiple-output）雷達是近年來提出的一類新型雷達。區別于其他工作方式，MIMO雷達在發射端發射正交信號，在接收端通過匹配濾波實現信號的通道分離，這樣的方式使它擁有更多的虛擬孔徑和更高的陣列自由度，極大地提高了雷達成像的性能[1-2]。為了更好地呈現實際目標的三維結構及特征，近年來，三維成像技術開始在MIMO雷達中得到研究及發展，現階段主要從優化陣元的排布方式[3]、提高成像精度和效率[4-5]等方面展開。但是這些研究主要圍繞空間中的靜態目標，在實際應用中，目標運動因多普勒頻移帶來的影響不可忽視。

首先，多普勒頻移的產生會使回波信號的中心頻率發生改變，導致估計結果與實際情況存在較大偏差。其次，對多普勒頻率的采樣使測量矩陣各列間的互相干性不再滿足Welch界限制條件[6]，稀疏重構的恢復概率難以得到保證。文獻[7]從速度帶來的頻移影響考慮，用迭代自適應估計方法（IAA,iterative adaptive approach）對距離-角度-多普勒的三維特征參數進行估計，但IAA的高運算復雜度使其難以應用于三維成像中。文獻[8]在回波中考慮了速度帶來的多普勒影響，并用變分貝葉斯的方法對對速度參數進行精確的估計，但多層貝葉斯模型使對速度求解和稀疏成像的迭代過程更加復雜。

對于多普勒頻率采樣帶來的強互相干性問題，文獻 [9-11]分析了各列之間的互相干性[12]對測量矩陣性能的影響，并以Welch界判定準則構造凸優化問題，找到最優投影矩陣，使測量矩陣的性能逼近最佳。文獻[13]則在此基礎上利用信號的先驗信息簡化了該凸優化問題的處理過程，預先給定加權矩陣，提高了算法效率。但是在這些方法中，投影矩陣及原始測量矩陣均參與運算。使整個優化過程尤為復雜，并且MIMO雷達稀疏成像中需要將區域每點的散射系數張成一個向量來處理，使整個模型的數據結構變得極大[14]，因此，運動目標三維成像亟需解決運算效率過低的問題。針對提升算法運算效率，文獻[15]利用矩陣完備性（MC, matrix completion）理論構造具有低秩和均勻隨機采樣性質的矩陣，來代替測量矩陣實現三維成像，成像效果良好，但對運算效率提升的效果一般。文獻[16]中將所得回波模型用基于克羅內克壓縮感知（KCS,Kronecker compressive sensing）的方式表示，并在算法的初始化階段利用匹配濾波得到的結果提取目標大致范圍，借此來優化測量矩陣，以降低 KCS成像的運算復雜度。文獻[17]則在初始化目標區域的基礎上實現了三維成像的張量分解，利用此結論為Tucker分解提供了先驗，在保證其唯一性條件的前提下提升了運算效率，但文獻[16-17]中的方法普遍受噪聲的影響比較大，不適合在偏低的信噪比下工作。

本文針對 MIMO雷達遠場條件下的運動目標提出一種高效的三維成像方法。以多普勒頻率作為一個新維度，在進行對目標區域的稀疏反演前，依據各維度的回波依次完成對目標的搜索，得到空間目標在各個維度的分布信息，并根據索引重構一個新的低維測量矩陣，縮小目標區域范圍，有效地提高了運算效率。同時分析多普勒頻率采樣對測量矩陣相干性的影響，并預先構造一個投影加權矩陣，應用高精度貝葉斯方法實現對該投影矩陣的優化，得到目標沿多普勒維度的精確分布信息。仿真結果表明，所提方法可以極大地提升運算效率，實現高效成像，具有精確穩定的成像性能。

2 MIMO雷達運動目標回波模型

MIMO雷達發射陣列是由M個陣元構成的均勻線性陣列，以Ti(i=1,…,M)表示第i個陣元，且在該發射陣元采用的是Hadamard正交相位編碼信號s(t) =p(t) ej2πfct+?i，其中，p(t)為信號的復包絡，iiifc為發射信號載波頻率，?i是信號的相位。此類編碼信號通過調整發射信號的相位?i來保證各發射信號間的正交性。接收陣列為由NR×NC個陣元構成的均勻分布二維陣列[16]，以Rj′l(j′ = 1 ,… ,NR；l=1,…,NC)表示接收陣元位于接收陣列的第j′行第l列。假設空間中的目標服從Sweiling II起伏模型，即對空間中K個遠場運動目標點，在一個脈沖周期內散射系數不變，并以目標對應的多普勒頻率fdk作為新的變量來表示目標運動帶來的影響[18-19]，則第i個發射信號經過空間中K個運動目標點的反射后，由Rj′l接收到的信號為

其中，σk為目標點k的散射系數，τk,i,j’,l為信號經由第k個散射點的時延，由目標及陣元位置決定。經過去載波處理后可得接收端回波信號，如式(2)所示。

其中，q=1,…,Q表示第q次采樣，tq為對應采樣時間。按照壓縮感知成像理論，對目標區域劃分網格點，沿目標的3個維度劃分的網格點數分別為U1、以表示網格點處的多普勒頻率，則式(2)的回波實際應為

根據式(3)的模型進行稀疏成像，需依據式中的旋轉矢量構造測量矩陣，就大多數研究而言，應用于CS成像方法中的測量矩陣是可以預先設計的，但式(3)中多普勒頻率的存在使測量矩陣中出現了未知項，無法實現對測量矩陣的精確構建，且矩陣中的不確定性因素會對稀疏目標的恢復產生嚴重的影響，因此需要在對目標的三維坐標進行估計的同時對目標的多普勒頻率進行精確的估計，故本文將目標的多普勒頻率作為X、Y、Z坐標軸外的另一變量進行處理，且對式(2)進一步處理，得

其中，Tik表示Ti發射陣元到第k個目標點的距離，Rj'lk表示第k個目標點到Rj'l接收陣元的距離。由于實際中目標速度遠小于信號傳播速度，即所以可近似處理因此有

在成像區域中選定點O作為參考中心，且令TiO和Rj’lO分別表示O點到發射陣列、接收陣元的距離，則由O點反射的回波信號可作為參考信號實現對目標回波信號的補償，即

如圖1所示，遠場目標區域中共有K個目標點，以O點為參考中心，在MIMO雷達收發陣列中，以T0和R00分別作為發射陣列、接收陣列的參考陣元。

圖1 MIMO雷達遠場成像系統幾何模型

在文獻[16,20]中對式(6)中所示的距離項做了等效近似處理，如式(7)所示。

綜上，可以構造運動目標稀疏回波模型為

其中，有

且對測量矩陣有

3 MIMO雷達運動目標的高效三維成像

傳統的 CS稀疏成像，需在各維度預先劃分網格點，并在此基礎上構成測量矩陣。通過稀疏恢復得到的稀疏解中包含了目標在區域內的分布信息，求解散射稀疏向量的稀疏恢復過程依賴測量矩陣每列與每個網格點的對應關系。對于稀疏成像而言，空間中大部分的位置散射強度認定為0，說明在稀疏恢復的過程中，測量矩陣的大部分列并不對應目標分布，卻仍然要進行算法運算，這增加了算法的運算復雜度。在式(12)～式(14)模型中，上述無關列的無用運算量問題更明顯，因此本文選擇對測量矩陣進行處理，從中僅選取稀疏相關列，將大部分的無關列剔除，那么稀疏求解的運算效率就會大大提升。同時，由于測量矩陣的互相干特性主要由多普勒維度的采樣所影響，基于本文方法可以單獨在多普勒維度實現對測量矩陣的優化，提升整體算法的優化效率。

3.1 目標區域優化搜索

首先要完成的是得到所需要的相關列的索引，這些索引對應空間中目標點的位置，因此可以利用回波從各個維度搜尋目標信息，得到在每個維度目標的粗略分布情況，以構建索引支撐集。那么各維度都可以得到一個典型的CS模型，如式(15)所示。

如此構成了各維度的稀疏恢復問題，以在X維的搜索為例，上述對應X維的模型可以轉化為對如式(18)所示的凸優化問題的求解。

各種方法（如正交匹配追蹤（OMP, orthogonal matching pursuit）算法、稀疏貝葉斯學習（SBL,sparse bayesian learning）方法等）均可以實現對以上凸優化問題的求解，得到目標區域沿X維的分布信息，在劃分的X維網格x1,x2,…,xUX中，可以得到對應目標在X維的位置索引構成的集合{ΛX|Λx1,Λx2,…,ΛxPx,Px≤K}，Px≤K是考慮到了空間中不同目標點可能會有相同的X維坐標，較易引起搜索過程中目標點的丟失，這Px個位置包含了空間中所有的K個目標點的X坐標，如此便完成利用接收回波在X維度上對目標的搜尋工作。

同理，經過類似的處理可以得到其余3個維度上的坐標索引集合分別為{ΛY|Λy1,Λy2,…,ΛyPy，Py≤K}、{ΛZ|Λz1,Λz2,…,ΛzPz，Pz≤K}和{Λf|Λf1,Λf2,…,ΛfPf，Pf≤K}，將所得的 4 個集合進行外積處理，可以得到Px×Py×Pz×Pf個組合指向成像區域內的Px×Py×Pz×Pf個點，這些目標點的集合Λ中必然包含了所有的真實目標點，即Px×Py×Pz×Pf≥K，于是可將此集合作為新支撐集的索引構造測量矩陣ψ∈ CMNRNC Q×(Px×Py×Pz×Pf)， 式(12)中的壓縮感知模型可轉換為如式(19)所示的簡化模型。

而P?U，故采用常用的稀疏恢復算法時能實現良好的提升運算效率的效果。此外，由于各維度的搜索過程是依靠稀疏恢復算法完成的，而大多數算法可以實現良好的抗噪性，因此本文所提方法在低信噪比下也能正常工作。但需要注意兩點，一是各維的CS模型需要保證存在唯一解；二是由式(8)可知，經過上述過程得到的稀疏解中，多普勒頻率是目標在成像區域真實的多普勒頻率，但其余坐標并不是目標在成像區域的真正坐標值，將其余坐標轉換為實際坐標還需進行相應處理。基于式(8)和式(9)，令則目標其余坐標的真實坐標和得到的坐標之間的關系如式(21)所示。

因此，根據式(21)可得目標其余坐標在成像區域的目標坐標值。

經過上述處理后，得到了目標的三維坐標及多普勒頻率，即實現了對運動目標三維成像算法運算效率的提升，但多普勒頻率的引入會對測量矩陣的性能造成影響，對此需要分析測量矩陣性能的變化，并選取合適的方法來消除此影響，改善測量矩陣的性能，獲得良好的稀疏恢復效果。

3.2 互相干性優化

本節對測量矩陣的性能進行優化，以實現良好的稀疏恢復。在常用的壓縮感知模型中，常用約束等距特性（RIP, restricted isometry property）[12]來衡量性能，但是其在實際中的應用比較困難，因此可以用運算較便捷的互相干性來代替 RIP作為判斷的準則。對一個M′×N′維的測量矩陣Θ，可以用Gram矩陣G=ΘHΘ來直觀描述各列之間的互相干性，在對測量矩陣Θ的每列進行歸一化處理后，最理想的情況下G應為單位矩陣，即只有對角線元素為1，其他位置元素均為0，代表各列向量之間是正交的。實際中這種完全理想的矩陣是極難得到的，因此能實現稀疏重構的較理想矩陣G[21]應該滿足對角線元素為1，并且對角線外其他位置的元素必須不超過Welch界[8]的限制，然而在式(19)的壓縮感知模型難以保證其性能。在X、Y、Z三維中，對目標信息的搜索沒有問題，而影響整體測量矩陣性能的主要因素出現在多普勒維度的搜索中，由式(16)可知，多普勒維度上Q個采樣點對應的采樣時間而遠場條件下采樣時間的選取存在一個弊端，如圖 2所示。

圖2 遠場采樣等效距離

圖2中，O為成像區域的中心，l1和l2表示成像區域與雷達間的最大距離和最小距離，R0表示成像中心到雷達的參考距離，Δl1、Δl2分別表示l1、l2與R0的距離差。信號發出后，在遠場條件下接收端最早到最晚接收到信號的時間差為

由此可見，采樣時間的取值與目標區域范圍大小及信號傳輸速度有關，但要保證探測的精度，目標范圍就不宜較大且實際中c認定為光速，這就導致Δt極小，而基于Q次采樣均是在該時間內完成的，所以各采樣時間取值均過小且在同一量級，使測量矩陣Af中的采樣數據比較相近，致其相關性較差，難以保證稀疏重構的可靠性。

其中，g表示格拉姆矩陣。

首先，Af先經過列歸一化處理，然后求解式(23)得到的Φ、g的精確聯合估計，使測量矩陣的性能最佳。在解決凸優化問題的常用算法中，Bayes方法依據先驗分布信息及后驗概率密度實現對參數的精確估計，可獲得較高的稀疏重構精度并且對噪聲的抑制效果較好。在本文中，依據Gram矩陣的理想單位陣形式及投影矩陣的結構，確定復高斯分布先驗信息，并且在多普勒維度的模型中，測量矩陣的維度較小，在處理過程中對應的乘法或取逆運算的運算量較小。因此，與文獻[9-11]中在 Welch界準則下直接求解投影矩陣的方法相比，應用Bayes方法可以實現精度和效率二者之間的均衡，依靠復高斯分布先驗信息可以獲得更高的精度，同時本文對測量矩陣的降維處理也降低了運算的復雜度，所以本文從構造Bayes模型的角度考慮優化測量矩陣。在第s次迭代中，有

在上述目標函數下，需要做的就是在每次迭代中 Gram矩陣已確定的情況下獲取對投影矩陣的最優解的精確求取。E中元素服從方差為γ的零均值復高斯分布，即Ea,b～CN(0,γI)，對應的概率密度函數為

投影矩陣Φ中元素也服從方差為β的零均值復高斯分布，即Φμ,κ～ CN(0,βI)，對應概率密度函數為

由此可得各參數的聯合概率密度函數為

其中，復高斯分布的方差在[0,∞)范圍內，即P(γ)∝1且P(β)∝1，因此式(26)可繼續化簡為

對式(28)取負對數處理可以得到各參數的目標函數，如式(29)所示。

經過推導，式(23)所示的原始問題轉換為目標函數L(Φ,γ,β)的極小值求解問題，因此需將L(Φ,γ,β)對各參數求偏導，且其聯合估計是以交替更新的方式實現的，以s表示第s次迭代，可得

3.3 運算復雜度分析

針對式(12)中的模型，若直接采用OMP算法，復雜度為O(I× (MNRNCQ) ×(UXUYUZUF))；若采用貝葉斯壓縮感知（BCS, Bayesian compressive sensing）方法，復雜度為O(I2 ×(UXUYUZ UF))，其中I是算法運行的迭代次數。采用文獻[16]的 DR-KCS方法時，其運算復雜度的高低與匹配濾波初始化得到的目標范圍有關，DR-KCS應用于OMP算法時的復雜度為O(I× (MNRNC Q) ×P0)，應用于貝葉斯壓縮感知方法時復雜度為O(I2×P0) ，其中P0與P的物理意義一樣，都是經優化縮減后的目標范圍，在高/低信噪比下P0的大小會有變化，進而會對運算復雜度產生影響。采用降維 Tucker分解[17]時，其復雜度近似為O(I× (MNRNCQ) ×(UXUYUZUF))。

本文方法的運算量分為兩部分，一部分是在各維的搜索過程中，以X維為例，若采用OMP算法，復雜度為O(I×M×UX)，若采用貝葉斯壓縮感知方法，其復雜度為O(I2×UX)，由此可見這個過程中的運算量極小，遠低于下一過程，可忽略不計；另一部分是在式(19)的稀疏重構過程中，若采用OMP算法，復雜度為O(I× (MNRNCQ) ×P)，若采用貝葉斯壓縮感知方法，復雜度為O(I2×P)。

本文方法與DR-KCS方法的復雜度相近，而與其他算法相比，由于P?U，本文方法的運算量比較小，對運算效率的提升比較明顯。

4 仿真實驗與分析

本節的仿真旨在證明本文方法可以有效地降低高運算量對成像性能的影響，實現高效的成像。實驗參數設置如下：發射陣元個數M=15，接收陣元行數與列數分別為NR=15，NC=15，發射陣元間距dt=2 m，接收陣元間距dr=1 m，信號載波頻率fc=10 GHz，信號帶寬B=50 MHz，發射脈沖寬度Tr=2 μs，采樣次數Q=10，成像區域參考中心到收發陣列的距離為 10 km，并且依據文獻[21]中對成像區域范圍的限定，取X、Y、Z維的最大成像范圍分別為-50～50 m，-50～50 m和 100～200 m，在該區域內設定由5個散射強度為1.0的強散射點構成，對應點的X、Y、Z三維坐標分別為(-20,10,150)、(0,-20,140)、(-30,0,130)、(20,-10,120)和(10,40,180)，5個目標點對應的多普勒頻率分別為 0、666.7 Hz、1 166.7 Hz、1 500 Hz和 166.7 Hz。

圖3描述的是在對測量矩陣優化前后各列之間的互相干性，其中，橫坐標是各列互相干性的絕對值，互相干性絕對值為0表示兩列完全不相干，其中互相干性絕對值最大值為 1，縱坐標表示與某一互相干性的值對應的列數。從圖 3(a)中可以看出，在測量矩陣優化前，互相干性的絕對值為0～1的分布比較雜亂，并且相對來說，絕對值比較大，說明此時測量矩陣的互相干特性比較差，這是由于未對多普勒維度進行針對性優化導致的。而在圖3(b)中，大多數列之間的相干性的絕對值為0，即大多數列之間互不相干，說明經過所提貝葉斯方法的優化處理后，測量矩陣的性能有了較大的提升，更有利于保證良好的稀疏恢復。

圖4和圖5是當信噪比SNR=30 dB時，本文所提方法的成像仿真結果，四維信息難以從圖像中直接呈現，因此直接選取如圖 4所示的XYFd、XZFd三維圖像及如圖5所示的XY、XZ、YZ這3個平面的投影來觀察成像效果。本文所提方法可以很精確地估計散射點的位置和幅值，有很高的恢復概率。

圖3 優化前后測量矩陣互相干性對比

圖4 三維成像結果

圖5 二維投影成像結果

圖6展示了SNR=30 dB且目標個數為5時本文方法的估計結果與迭代次數的關系，其中實驗參數與上述幾次實驗相同，針對性地從成像RMSE（root mean squared error）隨迭代次數的變化關系來分析本文方法的收斂性，并且從圖6中可以看出，本文方法在迭代5次左右基本完成收斂且誤差逼近于0，說明有較好的收斂特性。

圖6 成像誤差隨迭代次數變化關系

圖 7(a)是描述成像誤差與 SNR之間關系的曲線，其中本文方法應用BCS方法進行稀疏恢復。在SNR從-20～20 dB的范圍內，對本文方法與OMP、MAP這 2種傳統成像算法的稀疏成像性能做了比較，可以看出本文方法的成像性能稍有下降，這是因為重構測量矩陣時舍棄了大量的數據量，成像精度勢必會受到影響。圖7(b)是在SNR=30 dB下比較本文方法與MAP方法運行效率與稀疏度變化的關系，明顯可看出本文方法效率極高，會節省很多的算法運行時間。

圖7(c)是SNR在-20～20 dB的范圍內對本文方法與 Matrix Completion (MC)[15]及降維壓縮感知（DR-KCS, dimension-reduction kronecker compressive sensing）[16]這2種運算效率提升算法的成像性能對比，其中，本文方法與 DR-KCS方法均應用BCS方法進行稀疏恢復。從圖7(c)中可看出，本文方法的精度更高，這是由于在各維的搜索本身就是稀疏恢復的過程，對噪聲有一定的抑制作用，最后再經對式(19)的稀疏重構，因此會有更準確的稀疏成像效果。

圖7(d)和圖7(e)分別在不同的稀疏度和信噪比下比較了本文方法與DR-KCS方法的運行效率，圖7(d)中SNR=30 dB，圖7 (e)中目標數為5個。從圖7(d)和圖7(e)中可以看出，本文方法在目標個數較少的應用背景下有良好的處理效率，并且基于其各維搜索時的噪聲抑制過程，使其在低信噪比下的性能明顯優于DR-KCS方法。但目標個數多時算法效率會有所下降，這是由于新構造的測量矩陣的維度與目標個數成倍數關系，在多目標情況下測量矩陣列數會增大，因此使效率下降。

綜合來看，由圖7(a)和圖7(b)對比可知，本文方法對一般的成像方法而言，運算效率提升的效果極為明顯，雖精度稍有下降但可以滿足實際的應用。在圖7 (c)的對比中，本文方法與各效率提升算法相比又有更高的成像精度。在圖7(d)和圖7(e)的對比中，能夠展現本文方法的運算高效性，尤其在目標個數較少及低信噪比下可以較好地實現對算法運算效率提升的目的。因此從整體來看，本文方法可以達到本文的設計初衷，即在盡量保證成像效果的條件下盡可能地提升算法的運算效率。

5 結束語

圖7 本文方法與各類算法的性能對比結果

本文主要的貢獻在于基于 MIMO雷達提出了一種針對遠場運動目標的高效三維成像方法。從對測量矩陣處理的角度出發，依據回波模型依次在各維完成對目標的搜索，以初步縮小目標區域范圍，實現對算法運算效率的提升。同時分析多普勒采樣對測量矩陣相干性的影響，構造合適的投影矩陣。實現對整體測量矩陣互相干特性的優化。仿真結果表明，該方法可以在保證成像效果的條件下提升算法的運算效率，具有精確穩定的成像性能。