王志偉，劉國(guó)林，陶秋香，于勝文，王 珂，周傳義

(山東科技大學(xué)測(cè)繪科學(xué)與工程學(xué)院，山東 青島 266590)

合成孔徑雷達(dá)干涉測(cè)量(interferometric synthetic aperture radar，InSAR)技術(shù)以其高精度、高空間分辨率的地表形變監(jiān)測(cè)能力，目前已經(jīng)被廣泛應(yīng)用于監(jiān)測(cè)各種因素導(dǎo)致的地表形變[1-7]。隨著越來(lái)越多的衛(wèi)星發(fā)射升空，同一地區(qū)可以獲取多種不同傳感器SAR影像，進(jìn)而可以獲取不同雷達(dá)視線向地表形變監(jiān)測(cè)結(jié)果(多視線向D-InSAR形變量)。將多視線向D-InSAR地表形變量作為觀測(cè)量，以衛(wèi)星軌道方位角和入射角等參數(shù)構(gòu)建系數(shù)矩陣，建立三維地表形變解算模型，可以獲取真實(shí)三維地表形變[12]。然而，由于受到不同誤差因素的影響，多視線向D-InSAR地表形變監(jiān)測(cè)結(jié)果的精度不盡相同。因此，基于多視線向D-InSAR形變量聯(lián)合平差，以獲取高精度三維形變解算結(jié)果尤為重要。其研究關(guān)鍵點(diǎn)是如何合理確定各觀測(cè)量的權(quán)重。

本文提出一種基于GPS約束的多視線向D-InSAR三維地表形變解算中定權(quán)方法。該方法以高精度GPS形變監(jiān)測(cè)結(jié)果作為約束，通過對(duì)多視線向D-InSAR觀測(cè)值進(jìn)行評(píng)價(jià)確定其權(quán)重。本文以2009年意大利L’Aquila地震三維地表形變解算為例，采用GPS約束的定權(quán)方法實(shí)現(xiàn)三維形變解算。試驗(yàn)結(jié)果表明，本文提出的GPS約束定權(quán)的多視線向D-InSAR三維形變解算方法可獲取更可靠的解算結(jié)果，可以解決不同精度觀測(cè)量聯(lián)合平差權(quán)重確定存在的定權(quán)難、定權(quán)適用性差、精度低等突出問題。

1 多視線向D-InSAR三維地表形變GPS約束定權(quán)方法

1.1 多視線向D-InSAR三維地表形變解算

SAR衛(wèi)星的成像幾何關(guān)系如圖1所示。AZ向表示衛(wèi)星飛行方向在地面上的投影，GR向表示衛(wèi)星視線向在地面的投影，α表示衛(wèi)星方位角，θ表示SAR影像入射角，U、E和N分別代表大地測(cè)量坐標(biāo)系中的垂直、東西和南北方向，dU，dE和dN分別代表對(duì)應(yīng)方向上的形變信息。

圖1 SAR衛(wèi)星三維成像示意圖

R=Bd

(1)

式中，B為三維形變模型系數(shù)矩陣，具體可以表示為

(2)

根據(jù)最小二乘原理，可以解算出三維形變信息

d=(BTPB)-1BTPR

(3)

1.2 基于GPS約束的定權(quán)方法

針對(duì)多視線向D-InSAR三維地表形變解算中存在的定權(quán)問題，本文提出以GPS精密形變觀測(cè)量為基準(zhǔn)數(shù)據(jù)，對(duì)多視線向D-InSAR觀測(cè)量進(jìn)行精度評(píng)價(jià)，綜合估算出軌道和地形誤差及噪聲等因素引起的誤差。由于GPS測(cè)站數(shù)量有限，無(wú)法逐像元進(jìn)行精度評(píng)價(jià)。顧及軌道和地形及噪聲等因素引起的誤差對(duì)形變監(jiān)測(cè)結(jié)果影響的特點(diǎn)，建立聯(lián)合擬合模型并計(jì)算擬合參數(shù)。根據(jù)確定出的擬合模型估算多視線向D-InSAR每個(gè)像元觀測(cè)量的方差，以計(jì)算出的方差的倒數(shù)確定權(quán)矩陣，建立加權(quán)最小二乘三維地表形變解算模型。具體步驟如下：

(1) 將GPS三維形變觀測(cè)量換算到D-InSAR視線方向，如下

(4)

(2) 根據(jù)雙觀測(cè)值之差標(biāo)準(zhǔn)差計(jì)算公式，計(jì)算出GPS測(cè)站對(duì)應(yīng)D-InSAR觀測(cè)量方差

(5)

(3) 根據(jù)步驟(2)計(jì)算出的雙觀測(cè)值之差標(biāo)準(zhǔn)差及對(duì)應(yīng)像元相干值，按照下式計(jì)算方差

(6)

式中，γi為第i個(gè)像元相干系數(shù)；E(γGPS)為GPS觀測(cè)站對(duì)應(yīng)的SAR影像像元相干系數(shù)均值。

(4) 以步驟(3)估算出的D-InSAR像元的方差為依據(jù)，充分考慮軌道和地形誤差對(duì)形變量影響的特點(diǎn)，建立方差與像元點(diǎn)的位置和高程聯(lián)合擬合模型

(7)

式中，等號(hào)右側(cè)前6項(xiàng)為軌道誤差項(xiàng)；x表示第i個(gè)像元方位向坐標(biāo)；y表示第i個(gè)像元距離向坐標(biāo)；后兩項(xiàng)表示地形誤差項(xiàng)，z表示坐標(biāo)為(x,y)像元的高程值，ai(i=1,2,…,5)和bj(j=1,2)分別表示待估參數(shù)。

估算出上式中系數(shù)參數(shù)后便可以確定出擬合模型，按照上述模型估算每個(gè)像元D-InSAR形變監(jiān)測(cè)結(jié)果方差。

(5) 根據(jù)D-InSAR形變監(jiān)測(cè)結(jié)果方差，按照下式定權(quán)

(8)

按照式(8)計(jì)算出權(quán)陣P之后，按照式(3)進(jìn)行三維地表形變解算。

2 L’Aquila地震多視線向D-InSAR地表形變監(jiān)測(cè)

2.1 試驗(yàn)區(qū)與數(shù)據(jù)源

2009年4月6日意大利中部城市L’Aquila發(fā)生了Mw6.3級(jí)地震，整個(gè)意大利中部地區(qū)均有震感。本文以該地震震中區(qū)域?yàn)檠芯繀^(qū)驗(yàn)證定權(quán)方法的有效性?？紤]干涉對(duì)的選取原則(足夠短的垂直基線和時(shí)間基線)，選取了4個(gè)干涉像對(duì)SAR數(shù)據(jù)，分別用P1、P2、P3和P4表示，影像覆蓋范圍如圖2所示。同期的GPS數(shù)據(jù)來(lái)源于文獻(xiàn)[13]。研究區(qū)范圍內(nèi)共有48個(gè)點(diǎn)，如圖2所示，其中黑點(diǎn)表示用來(lái)估算聯(lián)合擬合模型參數(shù)的GPS數(shù)據(jù)，白點(diǎn)用來(lái)驗(yàn)證解算精度。

表1 選取的SAR數(shù)據(jù)基本信息

圖2 研究區(qū)概況

2.2 數(shù)據(jù)處理

以震前、后獲取的SAR影像分別作為主、輔影像，采用雙軌D-InSAR技術(shù)對(duì)數(shù)據(jù)進(jìn)行處理，對(duì)主、輔影像進(jìn)行配準(zhǔn)、重采樣和共軛相乘處理，獲取干涉相位圖。采用30 m分辨率的SRTM DEM數(shù)據(jù)去除地形相位的影響。采用Goldstein濾波方法對(duì)干涉圖進(jìn)行噪聲濾除，獲取最終的干涉圖。采用最小費(fèi)用流的方法進(jìn)行相位解纏。解纏之后的干涉圖只包含形變相位、大氣相位和殘余噪聲相位。其中，大氣相位對(duì)干涉圖的影響主要體現(xiàn)在對(duì)流層延遲方面，往往跟地形高程有一定的相關(guān)性，通過GAMMA軟件的相關(guān)算法可以一定程度上去除大氣誤差相位。

圖3為地理編碼后的形變圖。從圖中可以看出，多視線向D-InSAR形變監(jiān)測(cè)結(jié)果在形變趨勢(shì)和量級(jí)上基本一致，震中和斷裂帶周邊形變?yōu)樨?fù)值表示該區(qū)域地面下沉(呈現(xiàn)深色)。多視線向D-InSAR形變監(jiān)測(cè)結(jié)果在雷達(dá)視線向上下沉的一致性說(shuō)明拉奎拉地震是以地表下沉為主。根據(jù)地質(zhì)調(diào)查的結(jié)果顯示，拉奎拉地區(qū)所在的阿爾卑斯地區(qū)斷層走向?yàn)槲鞅?東南方向，大小在40～50之間的西南傾向，正斷層。

圖3 不同軌道SAR影像獲取的形變圖

3 L’Aquila地震三維形變解算

3.1 三維地表形變解算

獲取了多視線向D-InSAR地表形變信息，按照第一部分介紹的數(shù)據(jù)處理方法進(jìn)行三維地表形變解算。值得注意的是兩種數(shù)據(jù)形變監(jiān)測(cè)結(jié)果具有不同的空間分辨率，本文將ASAR數(shù)據(jù)監(jiān)測(cè)結(jié)果重采樣至ALOS數(shù)據(jù)分辨率。入射角和方位角信息見表1。用來(lái)計(jì)算擬合參數(shù)的GPS數(shù)據(jù)具體分布如圖2中白色點(diǎn)所示。解算出的三維地表形變結(jié)果如圖4所示。

圖4 三維地表形變解算結(jié)果

從解算出的三維地表形變結(jié)果中可以看出，在震中主形變區(qū)內(nèi)，地表發(fā)生了東南方向的水平錯(cuò)動(dòng)。在斷層周邊的三維形變場(chǎng)中表現(xiàn)出了比較明顯的形變特征。在東西方向形變場(chǎng)的下半部分，地表向西發(fā)生了移動(dòng)，上側(cè)發(fā)生了向東的移動(dòng)，表現(xiàn)出比較明顯的形變錯(cuò)動(dòng)。在南北方向下半部分，地表向北發(fā)生了移動(dòng)，上側(cè)發(fā)生了向南的移動(dòng)，同樣表現(xiàn)出比較明顯的形變錯(cuò)動(dòng)。在垂直方向的右上方，地表表現(xiàn)為上升的運(yùn)動(dòng)，而在下方表現(xiàn)為沉降，垂直方向上也發(fā)生了明顯了錯(cuò)動(dòng)。

3.2 精度驗(yàn)證

為了說(shuō)明本文提出定權(quán)方法的有效性，通過解算出的三維結(jié)果與真實(shí)GPS數(shù)據(jù)進(jìn)行對(duì)比。GPS點(diǎn)的位置如圖2白點(diǎn)位所示。將解算結(jié)果與GPS數(shù)據(jù)的作差值比較。圖5給出了對(duì)比結(jié)果，曲線越接近于0說(shuō)明解算結(jié)果越準(zhǔn)確。

圖5 結(jié)果對(duì)比

從圖5可以看出，相對(duì)于南北方向來(lái)說(shuō)，東西方向和垂直方向2個(gè)方向的縱坐標(biāo)取值范圍小，這說(shuō)明東西和垂直向的解算精度比南北方向高(縱軸表示與GPS監(jiān)測(cè)出的三維形變的差值)。導(dǎo)致這一結(jié)果的原因有兩點(diǎn)：①衛(wèi)星近南北方向飛行，導(dǎo)致雷達(dá)視線向?qū)δ媳狈较蛐巫儾幻舾?，在解算過程中觀測(cè)值中的噪聲被放大，從而導(dǎo)致南北方向解算精度差[17]；②GPS監(jiān)測(cè)結(jié)果出垂直向精度比水平方向監(jiān)測(cè)精度要差[18]，這會(huì)對(duì)GPS約束改進(jìn)定權(quán)方法解算結(jié)果產(chǎn)生影響。

為了進(jìn)一步說(shuō)明本文提出定權(quán)方法的有效性，表2計(jì)算解算出的三維地表形變與收集的GPS三維地表形變對(duì)比的均方根誤差(RMSE)及統(tǒng)計(jì)了對(duì)比的最大差值。

4 結(jié) 論

本文針對(duì)多視向D-InSAR三維地表形變解算定權(quán)問題，提出了一種基于GPS約束的多視線向D-InSAR三維地表形變監(jiān)測(cè)定權(quán)方法。以2009年拉奎拉地震實(shí)測(cè)數(shù)據(jù)為例進(jìn)行了可行性驗(yàn)證，可以得出以下結(jié)論：

(1) 兩種數(shù)據(jù)監(jiān)測(cè)出的研究區(qū)形變結(jié)果基本一致，但兩種數(shù)據(jù)獲取的時(shí)間間隔、數(shù)據(jù)軌道和波段會(huì)導(dǎo)致監(jiān)測(cè)結(jié)果在形變位置、范圍、形變量等方面略有不同。

表2 對(duì)比結(jié)果定量統(tǒng)計(jì)

(2) 從解算出的三維地表形變分布圖來(lái)看，垂直方向和東西方向的解算結(jié)果相對(duì)較高，南北方向解算結(jié)果受噪聲影響較大。

(3) 相對(duì)于東西方向來(lái)說(shuō)，垂直方向的均方誤差較大，解算精度差一些，這表明GPS低精度的垂直方向形變監(jiān)測(cè)結(jié)果會(huì)對(duì)三維形變解算精度產(chǎn)生影響。

(4) 通過收集到的GPS監(jiān)測(cè)數(shù)據(jù)對(duì)解算結(jié)果進(jìn)行精度分析與驗(yàn)證，對(duì)比發(fā)現(xiàn)解算出的東西方向均方根誤差為0.036 4 m，南北方向均方根誤差為0.513 5 m，垂直方向均方根誤差為0.439 m。

綜上所述，利用本文所提出的GPS約束定權(quán)方法能正確確定多視線向D-InSAR觀測(cè)量的權(quán)，且計(jì)算簡(jiǎn)單。采用GPS約束改進(jìn)定權(quán)方法能提高三維解算精度。