趙子恒 林 赟 王彥平 李 洋 申文杰

（北方工業(yè)大學(xué)信息學(xué)院雷達(dá)監(jiān)測(cè)技術(shù)實(shí)驗(yàn)室，北京 100144）

1 引言

地基合成孔徑雷達(dá)（Ground-based Synthetic Aperture Radar，GBSAR）作為一種新型高精度測(cè)量技術(shù)，具有非接觸、全天候、設(shè)站靈活等優(yōu)點(diǎn)，能夠?qū)崟r(shí)監(jiān)測(cè)場(chǎng)景整體形變情況，已成為國(guó)內(nèi)外各監(jiān)測(cè)領(lǐng)域研究的熱點(diǎn)［1-2］。傳統(tǒng)GBSAR 系統(tǒng)通常采用直線軌道，該系統(tǒng)在滑坡監(jiān)測(cè)預(yù)警和礦山橋梁形變監(jiān)測(cè)等領(lǐng)域的技術(shù)適應(yīng)性已經(jīng)得到充分的應(yīng)用檢驗(yàn)。然而該系統(tǒng)只能獲取二維圖像和一維形變，且在地形起伏區(qū)域存在疊掩現(xiàn)象。為了解決上述問(wèn)題，學(xué)者們提出了多基線直線軌道GBSAR 和圓周掃描GBSAR（Ground-based Circular Synthetic Aperture Radar，GBCSAR）兩種具備三維分辨能力的GBSAR模式［3-5］，國(guó)內(nèi)外研究機(jī)構(gòu)經(jīng)過(guò)大量實(shí)驗(yàn)工作后證明了GBCSAR 具備三維分辨能力，且可以有效解決二維成像中的疊掩現(xiàn)象。GBCSAR的合成孔徑通過(guò)附在旋臂末端上的相位中心旋轉(zhuǎn)生成，在保證系統(tǒng)分辨率的前提下，能夠在一次旋轉(zhuǎn)觀測(cè)中覆蓋監(jiān)測(cè)區(qū)域的全部場(chǎng)景，實(shí)現(xiàn)三維成像。然而GBCSAR 獨(dú)特的運(yùn)動(dòng)軌跡給成像算法的實(shí)現(xiàn)增加了難度，目前對(duì)于GBCSAR，學(xué)者們通常采用后向投影（Back Projection，BP）成像算法。該算法成像精度高，且可以應(yīng)用于多種架構(gòu)的雷達(dá)而不受天線運(yùn)動(dòng)軌跡所限制。文獻(xiàn)［5］中提出了GBCSAR 模式下的三維層析BP算法，該算法通過(guò)疊加距離向上不同距離處的二維圖像實(shí)現(xiàn)三維成像［6］。但是該成像算法未曾應(yīng)用于大場(chǎng)景成像，且面對(duì)多目標(biāo)點(diǎn)的場(chǎng)景時(shí)，層析BP算法數(shù)據(jù)量巨大，存在冗余現(xiàn)象，成像效率較低。在連續(xù)監(jiān)測(cè)時(shí)數(shù)據(jù)量更大，無(wú)法滿足實(shí)時(shí)形變監(jiān)測(cè)的需求［7-8］。

因此本文提出了一種地形表面BP 成像算法，該算法首先通過(guò)層析BP 算法獲取場(chǎng)景目標(biāo)點(diǎn)的幅值三維矩陣，利用幅值矩陣提取出地形有效點(diǎn)的DEM 數(shù)據(jù)，而后在有效點(diǎn)三維曲面上實(shí)現(xiàn)BP 成像，不再借助參考平面投影成像［9］。在連續(xù)地形形變監(jiān)測(cè)時(shí)，同一位置不同高度存在很多目標(biāo)點(diǎn)，它們的形變趨勢(shì)基本一致，因此可以將這些同一位置不同高度的多個(gè)目標(biāo)點(diǎn)用一個(gè)有效點(diǎn)代替，這個(gè)有效點(diǎn)的成像數(shù)據(jù)可以有效地反映該位置不同高度處多個(gè)目標(biāo)點(diǎn)的成像數(shù)據(jù)［10-12］。這樣將場(chǎng)景中所有目標(biāo)點(diǎn)的三維矩陣降維成一個(gè)由各個(gè)有效點(diǎn)接連的三維曲面，極大地減少了數(shù)據(jù)量，且三維曲面上的BP成像結(jié)果可以反映出監(jiān)測(cè)場(chǎng)景的形變信息。

本文其余部分安排如下：第2部分介紹GBCSAR的幾何模型和信號(hào)模型。第3部分對(duì)算法的原理進(jìn)行論述，并具體介紹地形表面有效點(diǎn)曲面的獲取方法以及如何在有效點(diǎn)曲面上實(shí)現(xiàn)BP 成像。第4 部分使用機(jī)載三維數(shù)據(jù)進(jìn)行仿真實(shí)驗(yàn)，對(duì)算法的可行性進(jìn)行驗(yàn)證，并與層析BP 算法進(jìn)行計(jì)算量比較。最后，結(jié)論總結(jié)。

2 GBCSAR的幾何模型和信號(hào)模型

GBCSAR 系統(tǒng)在旋臂的末端裝有收發(fā)天線，該系統(tǒng)通過(guò)旋臂的旋轉(zhuǎn)形成二維合成孔徑。旋臂在YOZ 平面沿逆時(shí)針旋轉(zhuǎn)，方位向和高度向都在YOZ平面內(nèi)，距離向與旋轉(zhuǎn)平面YOZ垂直。

GBCSAR 的幾何模型如圖1 所示。相位中心位于旋臂尾端a處，旋臂以坐標(biāo)系原點(diǎn)O為中心在YOZ 平面旋轉(zhuǎn)，監(jiān)測(cè)區(qū)域位于與YOZ 平面垂直方向。在圖1 中，a為相位中心位置，θ為相位中心的瞬時(shí)旋轉(zhuǎn)角，a坐標(biāo)隨著θ變化。因此a坐標(biāo)可以表示為(0，rcosθ，rsinθ)，其 中r為旋臂長(zhǎng)度。p1，p2，…pn分別代表不同目標(biāo)點(diǎn)，根據(jù)相位中心a坐標(biāo)(xa，ya，za)和目標(biāo)點(diǎn)pn坐標(biāo)(xn，yn，zn)可以得到目標(biāo)點(diǎn)與相位中心a之間的斜距Rn的表達(dá)式

為了方便計(jì)算，式（1）結(jié)合相位中心a坐標(biāo)可得到瞬時(shí)斜距Rn表達(dá)式

已知了GBCSAR 的幾何模型，接下來(lái)給出系統(tǒng)信號(hào)模型。這里以線性調(diào)頻信號(hào)（LFM）作為發(fā)射信號(hào)，通常距離向脈沖壓縮是在頻域執(zhí)行的，因此直接給出經(jīng)過(guò)頻域脈沖壓縮后回波的頻域表達(dá)式

其 中，σ是散射系數(shù)，Br是發(fā)射信號(hào)帶寬，f∈[-Br/2，Br/2]，fc是中心頻率，c是光速。在頻域脈沖壓縮后，需要對(duì)回波信號(hào)進(jìn)行逆傅里葉變換，由于之后的BP 算法是在時(shí)域處理信號(hào)，因此這里給出距離脈沖壓縮后信號(hào)的時(shí)域表達(dá)式

GBCSAR具備三維分辨能力，它的距離向分辨率為

其中Tr為雷達(dá)脈寬，又因?yàn)槊}寬Tr與帶寬Br互為倒數(shù)，因此（5）可表示為

由于GBCSAR 的方位向和高度向都在旋轉(zhuǎn)平面YOZ 中，因此在Y方向和Z方向上，GBCSAR 的分辨率相同。這里推導(dǎo)方位向分辨率δy，定義波數(shù)K用來(lái)表示頻率

定義比值αa來(lái)表示帶寬與頻率比，表達(dá)式如下

由此可以得到GBCSAR 在Y方向和Z方向上的分辨率

3 GBCSAR地形表面成像算法

對(duì)于處在同一(X，Y)坐標(biāo)上的多個(gè)目標(biāo)點(diǎn)，所在高度Z不同，導(dǎo)致各目標(biāo)點(diǎn)與雷達(dá)之間的斜距Rn和幅值大小f(X，Y，Z) 都不相同。其中幅值f(X，Y，Z)最大的目標(biāo)點(diǎn)Pn通常在對(duì)應(yīng)(X，Y)坐標(biāo)上多個(gè)目標(biāo)點(diǎn)中散射系數(shù)最大，這個(gè)點(diǎn)的形變信息可以基本反映出該(X，Y)坐標(biāo)上所有目標(biāo)點(diǎn)的形變趨勢(shì)。因此將同一(X，Y)坐標(biāo)上多個(gè)目標(biāo)點(diǎn)降維成一個(gè)目標(biāo)點(diǎn)，即該(X，Y)坐標(biāo)上的有效點(diǎn)，對(duì)這個(gè)有效點(diǎn)進(jìn)行相應(yīng)的成像就可以有效地實(shí)現(xiàn)這個(gè)位置的形變監(jiān)測(cè)。而將整個(gè)場(chǎng)景所有(X，Y)坐標(biāo)對(duì)應(yīng)的有效點(diǎn)提取出來(lái)后，在有效點(diǎn)曲面上進(jìn)行BP 成像，即可實(shí)現(xiàn)對(duì)整片場(chǎng)景的形變監(jiān)測(cè)。本文提出的方法兩個(gè)主要步驟分別是獲取地形有效點(diǎn)所構(gòu)成的三維曲面和在該曲面上實(shí)現(xiàn)BP成像。

本節(jié)介紹了如何通過(guò)目標(biāo)點(diǎn)三維坐標(biāo)獲取地形有效點(diǎn)DEM數(shù)據(jù)和有效點(diǎn)三維曲面，以及如何在有效點(diǎn)三維曲面上BP成像。

3.1 地形表面有效點(diǎn)曲面提取

在采用層析BP 算法獲取所有目標(biāo)點(diǎn)的幅值前需要對(duì)來(lái)自場(chǎng)景的SAR 回波數(shù)據(jù)進(jìn)行預(yù)處理，對(duì)接收到的實(shí)信號(hào)進(jìn)行正交解調(diào)得到算法需要的復(fù)信號(hào)。而后采用層析BP 算法獲取所有目標(biāo)點(diǎn)的幅值，得到目標(biāo)點(diǎn)幅值矩陣。在連續(xù)監(jiān)測(cè)時(shí)，只需在開(kāi)始時(shí)使用一次層析BP 算法提取有效點(diǎn)曲面，之后的連續(xù)監(jiān)測(cè)都在提取得到的有效點(diǎn)曲面上成像。

提取有效點(diǎn)時(shí)需要根據(jù)場(chǎng)景大小和目標(biāo)點(diǎn)的分布對(duì)幅值矩陣進(jìn)行插值或等間隔采樣，以確保目標(biāo)點(diǎn)密集分布時(shí)不會(huì)產(chǎn)生數(shù)據(jù)冗余，同時(shí)不會(huì)出現(xiàn)目標(biāo)點(diǎn)稀疏分布時(shí)提取到的有效點(diǎn)不連續(xù)，并且要確保有效點(diǎn)的高度差大于系統(tǒng)理論分辨率δz，否則會(huì)出現(xiàn)混疊現(xiàn)象。而后按照目標(biāo)點(diǎn)的(X，Y)坐標(biāo)遍歷查找對(duì)應(yīng)(X，Y)坐標(biāo)上的有效點(diǎn)，此時(shí)得到的有效點(diǎn)坐標(biāo)為(X，Y，Z(x，y))，即有效點(diǎn)的高度Z是根據(jù)對(duì)應(yīng)(X，Y)坐標(biāo)一一對(duì)應(yīng)的。當(dāng)所有(X，Y)遍歷結(jié)束后，整個(gè)場(chǎng)景的有效點(diǎn)會(huì)連接成一個(gè)三維曲面，這個(gè)曲面就是最終BP成像所在的投影曲面。

如圖3 所示，右側(cè)的三維曲面就是通過(guò)提取地形有效點(diǎn)得到的，需要注意的是這個(gè)曲面并不是場(chǎng)景的頂面，而是根據(jù)有效點(diǎn)的分布位置連接而成的。相比于左側(cè)的三維立方場(chǎng)景，三維曲面極大地減少了目標(biāo)點(diǎn)個(gè)數(shù)，減少了計(jì)算量。

3.2 有效點(diǎn)三維曲面BP成像

與傳統(tǒng)的二維BP算法和層析三維BP算法不同的是，傳統(tǒng)二維BP算法設(shè)定像素點(diǎn)都在同一高度平面上，通常選擇地平面。將整個(gè)場(chǎng)景中某一個(gè)高度的回波數(shù)據(jù)都投影在一個(gè)平面，因此傳統(tǒng)二維BP算法只能實(shí)現(xiàn)單一高度的成像，無(wú)法獲取整個(gè)場(chǎng)景信息。而層析三維BP 算法是將各個(gè)高度的回波數(shù)據(jù)分別投影在各個(gè)平面，需要在三維空間中每個(gè)高度上都設(shè)定等量的像素點(diǎn)，因此計(jì)算量龐大，無(wú)法滿足連續(xù)監(jiān)測(cè)的需求。層析BP算法示意圖如圖4所示。

本文提出的地形表面有效點(diǎn)成像則是依據(jù)有效點(diǎn)的高度在對(duì)應(yīng)高度設(shè)定像素點(diǎn)，沒(méi)有有效點(diǎn)的高度也就沒(méi)有像素點(diǎn)。最終像素點(diǎn)會(huì)連接成有起伏的三維曲面，并在該三維曲面上進(jìn)行投影成像。地形表面有效點(diǎn)成像示意圖如下圖5所示。

三維曲面上的BP 成像首先需要借助距離脈沖壓縮后的回波獲得目標(biāo)點(diǎn)幅值矩陣，借助幅值矩陣獲取有效點(diǎn)坐標(biāo)，然后根據(jù)有效點(diǎn)坐標(biāo)和DEM數(shù)據(jù)構(gòu)建成像網(wǎng)格。計(jì)算成像網(wǎng)格中每個(gè)像素點(diǎn)距離相位中心a的瞬時(shí)斜距，得到每個(gè)像素點(diǎn)的雙程延時(shí)以提取對(duì)應(yīng)像素點(diǎn)的幅值，最后疊加各個(gè)方位向在各像素點(diǎn)的幅值得到SAR 圖像。具體實(shí)現(xiàn)步驟如下：

步驟1：根據(jù)系統(tǒng)理論分辨率和提取到的三維曲面劃分成像網(wǎng)格；

步驟2：計(jì)算各個(gè)像素點(diǎn)距相位中心a當(dāng)前位置的雙程延時(shí)；

步驟3：通過(guò)插值從雷達(dá)回波中依據(jù)雙程延時(shí)提取當(dāng)前方位向上像素點(diǎn)的回波值；

步驟4：對(duì)得到的回波值進(jìn)行相位補(bǔ)償?shù)玫綄?duì)應(yīng)方位向上的子圖像；

步驟5：遍歷各個(gè)方位向，重復(fù)步驟3、4，疊加各個(gè)方位向的子圖像；

步驟6：得到最終的SAR圖像。

根據(jù)式（6）和（9）得到GBCSAR 的理論分辨率，X和Y方向設(shè)定的網(wǎng)格間隔應(yīng)略小于δx和δy，在提取有效點(diǎn)時(shí)已經(jīng)保證了有效點(diǎn)在Z方向的間隔大于理論分辨率δz。X和Y方向設(shè)定的成像網(wǎng)格個(gè)數(shù)根據(jù)場(chǎng)景大小設(shè)定，最終X和Y方向的成像網(wǎng)格范圍包含場(chǎng)景所有目標(biāo)點(diǎn)即可，成像網(wǎng)格在Z方向的分布和有效點(diǎn)在Z方向的分布一致。設(shè)定好的像素點(diǎn)坐標(biāo)可以表示為(xi，yj，z(xi，yj))，其中i，j表示XOY 平面中第i行第j列的像素點(diǎn)，像素點(diǎn)所在的高度z(xi，y j)是從有效點(diǎn)DEM 中遍歷出的，最終像素點(diǎn)都是分布在提取出的有效點(diǎn)附近。

根據(jù)相位中心a的坐標(biāo)(0，ya，za)和像素點(diǎn)的坐標(biāo)(xi，yj，z(xi，yj))計(jì)算出各個(gè)像素點(diǎn)到相位中心的距離

像素點(diǎn)到相位中心的雙程延時(shí)表示為：

由于回波數(shù)據(jù)是離散地存儲(chǔ)在二維矩陣中，回波值在矩陣中對(duì)應(yīng)的都是整數(shù)索引，而雙程延時(shí)在回波矩陣中對(duì)應(yīng)的回波值索引大多數(shù)不是整數(shù)，想要得到離散信號(hào)非整數(shù)點(diǎn)的值就需要對(duì)信號(hào)進(jìn)行插值。因此先對(duì)回波數(shù)據(jù)插值，再利用延時(shí)和回波的對(duì)應(yīng)關(guān)系從回波數(shù)據(jù)中取得對(duì)應(yīng)方位向的回波值，本文采用的是sinc插值，插值核為8，插值公式為：

插值后得到各像素點(diǎn)的回波值需要進(jìn)行相位補(bǔ)償，以實(shí)現(xiàn)各方位向的聚焦。將回波值與對(duì)應(yīng)像素點(diǎn)的相位補(bǔ)償因子相乘即可實(shí)現(xiàn)補(bǔ)償：

其中R(i，j)為當(dāng)前像素點(diǎn)與相位中心的距離。最終得到像素點(diǎn)P處的幅值

3.3 計(jì)算量分析

層析BP 算法逐個(gè)平面進(jìn)行插值和相位補(bǔ)償，成像結(jié)果準(zhǔn)確且立體，但在大場(chǎng)景下高度向Z的采樣點(diǎn)數(shù)可能很大，因此需要操作的點(diǎn)數(shù)龐大，計(jì)算效率較低。本文提出的地形表面有效點(diǎn)成像實(shí)現(xiàn)了降維成像，減少了像素點(diǎn)個(gè)數(shù)和插值及補(bǔ)償?shù)拇螖?shù)，能夠?qū)崿F(xiàn)更高效的成像。假設(shè)算法中一個(gè)像素點(diǎn)的插值和補(bǔ)償操作耗時(shí)t1，一次有效點(diǎn)DEM 提取操作耗時(shí)t2。

（1）層析BP 算法：設(shè)定X，Y，Z方向像素點(diǎn)個(gè)數(shù)分別為Nx，Ny，Nz，算法耗時(shí)為

（2）地形有效點(diǎn)BP 算法：設(shè)定提取的有效點(diǎn)共有M個(gè)不同的高度，那么像素點(diǎn)也是分布在不同高度的，在M個(gè)高度共分布Nx×Ny個(gè)像素點(diǎn)，算法耗時(shí)為

可以得到層析BP算法和地形有效點(diǎn)BP算法的耗時(shí)比為

由于DEM 提取耗時(shí)t2非常少，即t2＜＜Nx×Ny×t1，因此耗時(shí)比可以近似為層析BP時(shí)成像區(qū)域高度向的網(wǎng)格點(diǎn)數(shù)Nz，在大場(chǎng)景下Nz通常較大。因此可以得出結(jié)論，在同一場(chǎng)景且成像網(wǎng)格間隔設(shè)定相同的條件下，地形表面有效點(diǎn)BP 算法效率提升明顯。

4 驗(yàn)證實(shí)驗(yàn)

4.1 成像結(jié)果分析

本節(jié)我們利用機(jī)載三維圖像數(shù)據(jù)進(jìn)行仿真實(shí)驗(yàn)，主要的實(shí)驗(yàn)參數(shù)如表1所示。

表1 仿真參數(shù)Tab.1 Simulation parameters

仿真實(shí)驗(yàn)采用的原始數(shù)據(jù)是929*938*141的三維矩陣，具有目標(biāo)點(diǎn)的三維坐標(biāo)以及散射系數(shù)。由于數(shù)據(jù)量太大不方便仿真實(shí)驗(yàn)，因此對(duì)數(shù)據(jù)進(jìn)行等間隔采樣，采樣后得到92*93*141 的三維矩陣。首先給出三維數(shù)據(jù)在中間高度的二維圖像。

利用層析BP算法獲取到目標(biāo)點(diǎn)的幅值矩陣后，提取出有效點(diǎn)的DEM數(shù)據(jù)，并根據(jù)有效點(diǎn)的在矩陣中的索引從散射系數(shù)矩陣中提取出對(duì)應(yīng)的散射系數(shù)值。以下是有效點(diǎn)DEM圖像和有效點(diǎn)RCS圖像。

可以從圖6和圖7（a）對(duì)比發(fā)現(xiàn)原始圖像和提取出的有效點(diǎn)DEM圖像基本一致，有效點(diǎn)的分布反映了整個(gè)場(chǎng)景目標(biāo)點(diǎn)的分布。

在三維數(shù)據(jù)圖像中進(jìn)行了高度向成像切片分析，選擇了目標(biāo)點(diǎn)高度變化比較明顯的位置Y=40 m 處，下面是Y=40 m 處目標(biāo)點(diǎn)的幅值與高度的關(guān)系以及Y=40 m處有效點(diǎn)的高度。

從圖8（a）可以看出在X=0 m 到35 m 處最大幅值交替分布在60 m 到70 m 高度，X=35 m 之后最大幅值穩(wěn)定在70 m 高度處。這和圖8（b）中Y=40 m 處的有效點(diǎn)DEM 分布相同，可以證明提取得到的目標(biāo)點(diǎn)DEM數(shù)據(jù)準(zhǔn)確。

接下來(lái)利用地形有效點(diǎn)BP算法在提取出的有效點(diǎn)曲面上進(jìn)行成像，并將有效點(diǎn)的散射系數(shù)附加在對(duì)應(yīng)有效點(diǎn)上生成回波。得到BP成像結(jié)果如下。

通過(guò)圖9可以發(fā)現(xiàn)當(dāng)回波附加上對(duì)應(yīng)的散射系數(shù)后，有效點(diǎn)曲面上的BP 成像結(jié)果和原始三維數(shù)據(jù)圖像一致，可以看出場(chǎng)景中目標(biāo)點(diǎn)的分布。在連續(xù)形變監(jiān)測(cè)時(shí)，根據(jù)前后兩張成像結(jié)果中像素點(diǎn)的幅值變化就可以分析出場(chǎng)景的形變信息。

圖10中可以看到有一個(gè)相對(duì)獨(dú)立的強(qiáng)散射點(diǎn)，該點(diǎn)周?chē)纳⑸潼c(diǎn)幅值較低，不會(huì)影響該點(diǎn)的旁瓣值，因此針對(duì)這個(gè)點(diǎn)進(jìn)行成像效果分析。

如圖11 所示，將實(shí)際場(chǎng)景成像圖按照5∶1 放大，可以得到場(chǎng)景中位于坐標(biāo)（521，527）處的獨(dú)立點(diǎn)目標(biāo)周?chē)鷧^(qū)域展示圖。再對(duì)該點(diǎn)的成像結(jié)果進(jìn)行一維成像分析，得到對(duì)應(yīng)的峰值旁瓣比和-3 dB處寬度，與理論值進(jìn)行比較。

圖12是該點(diǎn)方位向的剖面圖，可以從圖中分析方位向成像指標(biāo)。

經(jīng)過(guò)分析后得到實(shí)際峰值旁瓣比為-8.08，-3 dB處寬度為27 m，與理論值基本一致。

如圖13所示，獨(dú)立點(diǎn)目標(biāo)在距離向上的一維成像結(jié)果是sinc 函數(shù)，符合圓周掃描地基SAR 的成像模型原理。

綜上，本文中對(duì)于原始三維機(jī)載圖像數(shù)據(jù)的成像仿真實(shí)驗(yàn)達(dá)到了預(yù)期的效果，且方位向和距離向的成像指標(biāo)都符合成像模型原理。

5 結(jié)論

本文提出了一種針對(duì)GBCSAR 的地形表面成像算法，通過(guò)將場(chǎng)景三維目標(biāo)點(diǎn)提取為三維曲面，并在該曲面上投影成像，實(shí)現(xiàn)了地形表面成像。本文利用BP 算法的優(yōu)勢(shì)避免了由于GBCSAR 弧形孔徑帶來(lái)的軌跡誤差，并且在曲面投影時(shí)省去了平面坐標(biāo)轉(zhuǎn)換等復(fù)雜操作。此外，與傳統(tǒng)的層析后向投影算法相比，地形表面后向投影算法通過(guò)降維的方式減小了計(jì)算量，提高了成像效率，這使得GBCSAR在地形形變連續(xù)監(jiān)測(cè)領(lǐng)域具有優(yōu)勢(shì)。