星載全方位SAR聯(lián)合圖像配準(zhǔn)和DEM提取方法

張金強(qiáng), 秦 強(qiáng), 宮 俊, 劉 偉, 蘇皎陽(yáng)

(1.上海無(wú)線電設(shè)備研究所,上海201109;2.上海目標(biāo)識(shí)別與環(huán)境感知工程技術(shù)研究中心,上海201109)

0 引言

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

1 聯(lián)合配準(zhǔn)和DEM提取方法

1.1 星載全方位SAR幾何構(gòu)型

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

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

圖1 星載全方位SAR幾何構(gòu)型

1.2 多角度SAR圖像序列聯(lián)合配準(zhǔn)方法

利用多角度SAR圖像序列提取觀測(cè)場(chǎng)景DEM時(shí),圖像對(duì)間方位夾角越大,配準(zhǔn)誤差導(dǎo)致的高程提取誤差越小。但是,圖像對(duì)間方位夾角越大,其相關(guān)性越低,配準(zhǔn)偏移量估計(jì)精度越低[5]。針對(duì)該問(wèn)題,基于星載全方位SAR獲取的多角度SAR圖像序列方位角連續(xù)這一特性,提出了多角度SAR圖像序列聯(lián)合配準(zhǔn)方法。該方法利用小方位夾角圖像對(duì)之間的配準(zhǔn)偏移量,估計(jì)大方位夾角圖像對(duì)之間的配準(zhǔn)偏移量。

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

則稱(chēng)之為小方位夾角圖像對(duì),反之則稱(chēng)之為大方位夾角圖像對(duì)。CΔfD需要根據(jù)具體的圖像分辨率、信噪比以及觀測(cè)場(chǎng)景內(nèi)地形起伏和地物類(lèi)型等參數(shù)確定。小方位夾角圖像對(duì)數(shù)量用K表示,并且K≥2N。

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

圖2 小方位夾角圖像對(duì)連接示意圖

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

式中:a為方位向坐標(biāo)向量;r為距離向坐標(biāo)向量;T表示矢量或矩陣的轉(zhuǎn)置操作。在對(duì)小方位夾角圖像對(duì)進(jìn)行配準(zhǔn)處理時(shí),需要指定主輔圖像,其編號(hào)可以用向量形式表示為

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

K個(gè)小方位夾角圖像對(duì)間方位向和距離向的配準(zhǔn)偏移量可以用向量表示為

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

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

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

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

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

利用該方法可以同時(shí)估計(jì)得到所有圖像相對(duì)于參考圖像的配準(zhǔn)偏移量。

1.3 聯(lián)合圖像配準(zhǔn)和DEM提取方法

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

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

式中:tn0表示圖像#n方位向第一個(gè)像素對(duì)應(yīng)的聚焦多普勒時(shí)刻;fp表示脈沖重復(fù)頻率;Rn0表示圖像#n距離向第一個(gè)像素對(duì)應(yīng)的雷達(dá)斜距;fs表示脈沖采樣頻率;c表示光速。

利用幾何法計(jì)算圖像對(duì)間配準(zhǔn)偏移量時(shí),需要輸入目標(biāo)高程h。若輸入先驗(yàn)?zāi)繕?biāo)高程存在誤差,將導(dǎo)致幾何法計(jì)算得到的配準(zhǔn)偏移量存在誤差。下面分析幾何法計(jì)算得到的配準(zhǔn)偏移量誤差隨輸入先驗(yàn)?zāi)繕?biāo)高程誤差的變化。

根據(jù)正向定位模型,將式(19)～式(21)中目標(biāo)三維位置矢量對(duì)目標(biāo)高程求一階偏導(dǎo)數(shù),得到

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

其中

式中:Pn、Vn和An分別表示T對(duì)應(yīng)的圖像#n的聚焦多普勒時(shí)刻的雷達(dá)APC位置矢量、速度矢量和加速度矢量。

根據(jù)式(22)～式(28),得到目標(biāo)在圖像#n上的二維坐標(biāo)相對(duì)于目標(biāo)高程的一階偏導(dǎo)數(shù)為

由于星載SAR雷達(dá)斜距較大,在局部范圍內(nèi)可以認(rèn)為,目標(biāo)在圖像對(duì)間的二維偏移量隨其高程線性變化。

將其他各圖像相對(duì)于參考圖像的幾何配準(zhǔn)偏移量用向量形式表示為

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

其中

因此,幾何配準(zhǔn)時(shí)所用先驗(yàn)?zāi)繕?biāo)高程誤差可以利用最小二乘算法估計(jì)得到

將其與先驗(yàn)?zāi)繕?biāo)高程相加,即為目標(biāo)估計(jì)高程。先驗(yàn)?zāi)繕?biāo)高程可以設(shè)置為0m,通過(guò)迭代處理最終實(shí)現(xiàn)高精度DEM提取。

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

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

因此,可以利用最小二乘算法同時(shí)估計(jì)得到所有圖像相對(duì)于參考圖像的配準(zhǔn)偏移量和先驗(yàn)?zāi)繕?biāo)高程誤差

2 仿真數(shù)據(jù)處理及結(jié)果分析

2.1 多角度SAR圖像序列仿真

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

表1 多角度SAR圖像仿真參數(shù)

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

圖3 衛(wèi)星星下點(diǎn)軌跡及觀測(cè)場(chǎng)景位置

圖4 仿真觀測(cè)場(chǎng)景分形插值后SRTM DEM

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

2.2 仿真數(shù)據(jù)處理及結(jié)果分析

圖5 多角度SAR圖像

利用本文方法對(duì)仿真得到的四個(gè)航過(guò)SAR圖像序列分別進(jìn)行處理,分別選擇方位角為0°,176.97°,95.98°和266.11°的圖像作為參考圖像。選擇小方位夾角圖像對(duì)時(shí),令式(2)中的α等于3。利用相關(guān)函數(shù)法逐像素估計(jì)小方位夾角圖像對(duì)間的精配準(zhǔn)偏移量。逐像素估計(jì)大方位夾角圖像對(duì)間精配準(zhǔn)偏移量時(shí),剔除相關(guān)系數(shù)低于0.5的小方位夾角圖像對(duì)間的配準(zhǔn)偏移量。如果剩余小方位夾角圖像對(duì)數(shù)量少于原數(shù)量的60%或?qū)е率?37)中的觀測(cè)矩陣為非列滿秩矩陣,則將參考圖像該像素對(duì)應(yīng)的目標(biāo)高程值置為無(wú)效值。利用幾何法計(jì)算圖像對(duì)間配準(zhǔn)偏移量時(shí),在仿真DEM上添加均值為0 m、標(biāo)準(zhǔn)差為17 m的隨機(jī)高程誤差,將添加高程誤差后的DEM作為先驗(yàn)?zāi)繕?biāo)高程輸入值。

利用四個(gè)航過(guò)的多角度SAR圖像序列估計(jì)得到四組觀測(cè)場(chǎng)景目標(biāo)高程,圖6(a)～圖6(d)分別給出了四個(gè)航過(guò)目標(biāo)高程估計(jì)誤差的分布。高程估計(jì)誤差為所估計(jì)DEM與仿真DEM的差值。統(tǒng)計(jì)得到四個(gè)航過(guò)目標(biāo)高程估計(jì)誤差范圍分別為(-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;無(wú)效值像素?cái)?shù)量分別為1,0,0,0。由上述結(jié)果可知,基于星載全方位SAR多角度圖像序列,利用本文方法可以高精度提取觀測(cè)場(chǎng)景DEM。

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

圖6 高程誤差分布

3 結(jié)論

針對(duì)星載全方位SAR獲取的多角度圖像序列,提出了多角度SAR圖像序列配準(zhǔn)和DEM提取聯(lián)合處理方法。首先,給出了星載全方位SAR幾何構(gòu)型;然后,提出了多角度SAR圖像序列聯(lián)合配準(zhǔn)方法,利用小方位夾角圖像對(duì)間配準(zhǔn)偏移量估計(jì)大方位夾角圖像對(duì)間配準(zhǔn)偏移量,解決了大方位夾角圖像對(duì)由于相關(guān)性低而難以精確配準(zhǔn)的問(wèn)題;最后,通過(guò)構(gòu)建圖像對(duì)間配準(zhǔn)偏移量與目標(biāo)高程的關(guān)系,提出了多角度SAR圖像序列配準(zhǔn)和DEM提取聯(lián)合處理方法,實(shí)現(xiàn)配準(zhǔn)偏移量和目標(biāo)高程的同時(shí)估計(jì)。仿真數(shù)據(jù)處理結(jié)果驗(yàn)證了本文方法的有效性和精確性。