星載全方位SAR聯合圖像配準和DEM提取方法

張金強, 秦 強, 宮 俊, 劉 偉, 蘇皎陽

(1.上海無線電設備研究所,上海201109;2.上海目標識別與環境感知工程技術研究中心,上海201109)

0 引言

星載合成孔徑雷達(Synthetic Aperture Radar,SAR)因其全天時、全天候的全球觀測能力,在軍事偵察、國民經濟建設和科學研究中得到了廣泛的應用[1]。為了實現目標的精細化探測、高精度定位和全方位描述,近年來人們提出了星載全方位SAR系統。SAR衛星按照四組不同軌道根數的軌跡繞同一地面場景飛行,通過控制雷達波束指向實現對目標360°全方位觀測[2-4]。每一次航過的雷達波束斜視角范圍在-45°～+45°之間,合成孔徑時間可達百秒量級。將全孔徑回波數據進行子孔徑劃分,對各子孔徑回波數據進行聚焦成像處理,獲取多角度SAR圖像序列。為了利用單顆SAR衛星單航過錄取數據實現地形數字高程模型(Digital Elevation Model,DEM)的高精度提取,本文基于星載全方位SAR獲取的多角度圖像序列,提出了多角度SAR圖像序列配準和DEM提取聯合處理方法,并進行了仿真驗證。

1 聯合配準和DEM提取方法

1.1 星載全方位SAR幾何構型

星載全方位SAR幾何構型如圖1所示。圖中所有參數均定義在地心地固坐標系中。雷達天線相位中心(Antenna Phase Center,APC)位置矢量和速度矢量分別用P和V表示。P在地球表面的投影點稱為其星下點,用P′表示,P至P′的連線是地球表面的法線。雷達APC與目標最接近的時刻稱為零多普勒時刻。零多普勒時刻雷達APC所在位置稱為最近位置,用P0表示,P0的星下點用P′0表示。T表示地面目標位置矢量,P到T的矢量R稱為斜距矢量,P0到T的矢量R0稱為最短斜距矢量。垂直于雷達APC速度矢量并且包含雷達APC的平面稱為零多普勒面。R與零多普勒面之間的夾角稱為斜視角,用φ表示。當R與V的夾角小于90°時,φ的值為正;反之,φ的值為負。

為了描述不同航過雷達APC位置間的關系,定義方位角和方位夾角這兩個參數。在地球表面建立如圖1(b)所示的二維直角坐標系,其中坐標原點o為觀測場景中心目標所在位置,x軸和y軸相互垂直。o至P′的矢量與x軸的夾角稱為方位角,用φ表示。o到任意兩個雷達APC位置的星下點P′i和P′j的矢量之間的夾角稱為方位夾角,用Δφij表示。

圖1 星載全方位SAR幾何構型

1.2 多角度SAR圖像序列聯合配準方法

利用多角度SAR圖像序列提取觀測場景DEM時,圖像對間方位夾角越大,配準誤差導致的高程提取誤差越小。但是,圖像對間方位夾角越大,其相關性越低,配準偏移量估計精度越低[5]。針對該問題,基于星載全方位SAR獲取的多角度SAR圖像序列方位角連續這一特性,提出了多角度SAR圖像序列聯合配準方法。該方法利用小方位夾角圖像對之間的配準偏移量,估計大方位夾角圖像對之間的配準偏移量。

不失一般性,假設所有圖像均聚焦成像于斜平面,聚焦多普勒頻率等于子孔徑回波數據的多普勒中心頻率,并且聚焦多普勒頻率均勻分布。設聚焦多普勒頻率間隔為ΔfD,0,圖像數量為2N+1,圖像編號用向量表示為[0,1,2,…,2N],則多角度圖像序列的聚焦多普勒頻率可以用向量表示為[-N,-(N-1),…,0,…,(N-1),N]×ΔfD,0。根據圖像對間聚焦多普勒頻率差值的不同(即圖像對間方位夾角的不同),將其分為小方位夾角圖像對和大方位夾角圖像對。對于聚焦多普勒頻率為fD,i的圖像#i和聚焦多普勒頻率為fD,j的圖像#j,如果兩者的聚焦多普勒頻率差與判別門限CΔfD之間的關系滿足

則稱之為小方位夾角圖像對,反之則稱之為大方位夾角圖像對。CΔfD需要根據具體的圖像分辨率、信噪比以及觀測場景內地形起伏和地物類型等參數確定。小方位夾角圖像對數量用K表示,并且K≥2N。

圖2給出了N=3,CΔfD=2ΔfD,0和CΔfD=3ΔfD,0時,小方位夾角圖像對連接示意圖及其數量。隨著CΔfD的增大,所選擇的小方位夾角圖像對的數量增加。小方位夾角圖像對間的配準偏移量可以利用相關函數法等傳統SAR圖像對配準方法進行估計[6]。

圖2 小方位夾角圖像對連接示意圖

下面根據小方位夾角圖像對間配準偏移量與大方位夾角圖像對間配準偏移量之間的關系構建線性方程組。不失一般性,選擇聚焦多普勒頻率等于0Hz的SAR圖像作為參考圖像,參考圖像編號為#0。在參考圖像中選擇任意像素,設其坐標為(aref,rref),該像素對應的目標在其他2N幅圖像上的像素坐標可以用向量表示為

式中:a為方位向坐標向量;r為距離向坐標向量;T表示矢量或矩陣的轉置操作。在對小方位夾角圖像對進行配準處理時,需要指定主輔圖像,其編號可以用向量形式表示為

式中:M為主圖像編號向量;S為輔圖像編號向量;k=1,…,K,表示小方位夾角圖像對編號。

K個小方位夾角圖像對間方位向和距離向的配準偏移量可以用向量表示為

式中:δa為方位向配準偏移量;δr為距離向配準偏移量。每一小方位夾角圖像對之間的方位向和距離向配準偏移量可以表示為

式(9)和式(10)定義了兩個方程組,每一個方程組包含K個等式、2N個未知數,可以用矩陣表示為

式中:Aa和Ar均為K×2N維的矩陣;aref和rref均為K×1維的向量。對于?k=1,…,K和?n=1,…,2N,Aa、Ar、aref和rref對應的元素可以表示為

式中:n表示矩陣Aa和Ar的列索引。

矩陣Aa和Ar均為稀疏矩陣,其元素值取決于所選擇的小方位夾角圖像對。由稀疏矩陣的特性以及K≥2N可知,Aa和Ar均為列滿秩矩陣。因此,其他各圖像相對于參考圖像的配準偏移量可以利用最小二乘算法[7]估計得到

利用該方法可以同時估計得到所有圖像相對于參考圖像的配準偏移量。

1.3 聯合圖像配準和DEM提取方法

所有圖像相對于參考圖像的配準偏移量,還可以根據SAR成像幾何,利用幾何法計算得到。幾何法流程包括兩步:首先,通過正向定位獲取參考圖像像素對應的目標三維位置矢量;然后,通過反向定位獲取上述目標在圖像#n上的坐標。針對參考圖像上某一像素(aref,rref),將成像參數和軌道數據代入距離-多普勒模型[8],求解目標三維位置矢量,即正向定位。距離-多普勒模型表達式為

利用目標三維位置矢量T、圖像#n成像參數和軌道數據,首先根據式(20)計算該目標對應的圖像#n的聚焦多普勒時刻tn,將不同時刻的Pref和Vref代入式(20)計算公式左右兩邊之間的差值,最小差值所對應的時刻即為tn;然后根據式(19)計算時刻tn目標到雷達的斜距Rn;最后計算該目標在圖像#n上的二維坐標(an,rn)。上述過程稱為反向定位。其中an和rn的計算式為

式中:tn0表示圖像#n方位向第一個像素對應的聚焦多普勒時刻;fp表示脈沖重復頻率;Rn0表示圖像#n距離向第一個像素對應的雷達斜距;fs表示脈沖采樣頻率;c表示光速。

利用幾何法計算圖像對間配準偏移量時,需要輸入目標高程h。若輸入先驗目標高程存在誤差,將導致幾何法計算得到的配準偏移量存在誤差。下面分析幾何法計算得到的配準偏移量誤差隨輸入先驗目標高程誤差的變化。

根據正向定位模型,將式(19)～式(21)中目標三維位置矢量對目標高程求一階偏導數,得到

根據反向定位模型,將式(22)和式(23)中的聚焦多普勒時刻tn和雷達斜距Rn分別對目標三維位置矢量T求一階偏導數,得到

其中

式中:Pn、Vn和An分別表示T對應的圖像#n的聚焦多普勒時刻的雷達APC位置矢量、速度矢量和加速度矢量。

根據式(22)～式(28),得到目標在圖像#n上的二維坐標相對于目標高程的一階偏導數為

由于星載SAR雷達斜距較大,在局部范圍內可以認為,目標在圖像對間的二維偏移量隨其高程線性變化。

將其他各圖像相對于參考圖像的幾何配準偏移量用向量形式表示為

將利用式(17)和式(18)估計得到的配準偏移量與幾何配準偏移量對比,可以得到

其中

因此,幾何配準時所用先驗目標高程誤差可以利用最小二乘算法估計得到

將其與先驗目標高程相加,即為目標估計高程。先驗目標高程可以設置為0m,通過迭代處理最終實現高精度DEM提取。

觀察式(11)、式(12)和式(33),可以利用一個線性方程組將其聯合

式中:Ea和Er均為2N階單位矩陣。

因此,可以利用最小二乘算法同時估計得到所有圖像相對于參考圖像的配準偏移量和先驗目標高程誤差

2 仿真數據處理及結果分析

2.1 多角度SAR圖像序列仿真

多角度SAR圖像序列仿真參數如表1所示。軌道數據利用STK軟件生成,四組軌道的星下點軌跡如圖3所示,圖中箭頭指向衛星的飛行方向。觀測場景地形在航天飛機雷達地形測繪任務(Shuttle Radar Topography Mission,SRTM)DEM庫[9]中提取,高程變化范圍為(1225～1537)m。利用分形插值方法[10]對其進行插值處理,地面格網間距設為0.5m,地物類型設為土壤或巖石,目標后向散射系數根據Ulaby模型[11]計算得到。插值后DEM如圖4所示。

表1 多角度SAR圖像仿真參數

航過一仿真得到57幅圖像,聚焦多普勒頻率在(-270～+290)kHz范圍內均勻分布;航過二仿真得到50幅圖像,聚焦多普勒頻率在(-220～+270)kHz范圍內均勻分布;航過三仿真得到50幅圖像,聚焦多普勒頻率在(-230～+260)kHz范圍內均勻分布;航過四仿真得到47幅圖像,聚焦多普勒頻率在(-230～+230)kHz范圍內均勻分布。四個航過圖像序列聚焦多普勒頻率間隔均為10kHz,圖像信噪比均為20dB。觀測場景中心目標對應的四個航過第一幅、最后一幅和斜視角為0°的圖像的聚焦多普勒時刻雷達APC位置的星下點,分別如圖3中黑色實點所示。在該仿真中,圖1(b)中所定義的坐標系的x軸,由觀測場景中心目標位置指向其航過一零多普勒時刻雷達APC位置的星下點,y軸方向與航過一雷達APC速度方向相反。

圖3 衛星星下點軌跡及觀測場景位置

圖4 仿真觀測場景分形插值后SRTM DEM

圖5給出了方位角φ分別為0°,95.98°,176.97°和266.11°的仿真圖像,圖像大小均為19 000像素×10 500像素。

2.2 仿真數據處理及結果分析

圖5 多角度SAR圖像

利用本文方法對仿真得到的四個航過SAR圖像序列分別進行處理,分別選擇方位角為0°,176.97°,95.98°和266.11°的圖像作為參考圖像。選擇小方位夾角圖像對時,令式(2)中的α等于3。利用相關函數法逐像素估計小方位夾角圖像對間的精配準偏移量。逐像素估計大方位夾角圖像對間精配準偏移量時,剔除相關系數低于0.5的小方位夾角圖像對間的配準偏移量。如果剩余小方位夾角圖像對數量少于原數量的60%或導致式(37)中的觀測矩陣為非列滿秩矩陣,則將參考圖像該像素對應的目標高程值置為無效值。利用幾何法計算圖像對間配準偏移量時,在仿真DEM上添加均值為0 m、標準差為17 m的隨機高程誤差,將添加高程誤差后的DEM作為先驗目標高程輸入值。

利用四個航過的多角度SAR圖像序列估計得到四組觀測場景目標高程,圖6(a)～圖6(d)分別給出了四個航過目標高程估計誤差的分布。高程估計誤差為所估計DEM與仿真DEM的差值。統計得到四個航過目標高程估計誤差范圍分別為(-17.697 8～+16.020 3)m,(-12.329 7～+12.337 3)m,(-15.730 8～+15.337 4)m,(-12.320 8～+14.504 4)m;均方根值分別為1.891 1,1.604 0,1.387 4,1.327 5 m;無效值像素數量分別為1,0,0,0。由上述結果可知,基于星載全方位SAR多角度圖像序列,利用本文方法可以高精度提取觀測場景DEM。

由于雷達波束照射方向不同,不同角度SAR圖像中疊掩區域或低信噪比區域不同,導致利用不同航過多角度SAR圖像序列所提取的目標高程的無效值像素以及估計誤差大小的分布不同。為了提高DEM提取精度、降低無效值像素的數量,基于相關系數的高低將四個航過獲取的DEM進行融合處理。圖6(e)給出了融合后目標高程估計誤差的分布。融合后目標高程估計誤差范圍為(-4.199 8～+3.699 7)m,均方根值為0.443 4 m,無效值像素數為0。由上述結果可知,通過融合基于不同角度SAR圖像序列獲取的DEM,可以提高觀測場景DEM提取精度,同時降低無效值像素的數量。

圖6 高程誤差分布

3 結論

針對星載全方位SAR獲取的多角度圖像序列,提出了多角度SAR圖像序列配準和DEM提取聯合處理方法。首先,給出了星載全方位SAR幾何構型;然后,提出了多角度SAR圖像序列聯合配準方法,利用小方位夾角圖像對間配準偏移量估計大方位夾角圖像對間配準偏移量,解決了大方位夾角圖像對由于相關性低而難以精確配準的問題;最后,通過構建圖像對間配準偏移量與目標高程的關系,提出了多角度SAR圖像序列配準和DEM提取聯合處理方法,實現配準偏移量和目標高程的同時估計。仿真數據處理結果驗證了本文方法的有效性和精確性。