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

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

引 言

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

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

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

1 基于SIFT子圖像融合直角坐標FBP算法

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

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

算法主要包括以下四個步驟：

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

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

③對子圖像進行相位梯度自聚焦處理；

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

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

1.1 直角坐標系的FBP算法

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

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

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

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

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

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

1.2 自聚焦處理

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

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

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

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

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

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

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

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

1.3 基于SIFT的子孔徑融合

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

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

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

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

②對兩幅圖像檢測到的特征點進行描述，以形成特征描述符；

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

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

①圖像極值點提取

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

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

②精確定位極值點

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

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

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

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

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

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

③生成特征描述符

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

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

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

④ SIFT 特征匹配與圖像融合

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

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

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

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

2.1 點目標仿真結(jié)果

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

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

圖7 點目標仿真結(jié)果Fig.7 Focusing results with simulated data

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

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

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

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

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

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

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

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

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

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

圖11 子圖像拼接后的強點的分辨率Fig.11 The resolution of strong character point after sub-image fusion

3 結(jié)束語

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

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