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

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

0 引言

星載合成孔徑雷達(dá)(Synthetic Aperture Radar,SAR)因其全天時(shí)、全天候的全球觀測能力,在軍事偵察、國民經(jīng)濟(jì)建設(shè)和科學(xué)研究中得到了廣泛的應(yīng)用[1]。為了實(shí)現(xiàn)目標(biāo)的精細(xì)化探測、高精度定位和全方位描述,近年來人們提出了星載全方位SAR系統(tǒng)。SAR衛(wèi)星按照四組不同軌道根數(shù)的軌跡繞同一地面場景飛行,通過控制雷達(dá)波束指向?qū)崿F(xiàn)對目標(biāo)360°全方位觀測[2-4]。每一次航過的雷達(dá)波束斜視角范圍在-45°～+45°之間,合成孔徑時(shí)間可達(dá)百秒量級。將全孔徑回波數(shù)據(jù)進(jìn)行子孔徑劃分,對各子孔徑回波數(shù)據(jù)進(jìn)行聚焦成像處理,獲取多角度SAR圖像序列。為了利用單顆SAR衛(wèi)星單航過錄取數(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)稱為其星下點(diǎn),用P′表示,P至P′的連線是地球表面的法線。雷達(dá)APC與目標(biāo)最接近的時(shí)刻稱為零多普勒時(shí)刻。零多普勒時(shí)刻雷達(dá)APC所在位置稱為最近位置,用P0表示,P0的星下點(diǎn)用P′0表示。T表示地面目標(biāo)位置矢量,P到T的矢量R稱為斜距矢量,P0到T的矢量R0稱為最短斜距矢量。垂直于雷達(dá)APC速度矢量并且包含雷達(dá)APC的平面稱為零多普勒面。R與零多普勒面之間的夾角稱為斜視角,用φ表示。當(dāng)R與V的夾角小于90°時(shí),φ的值為正;反之,φ的值為負(fù)。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

利用幾何法計(jì)算圖像對間配準(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)三維位置矢量對目標(biāo)高程求一階偏導(dǎo)數(shù),得到

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

其中

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

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

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

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

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

其中

因此,幾何配準(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,通過迭代處理最終實(shí)現(xiàn)高精度DEM提取。

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

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

因此,可以利用最小二乘算法同時(shí)估計(jì)得到所有圖像相對于參考圖像的配準(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)星的飛行方向。觀測場景地形在航天飛機(jī)雷達(dá)地形測繪任務(wù)(Shuttle Radar Topography Mission,SRTM)DEM庫[9]中提取,高程變化范圍為(1225～1537)m。利用分形插值方法[10]對其進(jìn)行插值處理,地面格網(wǎng)間距設(shè)為0.5m,地物類型設(shè)為土壤或巖石,目標(biāo)后向散射系數(shù)根據(jù)Ulaby模型[11]計(jì)算得到。插值后DEM如圖4所示。

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

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

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

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

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

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

圖5 多角度SAR圖像

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

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

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

圖6 高程誤差分布

3 結(jié)論

針對星載全方位SAR獲取的多角度圖像序列,提出了多角度SAR圖像序列配準(zhǔn)和DEM提取聯(lián)合處理方法。首先,給出了星載全方位SAR幾何構(gòu)型;然后,提出了多角度SAR圖像序列聯(lián)合配準(zhǔn)方法,利用小方位夾角圖像對間配準(zhǔn)偏移量估計(jì)大方位夾角圖像對間配準(zhǔn)偏移量,解決了大方位夾角圖像對由于相關(guān)性低而難以精確配準(zhǔn)的問題;最后,通過構(gòu)建圖像對間配準(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)證了本文方法的有效性和精確性。