本發(fā)明涉及的是探地雷達(dá)信號處理技術(shù)、合成孔徑成像技術(shù)、災(zāi)害地質(zhì)調(diào)查、考古調(diào)查、公路工程質(zhì)量檢測以及地質(zhì)勘探等領(lǐng)域，特別是涉及一種基于主成分分析的加權(quán)后向投影算法的成像方法。

背景技術(shù)：

探地雷達(dá)(Ground Penetrating Radar,GPR)是用來對地下目標(biāo)或場景進(jìn)行探測的雷達(dá)系統(tǒng)。利用電磁波在地下介質(zhì)電磁特性不連續(xù)之處產(chǎn)生的反射和散射現(xiàn)象，GPR可以獲取地下目標(biāo)或場景的幾何和物理信息。與傳統(tǒng)的對地探測方式相比，GPR具有不破壞探測場景、可重復(fù)性強(qiáng)、采樣密集、安全可靠、高分辨率、操作靈活、經(jīng)濟(jì)性強(qiáng)等優(yōu)勢。直接使用一維或二維回波來完成檢測和識別工作不方便直觀，因此GPR成像技術(shù)應(yīng)運而生，但是受到地下干擾嚴(yán)重，成像效果不佳。

在探地雷達(dá)成像方面，文獻(xiàn)《基于Stolt偏移的探地雷達(dá)合成孔徑成像研究》主要是通過一種新的stolt偏移插值實現(xiàn)方法，用于探地雷達(dá)合成孔徑成像，克服傳統(tǒng)方法中偏移能量不集中的缺點并保持處理速度快的優(yōu)點，文獻(xiàn)《實現(xiàn)淺地層探地雷達(dá)快速合成孔徑成像的一種有效方法》主要是通過基于A-scan能量的特征來確定目標(biāo)回波數(shù)據(jù)在整個回波數(shù)據(jù)中的分布區(qū)域，并只用該區(qū)域數(shù)據(jù)參與探地雷達(dá)合成孔徑成像運算，從而排除非目標(biāo)回波數(shù)據(jù)對合成孔徑成像的影響，減少合成孔徑運算量，實現(xiàn)快速合成孔徑成像，文獻(xiàn)《基于因式分解BP的車載前視地表穿透SAR快速成像算法》主要是針對雙站SAR模型分析FFBP的距離誤差，給出了極角采樣率限制條件，提出了雙站FFBP算法，通過合理設(shè)置子孔徑中心提高聚焦效果，并結(jié)合距離相位補(bǔ)償和線性差值技術(shù)，進(jìn)一步提高聚焦精度。文獻(xiàn)《基于幅度加權(quán)的探地雷達(dá)后向投影算法》主要是利用邊緣提取獲得雙曲線回波上的點對應(yīng)的響應(yīng)幅值構(gòu)造窗函數(shù)，然后對各通道的散射響應(yīng)幅值進(jìn)行加權(quán)成像，文獻(xiàn)《一種基于Hilbert變換的探地雷達(dá)成像方法》主要是通過EEMD方法分解出本征模量，然后進(jìn)行Hilbert-Huang變換，求出順時振幅進(jìn)行疊加成像，文獻(xiàn)《一種用于探地雷達(dá)成像的散射強(qiáng)度加權(quán)處理方法》利用成像區(qū)域中對應(yīng)的時延曲線處的散射數(shù)據(jù)求加權(quán)系數(shù)，而本發(fā)明利用主成分分析法求加權(quán)系數(shù)去提高成像效果，起到較好的旁瓣和雜波抑制效果，便于對地下目標(biāo)進(jìn)行檢測和識別。

技術(shù)實現(xiàn)要素：

本發(fā)明的目的在于克服現(xiàn)有技術(shù)的缺點與不足，提供了一種基于主成分分析的加權(quán)后向投影算法的成像方法，起到較好的旁瓣和雜波抑制效果。

本發(fā)明的目的通過以下技術(shù)方案實現(xiàn)：一種基于主成分分析的加權(quán)后向投影算法的成像方法，包括如下步驟：

步驟1、探地雷達(dá)采集數(shù)據(jù)，對采集到的數(shù)據(jù)信號進(jìn)行預(yù)處理；

步驟2、將預(yù)處理后的結(jié)果通過時延估計得到時延坐標(biāo)；將時延坐標(biāo)位置取窗得到時延位置的幅值區(qū)域；

步驟3、利用主成分分析法求各道時延位置幅值區(qū)域的能量和作為加權(quán)系數(shù)；最后利用加權(quán)系數(shù)和后向投影算法疊加得到成像結(jié)果。

進(jìn)一步的，所述預(yù)處理包括噪聲抑制、天線耦合波、直達(dá)波抑制和射頻干擾抑制。

進(jìn)一步的，所述步驟2具體為：

1)直線x＝0將場景分為兩部分，上半部分為空氣層，波速為光速c，下半部分為各向同性的均勻介質(zhì)，波速其中ε_r代表媒質(zhì)的相對介電常數(shù)；

2)收發(fā)天線的位置(x_k,-h)，A(x_A,z_A)為目標(biāo)所在位置，(x_r,0)表示折射點位置，根據(jù)Snell折射定律，有如下關(guān)系式成立：

其中θ_i和θ_r分別表示入射角和折射角，利用幾何關(guān)系，可知：

故可得到下面一元四次方程：

通過二分法求解此一元四次方程組，求得折射點的x_r值，則各道回波上對應(yīng)于A點的時延分別為：

求解一組時間單元方向上的時延坐標(biāo)，利用平移不變性原理便可完成整個場景中的時延計算；

3)利用發(fā)射信號的中心頻率f₀和等效采樣頻率Fs，求得散射響應(yīng)范圍利用此范圍結(jié)合時延位置，得到時延點對應(yīng)的散射響應(yīng)矩陣P：

K表示合成孔徑的長度,其中表示第k個等效陣元位置到A點對應(yīng)在時延位置的幅值，表示在第k道數(shù)據(jù)中距離第i個位置的幅值，若下坐標(biāo)超過原數(shù)據(jù)的范圍，則將其設(shè)為0。

進(jìn)一步的，所述步驟3具體為：利用主成分分析法可得P矩陣中相關(guān)性強(qiáng)的能量為其最近秩1矩陣，計算最近秩1矩陣，并將該矩陣的能量和作為加權(quán)系數(shù)，具體為：

對P矩陣進(jìn)行奇異值分解：

其中U＝{U_m,i}＝{u₁,u₂,…,u_S+1}∈R^(S+1)×(S+1),V＝{V_n,i}＝{v₁,v₂,…,v_K}∈R^K×K，U、V均為正交矩陣，U矩陣是屬于大小為(S+1)*(S+1)的實數(shù)矩陣，V矩陣同理，T為轉(zhuǎn)置符號，R為實數(shù)，D∈R^S+1×K為對角陣，對角元素即為奇異值，按由大到小排列，可以表示為D_1,1≥D_2,2≥…≥D_L,L；

令σ_i＝D_i,i，由奇異值分解的性質(zhì)可知為矩陣P的“最近”秩l矩陣；故其最近秩1矩陣可表示為：

P⁽¹⁾＝σ₁W₁

利用最近秩1矩陣內(nèi)所有值的能量和作為加權(quán)系數(shù)，即

利用加權(quán)系數(shù)和后向投影算法將點A在各道上的散射響應(yīng)幅值進(jìn)行疊加，完成對A點的成像E_A：

附圖說明

圖1本發(fā)明的方法流程圖；

圖2探地雷達(dá)原數(shù)據(jù)B-scan圖；

圖3探地雷達(dá)信號預(yù)處理數(shù)據(jù)B-scan圖；

圖4 GPR成像BP算法場景；

圖5成像點位置關(guān)系的平移不變性；

圖6時延曲線的平移不變性

圖7探地雷達(dá)后向投影算法成像圖；

圖8主成分分析成像圖；

圖9基于主成分分析的加權(quán)后向投影算法成像圖；

具體實施方式

下面將結(jié)合本發(fā)明實施例中的附圖對本發(fā)明實施例中的技術(shù)方案進(jìn)行清楚、完整地描述，顯然，所描述的實施例僅僅是本發(fā)明一部分實施例，而不是全部的實施例?；诒景l(fā)明中的實施例，本領(lǐng)域普通技術(shù)人員在沒有做出創(chuàng)造性勞動前提下所獲得的所有其他實施例，都屬于本發(fā)明保護(hù)的范圍。

結(jié)合圖1，一種基于主成分分析的加權(quán)后向投影算法的成像方法，包括如下步驟：

步驟1、探地雷達(dá)采集數(shù)據(jù)，對采集到的數(shù)據(jù)信號進(jìn)行預(yù)處理；圖2給出了預(yù)處理前的探地雷達(dá)原數(shù)據(jù)，圖3給出了預(yù)處理后的探地雷達(dá)數(shù)據(jù)；

步驟2、將預(yù)處理后的結(jié)果通過時延估計得到時延坐標(biāo)；將時延坐標(biāo)位置取窗得到時延位置的幅值區(qū)域；

步驟3、利用主成分分析法求各道時延位置幅值區(qū)域的能量和作為加權(quán)系數(shù)；最后利用加權(quán)系數(shù)和后向投影算法疊加得到成像結(jié)果。

其中預(yù)處理后通過計算各道數(shù)據(jù)的時延和計算中心頻率對應(yīng)時間單元數(shù)目，保存時延位置對應(yīng)的幅值和幅值區(qū)域值，利用主成分分析法計算互相關(guān)系數(shù)作為后向投影算法加權(quán)系數(shù)，從而得到成像結(jié)果。

所述預(yù)處理包括噪聲抑制、天線耦合波、直達(dá)波抑制和射頻干擾抑制。

結(jié)合圖4，所述步驟2具體為：

1)直線x＝0將場景分為兩部分，上半部分為空氣層，波速為光速c，下半部分為各向同性的均勻介質(zhì)，波速其中ε_r代表媒質(zhì)的相對介電常數(shù)；

2)收發(fā)天線的位置(x_k,-h)，A(x_A,z_A)為目標(biāo)所在位置，(x_r,0)表示折射點位置，根據(jù)Snell折射定律，有如下關(guān)系式成立：

其中θ_i和θ_r分別表示入射角和折射角，利用幾何關(guān)系，可知：

故可得到下面一元四次方程：

通過二分法求解此一元四次方程組，求得折射點的x_r值，則各道回波上對應(yīng)于A點的時延分別為：

結(jié)合圖5、圖6，求解一組時間單元方向上的時延坐標(biāo)，利用平移不變性原理便可完成整個場景中的時延計算；

3)利用發(fā)射信號的中心頻率f₀和等效采樣頻率Fs，求得散射響應(yīng)范圍利用此范圍結(jié)合時延位置，得到時延點對應(yīng)的散射響應(yīng)矩陣P：

K表示合成孔徑的長度,其中表示第k個等效陣元位置到A點對應(yīng)在時延位置的幅值，表示在第k道數(shù)據(jù)中距離第i個位置的幅值，若下坐標(biāo)超過原數(shù)據(jù)的范圍，則將其設(shè)為0。

結(jié)合圖7、圖8、圖9，所述步驟3具體為：利用主成分分析法可得P矩陣中相關(guān)性強(qiáng)的能量為其最近秩1矩陣，計算最近秩1矩陣，并將該矩陣的能量和作為加權(quán)系數(shù)，具體為：

對P矩陣進(jìn)行奇異值分解：

其中U＝{U_m,i}＝{u₁,u₂,…,u_S+1}∈R^(S+1)×(S+1),V＝{V_n,i}＝{v₁,v₂,…,v_K}∈R^K×K，U、V均為正交矩陣，U矩陣是屬于大小為(S+1)*(S+1)的實數(shù)矩陣，V矩陣同理，^T為轉(zhuǎn)置符號，R為實數(shù)，D∈R^S+1×K為對角陣，對角元素即為奇異值，按由大到小排列，可以表示為D_1,1≥D_2,2≥…≥D_L,L；

令σ_i＝D_i,i，由奇異值分解的性質(zhì)可知為矩陣P的“最近”秩l矩陣；故其最近秩1矩陣可表示為：

P⁽¹⁾＝σ₁W₁

利用最近秩1矩陣內(nèi)所有值的能量和作為加權(quán)系數(shù)，即

利用加權(quán)系數(shù)和后向投影算法將點A在各道上的散射響應(yīng)幅值進(jìn)行疊加，完成對A點的成像E_A：

利用主成分分析方法加權(quán)后向投影算法的成像結(jié)果如圖9所示，圖7為探地雷達(dá)后向投影算法成像圖，圖8為主成分分析方法成像圖，比較圖7和圖9可以發(fā)現(xiàn)，通過加權(quán)的后向投影算法可以有效抑制干擾和噪聲，提高信噪比和目標(biāo)檢測概率。

以上對本發(fā)明所提供的一種一種基于主成分分析的加權(quán)后向投影算法的成像方法，進(jìn)行了詳細(xì)介紹，本文中應(yīng)用了具體個例對本發(fā)明的原理及實施方式進(jìn)行了闡述，以上實施例的說明只是用于幫助理解本發(fā)明的方法及其核心思想；同時，對于本領(lǐng)域的一般技術(shù)人員，依據(jù)本發(fā)明的思想，在具體實施方式及應(yīng)用范圍上均會有改變之處，綜上所述，本說明書內(nèi)容不應(yīng)理解為對本發(fā)明的限制。