張佳琪，李涼海，張振華，左紹山

（北京遙測技術(shù)研究所 北京 100076）

引 言

合成孔徑雷達(dá)SAR（Synthetic Aperture Radar）具有全天時(shí)、全天候、高分辨率、大覆蓋面積、成像不受天氣和光照強(qiáng)度等影響的優(yōu)點(diǎn)，廣泛應(yīng)用于軍事偵察、地球遙感、資源勘探、農(nóng)作物估產(chǎn)、自然災(zāi)害監(jiān)測、環(huán)境保護(hù)等諸多國防和國民經(jīng)濟(jì)的重要領(lǐng)域。針對(duì)SAR不同的成像模式，國內(nèi)外學(xué)者先后研究了各種不同的距離多普勒算法、時(shí)域成像算法和二維頻域成像算法[1,2]。

后向投影算法BPA（Back Projection Algorithm）是SAR時(shí)域成像算法的典型，適用于任意場景、任意雷達(dá)工作模式成像，并且具有很好的并行性，易于工程實(shí)現(xiàn)。傳統(tǒng)的BPA是基于極坐標(biāo)系的，能直接從回波數(shù)據(jù)搜索出目標(biāo)點(diǎn)的后向散射系數(shù)，對(duì)含有運(yùn)動(dòng)誤差的SAR成像數(shù)據(jù)進(jìn)行精確聚焦[3,4]，但是其需要逐點(diǎn)遍歷，計(jì)算量過大，效率低，不適用于SAR實(shí)時(shí)成像處理。后來有學(xué)者提出的快速后向投影算法FBPA（Fast Back Projection Algorithm）采用極坐標(biāo)系，通過對(duì)子孔徑進(jìn)行相干積累來提高運(yùn)算效率，但是由于不同極坐標(biāo)系之間的關(guān)系是非線性的，因此需要通過逐像素投影來實(shí)現(xiàn)子圖像融合，這就增加了計(jì)算量，并且在融合過程中為了提高融合效率，通常會(huì)采用一定的近似處理，導(dǎo)致成像分辨率、峰值旁瓣比等成像指標(biāo)下降。左紹山等人提出一種改進(jìn)的直角坐標(biāo)系下的FBP算法[5]，只需坐標(biāo)平移即可實(shí)現(xiàn)直角坐標(biāo)系下的子圖像融合，既減少了計(jì)算量，又避免了極坐標(biāo)下對(duì)融合的近似，提高了融合精度。同時(shí)該算法采用方位譜壓縮技術(shù)解決了文獻(xiàn)[6]所提到的局部直角坐標(biāo)系下子圖像的方位譜寬度要比局部極坐標(biāo)系下子圖像的方位譜寬度大從而造成子孔徑信號(hào)后向投影成像過程計(jì)算量大的問題。但是隨著分辨率的提高，天線相位中心與成像目標(biāo)之間的運(yùn)動(dòng)補(bǔ)償精度也越來越高，因此即使是裝有先進(jìn)慣性導(dǎo)航測量單元的SAR系統(tǒng)，要想精確補(bǔ)償載機(jī)的運(yùn)動(dòng)也是很困難的，另外由于大氣傳播效應(yīng)產(chǎn)生的相位誤差也會(huì)使SAR高分辨圖像散焦。而文獻(xiàn)[5]算法沒有對(duì)運(yùn)動(dòng)誤差進(jìn)行補(bǔ)償，因此不適用于實(shí)測數(shù)據(jù)。針對(duì)實(shí)測數(shù)據(jù)，張磊等人提出了一種多孔徑圖像偏移自聚焦算法與時(shí)域快速后向投影算法相結(jié)合的算法，該算法仍是基于極坐標(biāo)系的。選擇方位角正弦為子孔徑坐標(biāo)系的方位坐標(biāo)，后向投影成像的極坐標(biāo)方位軸對(duì)應(yīng)多普勒頻率軸，對(duì)各子孔徑極坐標(biāo)成像直接進(jìn)行逆傅里葉變換到方位時(shí)間域，在時(shí)間域?qū)崿F(xiàn)各個(gè)子孔徑數(shù)據(jù)的調(diào)頻率估計(jì)，繼而進(jìn)行相位誤差估計(jì)，實(shí)現(xiàn)快速后向投影成像過程中的自聚焦處理。

本文提出一種將直角坐標(biāo)系下的子孔徑FBP算法與基于尺度不變特性變換SIFT（Scale Invariant Feature Transform）的SAR圖像配準(zhǔn)方法相結(jié)合的算法，在減小計(jì)算量的同時(shí)，提高融合精度和算法的并行性，另外通過合理劃分子孔徑，達(dá)到實(shí)時(shí)處理的效果。本文將原始數(shù)據(jù)劃分子孔徑后，經(jīng)過后向投影得到子孔徑圖像是基于直角坐標(biāo)系的，在融合過程中無需將極坐標(biāo)與直角坐標(biāo)進(jìn)行轉(zhuǎn)換，亦無需先插值再融合，減小了計(jì)算量。然后通過特征點(diǎn)匹配完成子圖像融合，達(dá)到提高成像效率與分辨率的目的。相比于文獻(xiàn)[5]，本文提出的方法進(jìn)行了自聚焦處理，能夠處理實(shí)測數(shù)據(jù)；相比于文獻(xiàn)[7]，本文提出的方法采用了一種新的融合方式，能夠提高分辨率與融合效率。本文算法具有良好的并行性，為后續(xù)子圖像融合而采用直角坐標(biāo)FBP結(jié)合SAR實(shí)時(shí)成像算法的FPGA實(shí)現(xiàn)做準(zhǔn)備。

1 基于SIFT子圖像融合直角坐標(biāo)FBP算法

直角坐標(biāo)系下的FBP算法，是將整個(gè)合成孔徑劃分為若干個(gè)子孔徑，在局部直角坐標(biāo)系下以子孔徑中心為原點(diǎn)重建子孔徑圖像，每幅子圖像具有方位向低分辨和距離向全分辨的特點(diǎn)，在直角坐標(biāo)系下將全部子圖像進(jìn)行一次方位平移，得到全孔徑分辨率的SAR圖像[8]。本文采用的融合方法與傳統(tǒng)方法不同，本文基于SIFT特征點(diǎn)匹配的，即將子孔徑圖像進(jìn)行特征點(diǎn)匹配，根據(jù)SIFT提取出來的特征因子找到不同子孔徑圖像之間的轉(zhuǎn)換關(guān)系，利用轉(zhuǎn)換關(guān)系完成子圖像的融合，得到全孔徑圖像。本文算法的具體流程如圖1所示。

圖1 基于SIFT子圖像融合的直角坐標(biāo)FBP算法流程Fig.1 Based on SIFT sub image fusion for rectangular coordinates FBPA flow chart

算法主要包括以下四個(gè)步驟：

①將原始數(shù)據(jù)進(jìn)行距離向脈壓，并對(duì)脈壓后的數(shù)據(jù)劃分子孔徑；

②對(duì)每個(gè)子孔徑的距離脈壓結(jié)果進(jìn)行后向投影，得到各個(gè)子孔徑回波的子圖像；

③對(duì)子圖像進(jìn)行相位梯度自聚焦處理；

④對(duì)處理后的子圖像進(jìn)行SIFT特征點(diǎn)提取，將相鄰子圖像的特征描述符作匹配，完成子圖像匹配融合，得到全孔徑圖像。

本文算法是基于直角坐標(biāo)系的，在完成子圖像SIFT特征點(diǎn)提取、特征點(diǎn)匹配之后，亦無需插值，無需逐點(diǎn)操作，只需按照坐標(biāo)轉(zhuǎn)換關(guān)系進(jìn)行方位向的平移即可完成子圖像融合。本文算法在子圖像融合過程中采用上采樣后的復(fù)圖像作融合，相比于傳統(tǒng)的基于實(shí)圖像的融合，分辨率能夠提高一倍。下面對(duì)算法的具體實(shí)現(xiàn)過程進(jìn)行介紹。

1.1 直角坐標(biāo)系的FBP算法

對(duì)已知場景中任意點(diǎn)目標(biāo)N(x,y)的回波基頻信號(hào)作距離向脈沖壓縮，可以得到脈壓的結(jié)果為

其中，c為電磁波傳播速度，B為信號(hào)帶寬，aa(tm)為雷達(dá)方位窗函數(shù)，tm為方位慢時(shí)間，為距離快時(shí)間，λ為雷達(dá)發(fā)射信號(hào)波長，R(tm;x,y)為任意tm時(shí)刻點(diǎn)目標(biāo)(x,y)到雷達(dá)的距離。

接下來完成子孔徑劃分，再將每個(gè)子孔徑信號(hào)經(jīng)過后向投影算法得到對(duì)應(yīng)的子圖像，第i個(gè)子孔徑對(duì)應(yīng)子圖像的表達(dá)式為

式中，v為雷達(dá)平臺(tái)飛行速度，(xn,yn)表示子圖像中網(wǎng)格點(diǎn)的坐標(biāo)。將瞬時(shí)斜距在vtm=xn處進(jìn)行泰勒級(jí)數(shù)展開，并將展開結(jié)果代入式（2）可以得到

式（3）最后一個(gè)指數(shù)項(xiàng)會(huì)導(dǎo)致方位向頻譜展寬[9]。對(duì)于第i個(gè)子孔徑，取近似refnR≈y，Rref為參考距離，可以得到降低方位譜寬度函數(shù)為

為了降低方位譜寬度，可將式（3）與式（4）相乘。如此一來，方位譜寬度變窄，成像網(wǎng)格方位間距變大為λyn/Lsub。因此在子孔徑后向投影成像過程中，將像素方位間隔設(shè)置為小于成像網(wǎng)格方位間距，使得方位采樣點(diǎn)數(shù)減少，進(jìn)而降低了計(jì)算量，提高了子孔徑成像效率。

1.2 自聚焦處理

自聚焦的目的是將相位誤差補(bǔ)償?shù)胶侠淼南薅?，便于?duì)實(shí)測數(shù)據(jù)進(jìn)行聚集成像處理。目前常用的自聚焦算法包含子孔徑相關(guān)算法MD（Map Drift），相位差算法PD（Phase Difference），相位梯度自聚焦算法PGA（Phase Gradient Autofocus）等。PGA算法是目前研究的最成熟的算法且較MD算法來說更適用于估計(jì)高階相位誤差，并且由于本文子圖像僅僅是經(jīng)過了距離向脈沖壓縮而沒有對(duì)方位向進(jìn)行處理，所以可以對(duì)子圖像采用相位梯度自聚焦算法作相位補(bǔ)償，以滿足SAR成像高分辨率的要求。PGA方法基于最大似然估計(jì)理論進(jìn)行最優(yōu)估計(jì)，并且利用了多個(gè)距離單元上的相位誤差信息，能夠提高估計(jì)結(jié)果的穩(wěn)健性[10]。

PGA算法流程如圖2所示，從復(fù)圖像域出發(fā)，以迭代的方式，逐漸改進(jìn)圖像的聚集，主要包含以下四個(gè)步驟：

①循環(huán)移位：首先需要在每個(gè)距離單元上選取最強(qiáng)散射點(diǎn)，并將其循環(huán)移位到孔徑中心對(duì)應(yīng)的位置。循環(huán)移位能夠保留相位誤差對(duì)所選取散射點(diǎn)的影響，還能去除與目標(biāo)相關(guān)的線性相位分量。對(duì)于包含較少孤立散射點(diǎn)的一幅散焦的圖像來說，循環(huán)移位是把各距離單元上的強(qiáng)點(diǎn)進(jìn)行排列，從而有利于確定窗寬，而正確的窗寬選擇能夠改善相位估計(jì)信號(hào)的信噪比。

②加窗：對(duì)經(jīng)過循環(huán)移位后的數(shù)據(jù)加矩形窗，能夠保持最強(qiáng)散射點(diǎn)所包含的模糊信息，同時(shí)能夠減少來自背景雜波或其它鄰近目標(biāo)的影響，使得用來估計(jì)相位誤差的數(shù)據(jù)具有較高的信噪比。

③相位誤差估計(jì)：由于相位誤差估計(jì)是在距離壓縮方位歷史域進(jìn)行，因此經(jīng)循環(huán)移位和加窗處理后的圖像數(shù)據(jù)首先要進(jìn)行方位向傅里葉反變換。Jakowatz等利用特征矢量的方法，進(jìn)行相位誤差的最大似然估計(jì)，并且證明了估計(jì)方差達(dá)到克拉美-羅限[11]。對(duì)于第n個(gè)距離單元，當(dāng)只采用相鄰的兩個(gè)方位向數(shù)據(jù)進(jìn)行估計(jì)時(shí)，相位差為

其中g(shù)*n(m)為距離壓縮相位時(shí)域數(shù)據(jù)，m=1,2,3,…,M，M為方位向脈沖數(shù)，n=1,2,3,…,N，N為距離向脈沖數(shù)。通過累加可算出全部孔徑上的相位誤差，即

④迭代相位校正：將距離壓縮方位歷史域的數(shù)據(jù)與相位誤差估計(jì)進(jìn)行共軛相乘，經(jīng)過傅立葉變換到圖像域，即可改善圖像聚焦。這一估計(jì)和校正過程需要重復(fù)迭代四至六次，以逐漸提高圖像的質(zhì)量[12]。

圖2 PGA流程Fig.2 The PGA work flow

1.3 基于SIFT的子孔徑融合

下面以兩個(gè)子圖像融合為例來說明基于SIFT圖像匹配算法進(jìn)行圖像融合的基本流程。本文采用的SIFT方法與傳統(tǒng)方法不同，它首先利用實(shí)圖像獲得特征點(diǎn)與融合移動(dòng)量，然后用復(fù)圖像進(jìn)行融合，其流程如圖3所示，可以概括為如下三個(gè)過程[13]：

圖3 SIFT匹配算法Fig.3 SIFT matching algorithm framework

圖4 構(gòu)造DOG尺度空間Fig.4 Construction of DOG scale space

①提取實(shí)數(shù)子圖像中的尺度不變特征點(diǎn)；

②對(duì)兩幅圖像檢測到的特征點(diǎn)進(jìn)行描述，以形成特征描述符；

③將實(shí)圖像1與實(shí)圖像2的SIFT特征描述符進(jìn)行匹配，對(duì)復(fù)數(shù)子圖像1和2按照匹配好的特征點(diǎn)進(jìn)行平移，完成子圖像的復(fù)值融合。

與傳統(tǒng)的快速后向投影算法的融合方式相比，本文基于SIFT的融合采用的是子圖像在時(shí)域的復(fù)數(shù)結(jié)果直接進(jìn)行累加，能夠提高分辨率，具體實(shí)現(xiàn)過程如下：

①圖像極值點(diǎn)提取

首先利用高斯函數(shù)G(x,y,σ)獲取子圖像GH1(x,y)的尺度空間L(x,y,σ)，其中σ為方差。具體過程為

式（8）是均值為0、方差為σ的標(biāo)準(zhǔn)二維高斯函數(shù)。高斯卷積可以得到相鄰尺度空間的圖像，然后對(duì)其進(jìn)行差分，得到高斯差分尺度空間DOG（Difference of Gaussian scale-space），過程如圖4所示。尺度空間極值檢測過程如圖5所示，在高斯差分尺度空間中選定一個(gè)圖像中的一個(gè)采樣點(diǎn)與其上下各一幅相鄰圖像中的8個(gè)鄰域點(diǎn)進(jìn)行比較，判斷其是否為極值點(diǎn)，如是則設(shè)置為粗選極值點(diǎn)。

②精確定位極值點(diǎn)

在獲得了所有的粗選極值點(diǎn)之后，還需要進(jìn)一步篩選，刪掉那些對(duì)比度較低的點(diǎn)，以提高算法匹配的穩(wěn)定性、增強(qiáng)最終極值點(diǎn)的抗噪性，具體通過擬合三維二次函數(shù)來實(shí)現(xiàn)。

將尺度空間函數(shù)D(x,y,σ)在待定點(diǎn)X處，進(jìn)行泰勒級(jí)數(shù)展開，得到

其中X=(x,y,σ)T。通過式（9）對(duì)X求導(dǎo)結(jié)果取0，可以得到極值點(diǎn)的偏移量為

若式（10）的極值點(diǎn)偏移量x與y中有一個(gè)數(shù)值超過0.5，可以判斷極值點(diǎn)偏向另一端的待選點(diǎn)。此時(shí)，再以該點(diǎn)為待選點(diǎn)重復(fù)上述過程，直到有較精確的點(diǎn)，即獲得最后的，該較精確的點(diǎn)與插值中心點(diǎn)的距離在任意一個(gè)方向上的偏差小于其到相鄰點(diǎn)的距離。如此得到的特征點(diǎn)的精度可達(dá)到亞像素級(jí)。

圖5 尺度空間極值檢測過程Fig.5 Scale space extreme detection

圖6 SIFT算法融合過程Fig.6 SIFT algorithm mix process

③生成特征描述符

在通過前面的處理獲得圖像特征算子的位置信息之后，能夠得到尺度不變特征算子。此外，由于每個(gè)特征點(diǎn)的鄰域像素內(nèi)梯度方向分布的不同，通過計(jì)算每個(gè)特征點(diǎn)的方向參數(shù)，使得特征算子具備旋轉(zhuǎn)不變性的特點(diǎn)。

根據(jù)特征點(diǎn)位置選擇離它最近的高斯平滑圖像L，對(duì)其中的每個(gè)采樣點(diǎn)L(x,y)，在其四鄰域內(nèi)利用下式計(jì)算獲取該點(diǎn)的梯度幅值m(x,y)和方向角θ(x,y)。

在實(shí)際中，為了提高圖像匹配的穩(wěn)健性，可以對(duì)每個(gè)特征點(diǎn)采用4×4個(gè)種子點(diǎn)進(jìn)行描述，得到128維SIFT特征描述符。

④ SIFT 特征匹配與圖像融合

在獲取圖像的特征描述符之后，還要對(duì)待匹配的兩幅圖像間的SIFT特征描述符進(jìn)行相似性度量，本文采用經(jīng)典的SIFT方法，即采用歐式距離進(jìn)行度量。

特征匹配完成后，需要對(duì)提取出的匹配特征點(diǎn)對(duì)進(jìn)行最小二乘算法擬合，得到最優(yōu)相似參數(shù)包括：平移量、尺度的伸縮變化量、旋轉(zhuǎn)角度，然后得到變化矩陣。以相鄰兩幅子圖像中特征匹配度高的區(qū)域?yàn)榛鶞?zhǔn)，完成子圖像的融合拼接，直至融合完所有的子圖像，得到一副拼接好的全孔徑圖像。由于本文采用直角坐標(biāo)系的FBP算法，能夠根據(jù)匹配特征點(diǎn)對(duì)直接進(jìn)行融合而無需坐標(biāo)轉(zhuǎn)換或插值，提高了融合效率。圖6舉例說明了采用SIFT算法的融合過程。

2 數(shù)據(jù)處理結(jié)果

為了驗(yàn)證本文所提出方法的有效性，本文首先采用條帶模式下點(diǎn)目標(biāo)仿真的方式，驗(yàn)證改進(jìn)后的FBP算法的有效性，然后結(jié)合機(jī)載實(shí)測數(shù)據(jù)驗(yàn)證基于SIFT子圖像融合的直角坐標(biāo)系下的FBP算法在實(shí)測數(shù)據(jù)成像中的有效性。以下是實(shí)驗(yàn)結(jié)果。

2.1 點(diǎn)目標(biāo)仿真結(jié)果

以條帶模式為例，點(diǎn)目標(biāo)驗(yàn)證結(jié)果如下。本文仿真了一個(gè)3×5的點(diǎn)陣目標(biāo)，仿真參數(shù)如表1所示，最終的成像結(jié)果如圖7所示，可以看到各點(diǎn)目標(biāo)均能得到良好聚集。

表1 條帶模式下的點(diǎn)目標(biāo)仿真參數(shù)Table 1 Stripmap SAR point target simulation parameters

圖7 點(diǎn)目標(biāo)仿真結(jié)果Fig.7 Focusing results with simulated data

為了定量分析，本文將左上角點(diǎn)目標(biāo)A的特性進(jìn)行具體分析如圖8所示，計(jì)算得到點(diǎn)A距離向和方位向的峰值旁瓣比分別為-13.44dB和-13.19dB。參考文獻(xiàn)[14]的評(píng)估方法，方位向的峰值旁瓣比能達(dá)到-13dB以上，因此本算法的成像結(jié)果在方位向具有良好的聚集效果。

圖8 點(diǎn)目標(biāo)A的距離和方位包絡(luò)Fig.8 Point spread functions of target point A

2.2 實(shí)測數(shù)據(jù)處理結(jié)果

實(shí)測數(shù)據(jù)采用的是某無人機(jī)SAR錄取的回波數(shù)據(jù)，原始回波數(shù)據(jù)大小為900×8192（距離×方位），劃分為8個(gè)子孔徑，每個(gè)子孔徑包含1024個(gè)脈沖，將原始回波數(shù)據(jù)進(jìn)行慣導(dǎo)補(bǔ)償與運(yùn)動(dòng)補(bǔ)償后采用本文的成像算法得到拼接前子孔徑圖像。本小節(jié)以第3個(gè)、第4個(gè)子孔徑融合為例說明采用SIFT方法的融合過程與效果。

圖9 子孔徑3、子孔徑4成像結(jié)果Fig.9 The results of sub-aperture 3 and sub-aperture 4

第3個(gè)、第4個(gè)子圖像的結(jié)果分別如圖9（a）、圖9（b）所示。對(duì)其中包含強(qiáng)點(diǎn)的區(qū)域放大可以得到圖9（c）、圖9（d）的結(jié)果，根據(jù)SIFT得到的坐標(biāo)關(guān)系，選取的強(qiáng)點(diǎn)如圖9（c）、圖9（d）所示，兩幅圖像的兩個(gè)強(qiáng)點(diǎn)具有相同的距離向坐標(biāo)和不同的方位向坐標(biāo)。分析這兩個(gè)強(qiáng)點(diǎn)的3dB方位分辨率特性如圖9（e）、圖9（f）所示，融合前第三個(gè)子圖像分辨率為2.2203m，融合前第四個(gè)子圖像分辨率為1.7128m。

完成SIFT子圖像拼接的結(jié)果如圖10所示。圖10（a）、圖10（b）分別表示SIFT融合前第3個(gè)、第4個(gè)子孔徑圖像，圖10（c）、圖10（d）分別為子圖像3、子圖像4的特征描述符，圖10（e）為子圖像3、子圖像4的SIFT匹配特征點(diǎn)對(duì)，圖10（f）為子圖像3、子圖像4融合后的圖像。

圖10 子圖像拼接過程Fig.10 The process of the sub-aperture mosaic

仍然選取圖9（c）位置的強(qiáng)點(diǎn)并分析其3dB的方位向分辨率，結(jié)果如圖11所示，該強(qiáng)點(diǎn)處分辨率為1.1419m，與圖9（e）、圖9（f）相比，分辨率均有提高。

表2為采用SIFT方法找到的每一幅子圖像的特征點(diǎn)。將全部子孔徑進(jìn)行拼接，可以得到拼接后的全孔徑圖像如圖10（g）所示。從圖10（g）可以看出，拼接結(jié)果良好，圖像無錯(cuò)位現(xiàn)象產(chǎn)生也沒有拼縫產(chǎn)生，而且采用SIFT融合比傳統(tǒng)的融合辦法得到的圖像分辨率高，實(shí)測數(shù)據(jù)處理結(jié)果驗(yàn)證了本文算法的有效性。

圖11 子圖像拼接后的強(qiáng)點(diǎn)的分辨率Fig.11 The resolution of strong character point after sub-image fusion

3 結(jié)束語

針對(duì)SAR實(shí)測數(shù)據(jù)計(jì)算量大的問題，本文提出了一種基于SIFT子圖像融合的SAR子孔徑成像算法，將基于尺度不變特性的融合方法與直角坐標(biāo)系下的快速后向投影算法相結(jié)合，解決了傳統(tǒng)后向投影算法計(jì)算量大的問題，另外進(jìn)行子孔徑劃分可以提高算法的并行性，為后期FPGA的實(shí)時(shí)SAR成像處理提供方便。本文采用的SIFT融合方法是對(duì)復(fù)數(shù)圖像的融合，能夠提高全孔徑圖像的分辨率。數(shù)據(jù)處理結(jié)果表明，該方法能夠提高SAR成像分辨率，并且很好地解決了BP算法計(jì)算量大的問題。后期研究將圍繞SAR子孔徑成像算法及其FPGA實(shí)現(xiàn)問題展開。

表2 各子圖像特征點(diǎn)數(shù)目Table 2 Number of every sub image features points