星載InSAR方位向逐行的順軌有效基線計(jì)算方法?

劉昆池, 苗洪利, 楊忠昊, 張佳輝

(中國(guó)海洋大學(xué)信息科學(xué)與工程學(xué)部, 山東 青島 266100)

海流的研究對(duì)海軍作戰(zhàn)、海上救援、海岸帶建設(shè)、漁業(yè)、船舶航運(yùn)路線的制定、海洋環(huán)境的監(jiān)測(cè)和預(yù)報(bào)等都有重要意義[1-2]。海流計(jì)[3]、浮標(biāo)、地波雷達(dá)[4-5]等難以實(shí)現(xiàn)全球范圍的觀測(cè),星載干涉合成孔徑雷達(dá)(Interferometric synthetic aperture radar,InSAR)是一種微波遙感技術(shù),具有全天時(shí)、全天候、數(shù)據(jù)處理精度高等優(yōu)點(diǎn)[6]。InSAR分為單星重復(fù)軌道、單星雙天線及雙星編隊(duì)伴飛等方式[7]。利用InSAR技術(shù)的順軌干涉測(cè)量(Along track interferometry,ATI)可以進(jìn)行地面運(yùn)動(dòng)目標(biāo)的速度監(jiān)測(cè),當(dāng)然也可應(yīng)用于海洋領(lǐng)域,如海表面流速、海面風(fēng)速、冰川運(yùn)動(dòng)及艦船速度等方面的觀測(cè)[8-11]。ATI反演海表面流速相比于單景SAR多普勒中心偏移法(DCA)[12-13]具有更高的時(shí)空分辨率和反演精度。而順軌干涉基線的準(zhǔn)確估算在InSAR數(shù)據(jù)影像配準(zhǔn)、平地相位去除[14]、地形相位計(jì)算等方面具有重要地位,是反演高精度海表面流場(chǎng)的關(guān)鍵參數(shù)[15]。目前,常用的方法是從衛(wèi)星影像的頭文件中直接讀取或基于影像相位中心計(jì)算獲得順軌有效基線長(zhǎng)度。在反演海流方面,使用單一值反演整幅影像對(duì)應(yīng)的徑向流速會(huì)有很大的誤差。本文在進(jìn)行軌道擬合計(jì)算基礎(chǔ)上,建立順軌有效基線值與方位向行號(hào)的函數(shù)關(guān)系式,得到方位向每一行的順軌有效基線長(zhǎng)度,以實(shí)現(xiàn)在海流反演方面的順軌有效基線的精細(xì)化使用,將會(huì)帶來反演精度的提高。

1 順軌干涉有效基線及海流速度反演

InSAR測(cè)量中的基線是指兩個(gè)天線對(duì)同一地面目標(biāo)觀測(cè)時(shí),主輔天線相位中心(Antenna phase centre, APC)的連線,也稱為物理基線[16]。如果是雙發(fā)雙收的模式,有效基線長(zhǎng)度等于物理基線;如果是單發(fā)雙收,有效基線長(zhǎng)度為物理基線的一半[10]。有效基線在方位向上的投影為順軌干涉有效基線,在距離向上的投影為交軌干涉有效基線。

ATI-SAR反演海流方法的基本原理為[8]:在兩個(gè)天線均工作在條帶模式下,在前后掃過海面上同一個(gè)區(qū)域時(shí),后天線到達(dá)前天線照射的同一海面區(qū)域的時(shí)間差ΔT為:

(1)

式中:V為衛(wèi)星平臺(tái)的速度;Be為ATI的順軌有效基線長(zhǎng)度。當(dāng)被前天線照射的區(qū)域再被后天線照射時(shí)徑向(信號(hào)視線方向)移動(dòng)的速度ULOS為:

(2)

式中:ΔR是前后兩個(gè)天線通過被照射區(qū)域所接收到的信號(hào)的路程差。

對(duì)于單發(fā)雙收模式:

(3)

式中:λ為信號(hào)波長(zhǎng);φATI為兩個(gè)天線所接收到信號(hào)的相位差,正是ATI-SAR所要獲得的兩幅圖像的干涉相位數(shù)據(jù)。

ULOS不但包括由風(fēng)、浪等因素造成的海流徑向速度,還包括長(zhǎng)波軌道速度和Bragg波相速度,將ULOS徑向速度做地距投影之后,再去除后兩項(xiàng)速度才能得到海流的距離向速度。

ATI-SAR獲取兩幅圖像的時(shí)間差應(yīng)該充分長(zhǎng),以便得到可靠的相位差,但同時(shí)也要比后向散射區(qū)域的去相干時(shí)間短[17]。海面后向散射區(qū)域的去相干時(shí)間τcorr可以近似表達(dá)為:

(4)

式中σorb為重力波軌道速度的標(biāo)準(zhǔn)差,可近似為σorb=0.068U10,其中U10是距離海面10 m高度處的風(fēng)速[18]。假設(shè)風(fēng)速大小為6～12 m/s,信號(hào)波長(zhǎng)為0.031 m,則海面后向散射區(qū)域的去相干時(shí)間為4.3～8.6 ms。

2 順軌有效基線計(jì)算方法

順軌有效基線的計(jì)算需要軌道方程的確定、圖像對(duì)偏移量的獲取、基線矢量的投影及有效基線方程的建立四個(gè)步驟。

2.1 軌道方程的確定

衛(wèi)星軌道如圖1所示意,由于軌道數(shù)據(jù)時(shí)間分辨率較高,采用多項(xiàng)式擬合,確定衛(wèi)星位置方程和衛(wèi)星速度方程[19]。

圖1 衛(wèi)星飛行軌道示意圖Fig.1 Schematic diagram of satellite flight orbit

在圖1中,A和B分別表示主輔衛(wèi)星對(duì)應(yīng)的軌道,從衛(wèi)星數(shù)據(jù)頭文件可知,每個(gè)衛(wèi)星提供了時(shí)間間隔為1 s 的軌道參數(shù),包括在WGS-84坐標(biāo)系下的空間位置和飛行速度。例如Pm為軌道A(主衛(wèi)星)的某一點(diǎn),其空間直角坐標(biāo)為(Xm,Ym,Zm),在三個(gè)坐標(biāo)軸上飛行速度為(Vxm,Vym,Vzm),Ps為軌道B(輔衛(wèi)星)的某一點(diǎn),其空間直角坐標(biāo)為(Xs,Ys,Zs),其飛行速度為(Vxs,Vys,Vzs),用一段時(shí)間內(nèi)的主輔軌道參數(shù)進(jìn)行擬合,采用三次多項(xiàng)式,以時(shí)間為自變量:

(5)

(6)

式中:(X,Y,Z)、(Vx,Vy,Vz)為衛(wèi)星的位置和飛行速度;ai,bi,ci,di,ei,fi為多項(xiàng)式系數(shù)(i=0,1,2,3)。

2.2 圖像間偏移量的獲取

在InSAR圖像中方位向第n行的成像時(shí)刻t(n)為:

(7)

式中:PRF為脈沖重復(fù)頻率;t0是成像起始時(shí)刻。由于主輔衛(wèi)星有各自的軌道方程(5)和方程(6),而原始的主輔圖像可能存在同一行號(hào)像元對(duì)應(yīng)的不是同一地面成像點(diǎn),因此,在使用輔衛(wèi)星軌道方程時(shí),使用的行號(hào)n需要和主衛(wèi)星的行號(hào)配準(zhǔn)。

解決主輔圖像行號(hào)的統(tǒng)一問題,首先需要進(jìn)行圖像對(duì)偏移量的確定。具體方法為,選取主衛(wèi)星圖像的場(chǎng)景中心為控制點(diǎn),利用相干系數(shù)作為匹配測(cè)度對(duì)輔衛(wèi)星圖像的偏移量進(jìn)行確定。相干系數(shù)γ的計(jì)算公式為:

(8)

式中:N和M為計(jì)算相干性的圖像塊尺寸大小;n′和m′為數(shù)據(jù)塊內(nèi)行列號(hào);μ1(n′,m′),μ2(n′,m′)為主輔影像數(shù)據(jù)塊內(nèi)影像坐標(biāo)(n′,m′)處的復(fù)數(shù)值;|·|2為數(shù)據(jù)的二階范數(shù)。

確定偏移量分為像元級(jí)和亞像元級(jí)。像元級(jí)偏移量是在主圖像中以控制點(diǎn)為中心選取一個(gè)匹配窗口,再在輔圖像中選取一個(gè)相對(duì)較大的搜索窗口,在搜索窗中遍歷匹配窗口,分別計(jì)算各自的相干系數(shù),以相干系數(shù)最大的遍歷位置確定圖像像元偏移量,即中心像元位置差,以此決定要調(diào)整輔圖像的行號(hào)數(shù)。

ATI-SAR成像間隔是毫秒量級(jí),在此期間海面隨機(jī)運(yùn)動(dòng)導(dǎo)致后天線記錄的海浪紋理相比于前天線已經(jīng)產(chǎn)生了一定的變化,圖像間海浪紋理的相關(guān)性會(huì)急劇下降[20],因此在利用相關(guān)系數(shù)法獲取圖像間偏移量時(shí),所選區(qū)的匹配窗口應(yīng)包含陸地區(qū)域,以減小偏移量計(jì)算的誤差。

為提高順軌有效基線計(jì)算精度,還可以進(jìn)一步計(jì)算亞像元級(jí)偏移量。首先對(duì)匹配窗口和搜索窗口進(jìn)行過采樣,采用sinc函數(shù)插值的方法,插值間隔取0.1個(gè)像元。用同樣的遍歷方法可以確定亞像元級(jí)偏移量,此時(shí)主輔圖像統(tǒng)一后的行號(hào)變?yōu)樾?shù)級(jí)。

2.3 基線矢量的投影

順軌有效基線幾何示意如圖2所示。

圖2 有效基線幾何示意圖Fig.2 Schematic diagram of effective baseline geometry

(9)

當(dāng)只有一個(gè)天線用于發(fā)射,兩個(gè)天線都用于接收時(shí),順軌有效基線Be是順軌物理基線的一半。

2.4 有效基線方程的建立

目前常用的方法是將成像區(qū)域中心行對(duì)應(yīng)的兩星位置及時(shí)刻和衛(wèi)星速度等信息獲得中心行對(duì)應(yīng)的有效基線,并將該基線值應(yīng)用于整個(gè)刈幅。然而,衛(wèi)星按照軌道運(yùn)行時(shí),由于各方面的原因?qū)е滦l(wèi)星姿態(tài)不穩(wěn)定,對(duì)方位向各行像素點(diǎn)成像時(shí)兩個(gè)衛(wèi)星之間并不是一個(gè)固定的基線值,使得獲取同一對(duì)干涉影像時(shí)順軌有效基線的長(zhǎng)度是一個(gè)隨方位向上行號(hào)的變量。

為了解決上述問題,需建立順軌有效基線與方位向上行號(hào)之間的關(guān)系式,以便獲得每一行的順軌有效基線值。

通過式(5)、(6)、(9)可推導(dǎo)出針對(duì)單發(fā)雙收模式順軌有效基線Be與成像行號(hào)的函數(shù)關(guān)系為:

(10)

式中tm和ts分別為主衛(wèi)星和輔衛(wèi)星對(duì)同一目標(biāo)點(diǎn)成像的時(shí)刻,即:

tm=[(n-1)/PRF]+tm0。

(11)

ts=[(n-1)/PRF]+ts0。

(12)

其中tm0為主圖像第一行成像時(shí)間,ts0為輔圖像與主圖像配準(zhǔn)后的第一行成像時(shí)間。

3 應(yīng)用結(jié)果分析

將以上有效基線的計(jì)算方法應(yīng)用于TanDEM-X/TerraSAR-X衛(wèi)星,時(shí)間為2012年3月19日6:41,地點(diǎn)為奧尼克群島地區(qū)。主圖像的大小為20 782×18 432,即行號(hào)從1到20 782。選取從2012年3月19日6:41:15—6:41:32時(shí)間段,采用三次多項(xiàng)式擬合主衛(wèi)星位置方程和速度方程分別為式(13)和式(14):

(13)

(14)

輔衛(wèi)星的位置方程為式(15):

(15)

使用相干系數(shù)法獲得方位向平均偏移量為0.141,把午夜零時(shí)設(shè)定為0 s,則主衛(wèi)星成像起始時(shí)刻為24 079.843 0,輔衛(wèi)星成像起始時(shí)刻為24 079.838 8。將行號(hào)(1—20782)逐個(gè)通過式(7)獲得所在行對(duì)應(yīng)的主輔衛(wèi)星成像的各自時(shí)刻,再通過軌道方程(13)、(14)、(15)獲得主輔衛(wèi)星的瞬時(shí)位置和主衛(wèi)星瞬時(shí)速度。通過公式(10)即可得到對(duì)應(yīng)各行的有效基線Be。結(jié)果如圖3(a)所示。

圖3 順軌有效基線與方位向行號(hào)關(guān)系圖Fig.3 The relationship between the effective baseline along the track and the azimuth line number

針對(duì)TanDEM-X/TerraSAR-X衛(wèi)星獲取的另一對(duì)順軌干涉數(shù)據(jù)(時(shí)間:2014年9月1日6:51,位置:蘇格蘭西南的艾萊島附近),主圖像的大小為:28 200×12 790。利用同樣方法獲得的順軌有效基線與主影像方位向行號(hào)關(guān)系如圖3(b)所示。

從圖3(a)可以看出,在整幅圖像上,順軌有效基線從41.7～43.6 m,變化范圍約1.9 m。從另一對(duì)數(shù)據(jù)的結(jié)果(見圖3(b))看,順軌有效基線的變化范圍也同樣接近2 m。而通過傳統(tǒng)方法計(jì)算的只是中心行對(duì)應(yīng)的順軌有效基線,兩對(duì)干涉數(shù)據(jù)的順軌有效基線分別為42.6和-17.7 m,與參考文獻(xiàn)[10]所述順軌有效基線大小接近。由于衛(wèi)星軌道高度距地都有幾百公里以上,無論是交軌干涉還是順軌干涉,干涉結(jié)果對(duì)基線長(zhǎng)度的變化十分敏感,將以上固定值應(yīng)用于整幅圖像無疑會(huì)帶來較大誤差。

從圖3(a)和3(b)還可看出,Be與圖像方位向行號(hào)呈線性變化,這不是必然結(jié)果,可能是由于兩衛(wèi)星在該幅圖像成像期間,其姿態(tài)穩(wěn)定,始終保持了平行。針對(duì)這樣的結(jié)果,在后期應(yīng)用時(shí),可以使用擬合的線性方程來替代公式(10)的繁瑣計(jì)算。兩對(duì)干涉數(shù)據(jù)的順軌有效基線與行號(hào)的簡(jiǎn)易線性方程如式(16)和式(17)所示:

Be=41.685 7+9.184 7×10-5n,

(16)

Be=-18.449 8+5.776 9×10-5n。

(17)

通過TanDEM-X/TerraSAR-X衛(wèi)星實(shí)測(cè)數(shù)據(jù)(時(shí)間:2012年3月19日6:41時(shí),位置:奧尼克群島地區(qū))進(jìn)行海面流速估計(jì),以判斷逐行順軌有效基線對(duì)海面流場(chǎng)精度的影響。逐行使用不同的順軌有效基線,海流流速反演結(jié)果如圖4所示。

圖4 奧尼克群島地區(qū)海流流速反演圖Fig.4 Inversion map of ocean current velocity in the Onik Islands area

分別使用傳統(tǒng)基線計(jì)算方法和改進(jìn)的基線計(jì)算方法進(jìn)行海流流速反演,計(jì)算兩幅圖像方位向上每一行的流速均值,其結(jié)果如圖5(a)所示,在方位向上每一行的海流流速差值如圖5(b)所示。由圖5(b)可知,在方位向上逐行使用不同的順軌有效基線值能夠一定程度提高海表面流速反演精度。

圖5 兩種算法流速反演結(jié)果對(duì)比Fig.5 Comparison of flow velocity inversion results between the two algorithms

4 結(jié)語

本文通過軌道方程的確定、圖像對(duì)偏移量的獲取、基線矢量的投影及有效基線方程的建立四個(gè)步驟完成順軌有效基線與方位向行號(hào)的函數(shù)關(guān)系式的建立,用以在主輔圖像干涉時(shí)順軌有效基線的精細(xì)化應(yīng)用是對(duì)傳統(tǒng)計(jì)算方法的改進(jìn),有效避免在整幅圖像上使用固定基線帶來的誤差。通過TanDEM-X/TerraSAR-X實(shí)測(cè)干涉數(shù)據(jù)對(duì)所述改進(jìn)方法進(jìn)行了應(yīng)用和檢驗(yàn),同時(shí)針對(duì)所使用干涉圖像對(duì)的特殊線性關(guān)系建立了快捷方便的線性方程,避免了繁瑣的計(jì)算,更有利于提高程序化運(yùn)行效率。通過在整幅圖像上精細(xì)化使用順軌有效基線會(huì)對(duì)諸如海流反演方面有效提高反演精度。