雷達(dá)成像近似二維模型及其超分辨算法

圖1　二維雷達(dá)目標(biāo)幾何圖

　　由于觀測(cè)角Δθ很小，取近似sin(nδθ)≈nδθ和cos(nδθ)≈1，則式(2)可近似寫(xiě)成:

　(3)

式中
　　式(3)指數(shù)項(xiàng)中的第三項(xiàng)是時(shí)頻耦合項(xiàng)，它是線頻調(diào)信號(hào)(其模糊函數(shù)為斜橢圓)所特有的，如果采用窄脈沖發(fā)射，則該項(xiàng)不存在.將該項(xiàng)忽略，則式(3)成為常用的回波二維正弦信號(hào)模型.
　　實(shí)際上，式(3)的第三項(xiàng)系“距離移動(dòng)”項(xiàng)，它與散射點(diǎn)的橫坐標(biāo)xk成正比，目標(biāo)區(qū)域大時(shí)必須考慮，而且這還遠(yuǎn)遠(yuǎn)不夠，散射點(diǎn)的多普勒移動(dòng)也必須考慮.為此，令sin(nδθ)≈nδθ和cos(nδθ)≈1-(nδθ)2/2，則式(2)較精確的近似式可寫(xiě)成：

　(4)

式(4)與式(3)相比較，指數(shù)中增加了兩項(xiàng)，其中前一項(xiàng)是“多普勒移動(dòng)”項(xiàng)，縱坐標(biāo)yk越大，影響也越大，這可以補(bǔ)充式(3)之不足；而后項(xiàng)是時(shí)頻耦合的多普勒移動(dòng)項(xiàng)，由于Mγ/Fs<

　(5)

　　需要指出，每個(gè)散射點(diǎn)的參數(shù)之間存在下述關(guān)系：ωk/μk=2γ/Fsfcδθ2和k/vk=fcFs/γδθ.由于雷達(dá)參數(shù)(fc,γ,Fs)和運(yùn)動(dòng)參數(shù)(δθ)均已知，所以待估計(jì)的五個(gè)參數(shù)中只有三個(gè)是獨(dú)立的.本文假設(shè)五個(gè)參數(shù)是獨(dú)立的，而在成像計(jì)算中已考慮參數(shù)之間的關(guān)系.
　　設(shè){ξk}Kk=1≡{αk,ωk,k,μk,vk}Kk=1，現(xiàn)在我們要從y(m,n)中估計(jì)參量{ξk}Kk=1.

三、二維推廣的RELAX算法
　　對(duì)于(5)式所示的信號(hào)模型，令：

Y=［y(m,n)］M×N

則　(6)
式中

　　設(shè)ξk估計(jì)值為，則ξk的估計(jì)問(wèn)題可通過(guò)優(yōu)化下述代價(jià)函數(shù)解決：

　(7)

式中‖.‖F(xiàn)表示矩陣的Frobenius范數(shù)，⊙表示矩陣的Hadamard積.
　　上式中C1的最優(yōu)化是一個(gè)多維空間的尋優(yōu)問(wèn)題，十分復(fù)雜.本文將RELAX［3］算法推廣以求解.為此，首先做以下準(zhǔn)備工作，令：

　(8)

　　即假定{i}i=1,2,…,K,i≠k已經(jīng)求出，則式(7)C1的極小化等效于下式的極小化：

C2(ξk)=‖Yk-αk(aM(ωk)bTN(k)Pk)⊙Dk(vk)‖2F　(9)

令：　　Zk=YkP-1k⊙Dk(-vk)　(10)
　　由于Pk為酉矩陣，矩陣Dk的每個(gè)元素的模|Dk(m,n)|=1,顯然矩陣Yk與Zk的F范數(shù)相同，故C2的極小化等效于下式的極小化：

C3=‖Zk-αkaM(ωk)bTN(k)‖2F　(11)

　　對(duì)上式關(guān)于αk求極小值就獲得αk的估計(jì)值k：

k=aHM(ωk)Zkb*N(k)/(MN)　(12)

　　從式(12)可以看出：是Zk歸一化的二維離散傅里葉變換在{ωk,k}處的值，所以只要得到估計(jì)值{k,k,k,k}，即可通過(guò)2D-FFT獲得k.
　　將估計(jì)值k代入式(11)后，估計(jì)值{k,k,k,k}可由下式尋優(yōu)得到：

　(13)

　　由上式可見(jiàn)，對(duì)于固定的{μk,vk}取值,估計(jì)值{k,k}為歸一化的周期圖|aHM(ωk)Zkb*N(k)|2/(MN)主峰處的二維頻率值.這樣，式(13)的優(yōu)化問(wèn)題歸結(jié)為：在(μk,vk)平面上可能的取值范圍內(nèi)尋找一點(diǎn){k,k}，在該點(diǎn)處周期圖|aHM(ωk)Zkb*N(k)|2/(MN)的主峰值比其余各點(diǎn)處的主峰值都大.所以，我們通過(guò)上述二維尋優(yōu)獲得{μk,vk}的估計(jì)值{k,k}，再由式(13)得到{ωk,k}的估計(jì)值{k,k}.
　　實(shí)際中，為了加快運(yùn)算速度，二維(μk,vk)平面的尋優(yōu)可以用Matlab中的函數(shù)Fmin()實(shí)現(xiàn).
　　在做了以上的準(zhǔn)備工作以后，基于推廣的RELAX算法的參量估計(jì)步驟如下：
　　第一步：假設(shè)信號(hào)數(shù)K=1,分別利用式(13)和式(12)計(jì)算1.
　　第二步(2)：假設(shè)信號(hào)數(shù)K=2，首先將第一步計(jì)算所得到的1代入式(8)求出Y2，再利用式(13)和式(12)計(jì)算2；將計(jì)算的2代入式(8)求出Y1，然后利用式(13)和式(12)重新計(jì)算1，這個(gè)過(guò)程反復(fù)疊代，直至收斂.
　　第三步:假設(shè)信號(hào)數(shù)K=3,首先將第二步計(jì)算所得到的1和2代入式(8)求出Y3，再利用式(13)和式(12)計(jì)算3；將計(jì)算的3和2代入式(8)求出Y1，然后利用式(13)和式(12)重新計(jì)算1；將計(jì)算的1和3代入式(8)求出Y2，然后利用式(13)和式(12)重新計(jì)算2，這個(gè)過(guò)程反復(fù)疊代，直至收斂.
　　剩余步驟：令K=K+1,上述步驟持續(xù)進(jìn)行，直到K等于待估計(jì)信號(hào)數(shù).
　　上述過(guò)程中的收斂判據(jù)與RELAX算法的收斂判據(jù)相同，即比較代價(jià)函數(shù)C1在兩次疊代過(guò)程中的變化值，如果這個(gè)變換值小于某個(gè)值，如ε=10-3，則認(rèn)為過(guò)程收斂.

四、數(shù)值模擬
　　1.算法參數(shù)估計(jì)性能模擬
　　模擬數(shù)據(jù)由式（5）產(chǎn)生，M=10,N=10，信號(hào)數(shù)K=2.信號(hào)參數(shù)和實(shí)驗(yàn)條件如表1所示,為復(fù)高斯白噪聲.注意兩信號(hào)的頻率差小于FFT的分辨率Δf=Δω/(2π)=0.1.表1給出了信號(hào)參數(shù)估計(jì)均方根誤差的統(tǒng)計(jì)結(jié)果及相應(yīng)情形時(shí)的C-R界，可見(jiàn)，估計(jì)均方根誤差與CR界十分接近.另外表中還給出了估計(jì)均值，與真實(shí)值也非常接近.

表1　二維信號(hào)的參數(shù)估計(jì)、CRB及與均方根差的比較

圖2　距離走動(dòng)誤差下的RELAX成像結(jié)果

圖3　距離走動(dòng)誤差下的

圖4　RELAX方法估計(jì)的信號(hào)強(qiáng)度推廣RELAX成像結(jié)果

圖5　推廣RELAX方法估計(jì)的信號(hào)強(qiáng)度

附 錄：參數(shù)估計(jì)的C-R界
　　下面我們給出式(5)所示的二維信號(hào)參量估計(jì)的C-R界表達(dá)式.同時(shí)假設(shè)式(5)中加性噪聲為零均值高斯色噪聲，其協(xié)方差矩陣未知.令：

y=vec(Y)　(A.1)
e=vec(E)　(A.2)
dk=vec(Dk)　(A.3)

式中vec(X)=(xT1,xT2，…，xTN)T，向量xn(n=1,2,…,N)為矩陣X的列向量.我們將式(5)改寫(xiě)為如下向量形式：

　(A.4)

式中表示Kronecker積,Ω=［{［P1bN(1)］aM(ω1)}⊙d1…{［PkbN(K)］aM(ωK)}⊙dK］,α=(α1，α2，…，αK)T.
　　令Q=E(eeH)為e的協(xié)方差矩陣，則對(duì)于由式(A.4)所示的二維信號(hào)模型，其Fisher信息陣(FIM)的第ij個(gè)元素推廣的Slepian-Bangs公式為［5,6］:
(FIM)ij=tr(Q-1Q′iQ-1Q′j)+2Re［(αHΩH)′iQ-1(Ωα)′j］　(A.5)
式中X′i表示矩陣X對(duì)第i個(gè)參數(shù)求導(dǎo)，tr(X)為矩陣的跡，Re(X)為矩陣的實(shí)部.由于Q與Ωα中的參量無(wú)關(guān)，而Ωα亦與Q的元素?zé)o關(guān)，顯然FIM為一塊對(duì)角陣.所以待估計(jì)參量的C-R界矩陣由(A.5)式的第二項(xiàng)得到.

令：η=(［Re(α)］T［Im(α)］TωTTμTvT)T　(A.6)

式中ω=(ω1，ω2，…，ωK)T,μ=(μ1,μ2,…,μK)T,=(1，2，…，K)T,v=(v1,v2,…,vK)T.
令：F=［Ω　jΩ　DωΘ　DΘ　DμΘ　DvΘ］　(A.7)
式中矩陣Dω、D、Dμ、Dv的第k列分別為：［{［PkbN(k)］aM(ωk)}⊙dk］/ωk、［{［PkbN(k)］aM(ωk)}⊙dk］/k、［{［PkbN(k)］aM(ωk)}⊙dk］/μk、［{［PkbN(k)］aM(ωk)}⊙dk］/vk，Θ=diag{α1　α2　…　αK}.則關(guān)于參量向量η的CRB矩陣為

CRB(η)=［2Re(FHQ-1F)］-1　(A.8)