馬 飛,隋立春,姚頑強
(1.長治學(xué)院 計算機系,山西 長治 046011; 2.長安大學(xué) 地質(zhì)與測繪工程學(xué)院,西安 710054;3.西安科技大學(xué) 測繪科學(xué)與技術(shù)學(xué)院,西安 710054)
在合成孔徑雷達(dá)(SAR)衛(wèi)星技術(shù)的應(yīng)用中,衛(wèi)星軌道參數(shù)是合成孔徑雷達(dá)干涉測量(InSAR, interferometric synthetic aperture radar)數(shù)據(jù)處理過程中重要參數(shù)之一,可直接影響主輔SAR圖像配準(zhǔn)、DEM測高、基線估算等精度。在干涉測量數(shù)據(jù)處理過程中,最初主影像和輔影像粗配準(zhǔn)過程中會引入衛(wèi)星軌道參數(shù)作為影像初始偏移量,去除平地相位過程中會利用衛(wèi)星軌道參數(shù)模擬計算平地相位,生成DEM過程中會利用衛(wèi)星軌道參數(shù)估算基線長度,因此衛(wèi)星軌道參數(shù)是干涉合成孔徑雷達(dá)的基礎(chǔ),它直接決定著合成孔徑雷達(dá)干涉測量系統(tǒng)的特性優(yōu)劣[1]。
大部分SAR衛(wèi)星運營商都能提供衛(wèi)星軌道精密參數(shù),例如歐空局發(fā)射的ERS衛(wèi)星、ENVISAT衛(wèi)星由荷蘭代夫特對地空間研究中心(DEOS, delft institute earth—oriented space research)提供的精密軌道數(shù)據(jù)精度優(yōu)于6 cm[2];哨兵衛(wèi)星POD精密定軌星歷數(shù)據(jù)定位精度優(yōu)于5厘米[3];TerraSAR X衛(wèi)星搭載全球定位系統(tǒng)(GPS)接收機和激光精密定軌反射器,軌道確定精度達(dá)到厘米級[4]。因此關(guān)于SAR衛(wèi)星軌道插值的文獻(xiàn)研究較少,很多學(xué)者專注于研究InSAR數(shù)學(xué)模型、軌道精度對InSAR測量結(jié)果的影響等,分析InSAR干涉測量的數(shù)學(xué)模型,包括軌道精度、距離向頻譜、干涉臨界基線距、模糊高度、差分相位對形變的敏感度等問題,但未就如何改進(jìn)軌道精度提出解決方案。針對衛(wèi)星軌道插值問題更多的文獻(xiàn)專注于研究全球定位系統(tǒng)(GPS, global positioning system)、 北斗導(dǎo)航衛(wèi)星定位系統(tǒng)(BDS, beiDou navigation satellite system)和伽利略(Galileo)衛(wèi)星定位系統(tǒng)的軌道插值問題,插值方法包括埃爾米特插值、切比雪夫插值、三次樣條插值、拉格朗日插值算法等。
針對未提供精密軌道文件的SAR衛(wèi)星數(shù)據(jù),嚴(yán)重制約著該衛(wèi)星影像數(shù)據(jù)處理精度。例如加拿大的RADARSAT衛(wèi)星,日本的JERS衛(wèi)星、ALOS衛(wèi)星和ALOS-2衛(wèi)星等等,其中ALOS衛(wèi)星搭載的PALSAR傳感器掃描地面影像的持續(xù)時間大概為13秒,在影像的參數(shù)文件中提供了衛(wèi)星從生成第一個像素到最后一個像素期間的所有衛(wèi)星狀態(tài)矢量,包括衛(wèi)星三維位置、三維移動速率等數(shù)據(jù),其狀態(tài)矢量采樣間隔為60秒;因此,主輔SAR影像在進(jìn)行干涉測量時只有一個衛(wèi)星狀態(tài)矢量參與了計算。這樣在利用軌道參數(shù)進(jìn)行去除橢球參考相位、像素點匹配、去除平地相位時都會遇到困難,干涉相位圖中會有殘余條紋,造成最終獲取的DEM或差分干涉圖中混有系統(tǒng)誤差。因此本文擬采用編程工具對衛(wèi)星軌道數(shù)據(jù)進(jìn)行模擬,利用埃爾米特插值法對其軌道參數(shù)進(jìn)行等距插值,從而估計其軌道參數(shù),縮小衛(wèi)星軌道狀態(tài)矢量采樣間隔,提高軌道精度。
插值算法是一種關(guān)于函數(shù)逼近的簡單又實用的數(shù)學(xué)方法,主要用于某函數(shù)在有限個點處取值狀況,估計其在其他點處的值。插值算法是數(shù)值分析中一種基本算法。根據(jù)被插值的函數(shù)自變量個數(shù),插值法可分為一維插值、二維插值、和多維插值。根據(jù)選用的分段直線、多項式或樣條函數(shù)插值函數(shù),插值又可分線性、多項式和樣條插值[5]。
目前較為常用的插值算法有Lagrange插值、牛頓插值、埃爾米特插值、三次樣條插值、高維插值等。Lagrange插值、牛頓插值主要對函數(shù)進(jìn)行直接插值,對函數(shù)的一階或多階導(dǎo)數(shù)沒有限制,沒有考慮插值點的導(dǎo)數(shù)值,即沒有考慮衛(wèi)星軌道的瞬時速度矢量,不適用于SAR衛(wèi)星軌道插值;埃爾米特插值及三次樣條插值不僅對原函數(shù)進(jìn)行限制,根據(jù)限制條件的不同,可以對插值函數(shù)的倒數(shù)進(jìn)行限制,如埃爾米特插值可以對函數(shù)的一階導(dǎo)數(shù)進(jìn)行限制,而三次樣條插值法可以根據(jù)不同的條件,可以對一階或二階導(dǎo)數(shù)進(jìn)行限制,從而可以保證插值函數(shù)的連續(xù)性,這兩種插值方法既可以擬合SAR衛(wèi)星軌道的位置矢量,也可以計算衛(wèi)星在某一位置的瞬時速度矢量,因此適用于SAR衛(wèi)星軌道插值。根據(jù)本次對衛(wèi)星軌道參數(shù)插值的條件,即在已知點位和速度矢量在X、Y、Z方向的分量條件下,在時間上進(jìn)行插值,故選擇的插值方法有埃爾米特插值和三次樣條插值法。武漢大學(xué)吳宏安教授[6]等人對InSAR數(shù)學(xué)模型研究發(fā)現(xiàn),利用原始軌道參數(shù)文件進(jìn)行埃爾米特插值,可以有效地加密軌道文件參數(shù),利用新的加密后的軌道參數(shù)文件進(jìn)行像素點配準(zhǔn)、基線估算、橢球參考相位估算等,可以提高干涉測量精度,去除平地殘余相位。本文利用埃爾米特插值法估算SAR衛(wèi)星軌道參數(shù),并進(jìn)行相關(guān)的實驗研究。
衛(wèi)星描述軌道模型的星歷數(shù)據(jù)包括衛(wèi)星的位置矢量(X,Y,Z)和瞬時的速度矢量(Vx,Vy,Vz),用下式來表示:
(1)
(2)
其中:(X,Y,Z)分別為固定笛卡爾參考坐標(biāo)系的三個坐標(biāo)軸,固定笛卡爾參考坐標(biāo)系統(tǒng)是以地心為原點,X、Y軸垂直位于赤道平面,X軸由地心指向格林威治子午線,Z軸垂直于赤道平面且符合右手定律,與地球橢球極軸重合指向北極。
衛(wèi)星軌道參數(shù)文件中的狀態(tài)矢量可以進(jìn)行多項式擬合,將(X,Y,Z)和(Vx,Vy,Vz)作為時間T的函數(shù),構(gòu)建矩陣如下式:
(3)
式中,a,b,c分別為X,Y,Z維的多項式系數(shù)。
當(dāng)求一個插值多項式H(x),使其滿足條件為:
H(xi)=yi,i=1,2,…,n
(4)
(5)
根據(jù)函數(shù)的插值條件,共有m+n+2個。所構(gòu)造函數(shù)H(x)次數(shù)一般小于或等于m+n+1次,可由牛頓插值的思想來構(gòu)造H(x)為:
H(x)=Nn(x)+Pm(x){(x-x0)(x-x1)…(x-xn)}
(6)
其中:Nn(X)為牛頓基本的差值多項式,Pm(X)為函數(shù)對應(yīng)的m次多項式。根據(jù)式(6)和已知數(shù)據(jù)可以得到:
H(xi)=Nn(xi)=yi,i=1,2,…,n
(7)
在確定Pm(X),需先對H(x)求導(dǎo):
(8)
其中:ωn+1(x)=(x-x0)(x-x1)…(x-xn),令x=xi,i=0,1,2,3,…,m,且把H′(xi)=y′,i=1,2,3,…,m代入式上式可得:
(9)
進(jìn)一步可得Pm(X):
(10)
對于求Pm(X)轉(zhuǎn)化成為已知Pm(x)的函數(shù)表,如表1所示。
表1 Pm(x)的數(shù)據(jù)表
當(dāng)Pm(X)不超過m的插值多項式。Pm(x)滿足式:
Pm(x)=Pm(x0)+Pm[x0,x1](x-x1)+…+
Pm[x0,…,xm](x-x0)…(x-xm-1)=
(11)
其中:Pm[x0,…,xm]表示為第m階差商,形式為:
(12)
再將Nn(X)、Pm(x)帶入式(6),可得到插值函數(shù)H(x)。
本論文采用埃爾米特等距插值法,即要求在給定節(jié)點處函數(shù)值和一階導(dǎo)數(shù)均相等[7]。滿足下式:
yi=f(xi)
(13)
(14)
以上為埃爾米特插值算法基本原理,本文基于軟件編程實現(xiàn)對SAR衛(wèi)星軌道參數(shù)插值計算。
InSAR技術(shù)測量是通過對地表目標(biāo)點的兩個視向的觀測實現(xiàn)的,其工作方式包括距離向干涉測量(Across-track Interferometry)、方位向干涉測量(Along-track Interferometry)和重軌干涉測量(repeat-pass Interferometry)。
距離向干涉測量通過在飛行平臺裝載兩個天線,一個天線用于發(fā)射并接收雷達(dá)波,另外一個僅用于接收雷達(dá)波。通過這種干涉測量方式,只要確定平臺位置,就較容易地獲得高質(zhì)量的干涉測量數(shù)據(jù)及較高精度的DEM。當(dāng)前這種干涉測量方式常用于航空平臺上,即為前面所提到了機載InSAR。在采用這種測量方式,要求飛行平臺較為穩(wěn)定,而且需要飛機姿態(tài)的影響。
方位向干涉測量,該方式將兩個天線安裝在平臺的前端和后端,主要適用于機載,已應(yīng)用在運動目標(biāo)的探測、海上波浪譜線等方面的測量工作。
本次用于研究的影像為基于衛(wèi)星的重軌干涉測量獲取的影像。重軌干涉測量,其作業(yè)方式僅需要一個天線,利用衛(wèi)星在短時間回到大致相同的軌道上,獲取同一范圍內(nèi)的主輔影像數(shù)據(jù)。該測量適合于幾百公里上的衛(wèi)星作業(yè)方式。一般而言,星載SAR要獲得較好的影像數(shù)據(jù),需滿足以下3個條件:首先是在兩次觀測期間,地面沒有發(fā)生形變或僅發(fā)生微小的形變;其次平臺與地表目標(biāo)間具有較穩(wěn)定的觀測幾何關(guān)系,這需要搭載平臺的穩(wěn)定和平臺的準(zhǔn)確定位參數(shù),在獲取相同的影像,當(dāng)衛(wèi)星有精密的軌道參數(shù)時,處理后的影像精度更高;最后為SAR處理器能夠?qū)ψ鲞\動的補償信號保留住內(nèi)在的相位信息。
數(shù)據(jù)預(yù)處理,該階段主要從主輔影像軌道參數(shù)文件等生成軌道參數(shù)文件和單視復(fù)影像。當(dāng)研究區(qū)域為全影像的一部分時,可以對原始影像進(jìn)行裁剪,將得到裁剪區(qū)域用于后面的處理研究。
數(shù)據(jù)配準(zhǔn),該過程可分為去噪、控制點自動搜索、重采樣與配準(zhǔn)質(zhì)量評價等,主要目的是將輔影像與主影像進(jìn)行匹配,使得兩景影像坐標(biāo)系統(tǒng)一致。在處理影像過程中,需采用多級配準(zhǔn)方案來提高影像的配準(zhǔn)精度。一般分為粗配準(zhǔn)、精配準(zhǔn)及子像元配準(zhǔn)。
干涉圖生成,需進(jìn)行相位干涉,即對兩景影像的匹配點的數(shù)值做共軛相乘,得到地面任意點的干涉后影像圖的數(shù)值。通常計算得到的相位差由干涉圖的灰度表示,并以條紋形狀表現(xiàn)。
重采樣與配準(zhǔn)質(zhì)量評價。重采樣指對對原始影像數(shù)據(jù)的實部和虛部做內(nèi)插,插值方法有雙線性和雙三次樣條方法。為了獲得高精度的相位值,對重采樣提出較高的要求。配準(zhǔn)結(jié)果的評價是對配準(zhǔn)結(jié)果的綜合評價,對干涉圖的生成有重要影響。質(zhì)量評價的量化標(biāo)準(zhǔn)為相干系數(shù),其取值范圍[-1,1],為配置質(zhì)量評價量化標(biāo)準(zhǔn)。當(dāng)相干系數(shù)越接近0,代表著兩景數(shù)據(jù)相關(guān)程度越小;當(dāng)相干系數(shù)的絕對值越接近1時,表示兩景數(shù)據(jù)的相關(guān)程度越大。
去平地效應(yīng),當(dāng)?shù)孛婺硟牲c高程不同,其相位差可能有兩部分造成的,其一為平地兩點的相位差;其二為兩點間存在的高差而造成相位差。通過去平地效應(yīng),使得到干涉圖的相位差只由地形引起,從而得到更為準(zhǔn)確的地形數(shù)據(jù)。
相位解纏,即從干涉圖得到的相位差,取值范圍為(-π,π],代表了各個相位差的主值,為了得到真實相位差需在這個主值的基礎(chǔ)上減去2π或加上2π的整數(shù)倍。相位解纏的算法可分為兩大類,一類為基于路徑積分相位解纏方法(主要為枝切法);另一類為基于最小范數(shù)的解纏算法。
地理編碼,由干涉圖得到高程值后,只是得到各點在SAR坐標(biāo)系下的高程值。為了得到對應(yīng)坐標(biāo)系下的數(shù)字高程模型,需將SAR坐標(biāo)系下的每一點高程值進(jìn)行坐標(biāo)轉(zhuǎn)換,得到最終的地面DEM。
巴姆地區(qū)位于伊朗的東南部,中心位置為北緯29.01°,東經(jīng)58.26°。本文利用覆蓋該地區(qū)的Envisat衛(wèi)星ASAR傳感器影像和DEOS提供的精密軌道文件對埃爾米特軌道插值實驗進(jìn)行分析驗證。根據(jù)雷達(dá)干涉測量原理[15],對兩幅Envisat衛(wèi)星主輔影像進(jìn)行干涉處理,先進(jìn)行主輔影像配準(zhǔn),計算主輔影像上所有的像素點在參考托球面上的投影點,對所有配準(zhǔn)的像素點進(jìn)行相位干涉,獲取各像素的相位值,該相位值中包含地形起伏引起的相位、地表形變引起的相位變化、橢球參考面相位、噪聲相位等等。本次數(shù)據(jù)處理過程中像素點配準(zhǔn)和橢球參考面相位都受到軌道文件精度的影響,因此干涉相位圖的參考面相位的條紋規(guī)律性和相干性值是軌道參數(shù)文件質(zhì)量定性和定量評判的重要標(biāo)準(zhǔn)。
如圖1所示,圖(a)為基于原始軌道文件進(jìn)行干涉測量得到的干涉圖,圖(b)為基于埃爾米特插值軌道文件進(jìn)行干涉測量得到的干涉圖,圖(c)為基于DEOS提供的精密軌道文件干涉測量得到的干涉圖。
圖1 基于不同軌道參數(shù)得到的干涉圖
對比圖(a),(b)兩幅干涉圖可知,兩幅圖中均出現(xiàn)了清晰且規(guī)律性很強的干涉條紋,其中心位置出現(xiàn)的蝴蝶狀條紋是地震形變引起的相位變化,屬于形變條紋;其他規(guī)律性條紋主要是參考面相位條紋;在圖中左下角的藍(lán)線區(qū)域的干涉條紋對比具有較大的差距,該區(qū)域地形起伏較大,屬于地形相位條紋和參考面相位條紋疊加的結(jié)果。圖(b)中經(jīng)軌道插值后得到的干涉圖干涉條紋更清晰,在地形復(fù)雜的地區(qū)仍具有較高的相干性,在形變區(qū)的干涉條紋具有較強的連貫性。分析圖(b),(c)可知,由插值得到的干涉圖和利用衛(wèi)星精密軌道得到的干涉圖有較強的一致性,尤其是在左下角地形起伏較大的區(qū)域,干涉條紋基本無差異,說明利用DEOS精密軌道文件和埃爾米特插值軌道文件進(jìn)行干涉測量處理后相位無差異,埃爾米特插值軌道文件的精度與DEOS精密軌道文件精度相當(dāng)。
通過對以上結(jié)果分析,可以定性地判斷基于埃爾米特衛(wèi)星軌道參數(shù)的插值算法,可以提高SAR影像的干涉測量效果,有效地提高軌道文件精度。
由于歐空局會提供Envisat衛(wèi)星的DEOS精密軌道文件,所以使用埃爾米特插值對原始文件插值計算精密軌道文件的研究價值并不在于為Envisat衛(wèi)星軌道插值,而在于為無精密軌道文件的SAR衛(wèi)星插值軌道數(shù)據(jù),例如日本ALOS衛(wèi)星。因此本文再次使用埃爾米特插值算法對ALOS衛(wèi)星軌道進(jìn)行插值實驗。
本文選取的ALOS/PALSAR數(shù)據(jù)為覆蓋陜西咸陽彬縣和長武縣的彬長礦區(qū)的影像。該礦儲煤豐富,質(zhì)量優(yōu)良,礦區(qū)的西北、北部以甘肅省為界,東西距離為46 km,南北長約為36 km,礦區(qū)地表覆蓋面積近1 000 km2。該礦區(qū)位于黃土高原的塬梁溝壑地帶,地表覆有黃土層。氣候干燥,降水量少,年平均氣溫在10°,高差多達(dá)200 m。該地區(qū)的屬于涇河流域,其中流經(jīng)該礦區(qū)的最大支流為黑河。彬長礦區(qū)的地表植被相對較少,位于開采沉陷區(qū)域內(nèi)的村落較少,有利于獲取高質(zhì)量的SAR影像。ALOS/PALSAR數(shù)據(jù)的原始參數(shù)文件中包含有11個軌道狀態(tài)矢量,但由于其軌道矢量時間間隔為60秒,每景影像像素掃描時間為13秒,因此干涉測量過程中有可能無軌道狀態(tài)矢量參與。本文采用埃爾米特插值法進(jìn)行等距插值,設(shè)置時間間隔為10秒、5秒、2秒分別計算插值后的軌道狀態(tài)矢量,在插值過程中以時間T為自變量,衛(wèi)星的位置矢量作為函數(shù)值,速度矢量作為時間的一階導(dǎo)數(shù)值。
具體實驗數(shù)據(jù)影像參數(shù)如表2所示,利用編程工具對實驗數(shù)據(jù)計算加密后的軌道參數(shù),再將新的軌道參數(shù)文件替換原始文件進(jìn)行干涉測量數(shù)據(jù)處理。
表2 軌道插值運算實驗數(shù)據(jù)
在ALOS/PALSAR數(shù)據(jù)的原始參數(shù)文件中,包含有影像起始時間(start_time)、結(jié)束時間(end_time)、中點位置時間(center_time),11個采樣時間間隔為60秒的軌道矢量(X,Y,Z,Vx,Vy,Vz,T),本文以這11個軌道狀態(tài)矢量為已知點構(gòu)建埃爾米特基函數(shù),以第一個軌道狀態(tài)矢量采樣時間為起始時間,以最后一個軌道狀態(tài)矢量采樣時間為結(jié)束,分別進(jìn)行10秒、5秒、2秒等間距插值,利用軌道矢量插值后干涉測量結(jié)果如圖2所示。
圖2 不同軌道參數(shù)采樣間隔的地理編碼影像干涉相位測量圖
通過觀察圖2中的(a),(b),(c),(d)四幅圖發(fā)現(xiàn)圖中的干涉條紋基本無差異,通過分析InSAR數(shù)學(xué)模型[7-8]可知,日本ALOS衛(wèi)星的PALSAR傳感器使用L波段對地掃描,L波段波長相對C波段較長,干涉測量后的干涉條紋對地形和形變的敏感度變低,地形高程變化178米才會引起一個周期的相位條紋變化。因此在(a),(b),(c),(d)四幅圖中的干涉條紋變化周期較圖1中的ENVISAT衛(wèi)星的C波段變化周期較少,四幅圖中的干涉條紋肉眼基本觀察不到差異。需要引入定量分析方法分析不同軌道采樣間隔下干涉測量的效果又和不同。
圖3 不同衛(wèi)星軌道參數(shù)采樣間隔的影像的相干圖
為了進(jìn)一步比較四種軌道參數(shù)采樣間隔對影像相干系數(shù)的影響,埃爾米特插值法得到的衛(wèi)星軌道參數(shù)采樣間隔為60秒、10秒、5秒、2秒的相干圖,選用編程工具軟件分別對該相干圖做灰度均值運算,求得每幅相干圖的平均系數(shù),如表3所示。
表3 不同采樣間隔下影像相干圖的平均相干系數(shù)表
通過表3可知當(dāng)衛(wèi)星軌道參數(shù)采樣間隔縮小時,兩幅SAR影像干涉得到平均相關(guān)系數(shù)比原始采樣間隔60秒的影像相干圖的平均相干系數(shù)均高,其代表對原始衛(wèi)星軌道參數(shù)進(jìn)行插值的方法對提高SAR影像的相干性具有可行性,并取得一定成效。再進(jìn)一步分析,隨著軌道參數(shù)采樣間隔的不斷減少,提高影像的相干性并非呈正相關(guān),本次實驗中,當(dāng)采樣間隔為5秒時,得到相干圖的平均相干系數(shù)最高,而其采樣間隔為2秒時,其相干圖的平均相關(guān)系數(shù)略低。這說明隨著軌道參數(shù)間隔的減小,其相干性呈先提高后降低的趨勢,以此估計得到礦區(qū)地表點也呈該沉降趨勢。本次研究中,原始軌道參數(shù)間隔為60秒,經(jīng)插值算法采樣間隔為10秒、5秒、2秒,初步判斷當(dāng)采樣間隔為5秒時,得到相干性最好的差分干涉影像圖。
1)本文基于埃爾米特插值多項式擬合SAR衛(wèi)星軌道,對SAR衛(wèi)星軌道模型的中的位置矢量和瞬時的速度矢量進(jìn)行插值計算,可以有效地加密衛(wèi)星軌道星歷數(shù)據(jù),提高衛(wèi)星軌道參數(shù)文件精度,減少干涉測量過程中的系統(tǒng)誤差。
2)通過對Envisat衛(wèi)星原始軌道文件進(jìn)行埃爾米特插值計算,得到加密后的衛(wèi)星軌道星歷數(shù)據(jù)。利用原始軌道文件、加密軌道文件、DEOS精密軌道文件分別進(jìn)行干涉測量,定性研究發(fā)現(xiàn),利用加密軌道文件和精軌文件得到的干涉圖中干涉條紋清晰且一致,說明埃爾米特插值算法可以有效加密SAR衛(wèi)星軌道星歷數(shù)據(jù),提高干涉測量精度。
3)通過對ALOS衛(wèi)星原始粗軌進(jìn)行埃爾米特插值計算,軌道星歷數(shù)據(jù)采樣間隔分別設(shè)置為10秒、5秒、2秒,結(jié)果表明進(jìn)行軌道插值后得到的干涉圖相干性比基于粗軌得到的干涉圖相干性顯著提升;隨著軌道參數(shù)間隔的縮小,其相干性呈先提高后降低的趨勢,采樣間隔為5秒時相干性最高。