王成

(重慶市勘測(cè)院，重慶 400020)

利用PALSAR數(shù)據(jù)提取地面沉降短基線集方法研究

王成?

(重慶市勘測(cè)院，重慶 400020)

短基線集DInSAR方法可以很好地克服數(shù)據(jù)相干性差的困難，本文在分析短基線集方法的基礎(chǔ)上，對(duì)大氣相位和軌道誤差相位的影響進(jìn)行處理得到較為精確的形變相位。采用天津?qū)嶒?yàn)區(qū)從2007年1月17日～2009年10月25日的17景PALSAR單視幅圖像提取地面相對(duì)沉降信息，將DInSAR與水準(zhǔn)測(cè)量數(shù)據(jù)進(jìn)行融合得到外部精度為3.8mm的沉降結(jié)果，同時(shí)確定了天津地面的沉降漏斗的分布和絕對(duì)形變量。

合成孔徑雷達(dá)干涉測(cè)量；短基線集；地面沉降；軌道誤差；大氣相位誤差

1 引 言

合成孔徑雷達(dá)差分干涉測(cè)量技術(shù)(Differential SAR Interferometry，DInSAR)是在InSAR技術(shù)基礎(chǔ)上進(jìn)一步發(fā)展而來的，它通過比較地表形變前后干涉圖之間的相位差來監(jiān)測(cè)地表目標(biāo)的位移[1]。具有全天時(shí)、全天候、廣覆蓋、高分辨率及高地表垂直形變靈敏度等優(yōu)點(diǎn)。但在長(zhǎng)時(shí)間微小地表形變監(jiān)測(cè)中，時(shí)空去相干和大氣相位延遲的影響制約了它的發(fā)展[2]。

Berardino[3]等提出的短基線集算法，可以很好地解決DInSAR長(zhǎng)時(shí)間序列干涉紋圖中缺乏大面積連續(xù)高相干區(qū)域的問題。它是在常規(guī)DInSAR基礎(chǔ)上限制時(shí)空基線的大小后自由組合干涉對(duì)，利用短基線干涉圖提取高相干區(qū)域，從而有效地進(jìn)行形變反演。

本文收集了天津地區(qū)自2007年1月17日～2009年10月25日的16景PALSAR單視復(fù)影像，基于GAMMA軟件平臺(tái)進(jìn)行常規(guī)二軌差分處理得到136組干涉像對(duì)。由于是對(duì)長(zhǎng)時(shí)間地面形變進(jìn)行監(jiān)測(cè)，因此引入了短基線集方法，分別選取其中干涉質(zhì)量好的垂直基線較短的78組干涉像對(duì)采用短基線集DInSAR方法進(jìn)行形變反演，通過基線精化處理和擬合最佳軌道面的方法來削弱軌道誤差影響，采用高斯低通濾波的方法來減弱大氣相位的影響，通過阻尼最小二乘算法分離平均線性形變速率和DEM殘余誤差的改正量，得到最終的垂直形變量。

最后利用短基線集DInSAR方法所得到的垂直形變量與7個(gè)已知水準(zhǔn)點(diǎn)同時(shí)段的沉降量進(jìn)行融合，使得短基線DInSAR方法監(jiān)測(cè)得到的相對(duì)形變量轉(zhuǎn)化為絕對(duì)形變量。PALSAR數(shù)據(jù)監(jiān)測(cè)結(jié)果的精度為3.8 mm，高相干點(diǎn)位密度為666個(gè)/km2。將監(jiān)測(cè)結(jié)果繪制成沉降等值線并疊加到Google Earth地圖上，由此分析確定天津地區(qū)的地面沉降漏斗中心位置分布及其絕對(duì)沉降量。

2 實(shí)驗(yàn)區(qū)概況

地面沉降是許多沿海城市面臨的頑癥，天津地區(qū)的地面沉降現(xiàn)象在中國(guó)具有代表性，其主要原因是地下水過量開采。由于地下流體資源過量開采和構(gòu)造沉降活動(dòng)，在寶坻斷裂和薊運(yùn)河斷裂以南約7 300 km2的廣大地區(qū)，出現(xiàn)不同程度的地面沉降現(xiàn)象，南部和濱海地區(qū)尤為明顯，并且與河北省沉降區(qū)連成一片，在這一大面積的沉降范圍內(nèi)，又形成了市區(qū)、塘沽、漢沽、大港及海河下游工業(yè)區(qū)形成了5個(gè)二級(jí)沉降中心。最近幾年隨著區(qū)縣、鄉(xiāng)鎮(zhèn)經(jīng)濟(jì)的發(fā)展，出現(xiàn)了武清、津南、靜海等新的沉降中心，其發(fā)展速度比當(dāng)年老的沉降區(qū)域還要快，局部沉降速率超過80 mm/y[4]。

常規(guī)地面沉降監(jiān)測(cè)技術(shù)主要有水準(zhǔn)測(cè)量和GPS測(cè)量方法，但二者的空間點(diǎn)位稀疏、周期長(zhǎng)的缺點(diǎn)大大限制其對(duì)區(qū)域地面沉降發(fā)展趨勢(shì)的研究。DInSAR技術(shù)具有大范圍、全天候、高精度、高分辨率地測(cè)量地表三維空間的微小變化，在地表變形的監(jiān)測(cè)方面顯示出前所未有的優(yōu)越性。

3 短基線集DInSAR方法與原理

設(shè)有按時(shí)間序列t0，……，tN獲取的N+1幅單視復(fù)數(shù)SAR影像，將它們以任意影像為主影像進(jìn)行配準(zhǔn)，設(shè)定垂直基線閾值(一般小于400 m)，垂直基線小于該閾值的影像作為一個(gè)基線集，并將其劃分為L(zhǎng)組，每組內(nèi)的影像進(jìn)行差分干涉處理，最終可得到M幅差分干涉圖，且有:

現(xiàn)以t0時(shí)刻為參考時(shí)刻，任意時(shí)刻ti(i=1，…，N)相對(duì)于參考時(shí)刻的差分相位φ(ti)為未知數(shù)，數(shù)據(jù)處理獲取的差分干涉相位δφ(tk)(i=1，…，M)為觀測(cè)量，不考慮大氣延遲相位和數(shù)字高程模型(DEM)誤差的影響，對(duì)第k幅(k=1，……，M)差分干涉相位圖中的任意像元(x，r)可以組成如下方程:

上式中λ為雷達(dá)波長(zhǎng)，d(ti，x，r)和d(tj，x，r)分別為ti和tj時(shí)刻像元(x，r)相對(duì)于參考時(shí)刻t0的雷達(dá)視線方向(line of sight，LOS)地表形變。將式中的φ(ti，x，r)變?yōu)橄噜徲跋裣裨豅OS方向的平均相位變化速率v，有:

其中IE和IS分別為按時(shí)間順序排列的主影像序列和從影像序列，滿足IEj＞ISj，?j=1，……，M化簡(jiǎn)即有:

其中B是M×N矩陣，B為奇異矩陣，可采用阻尼最小二乘算法進(jìn)行求解(Levenberg，1944)，可得到線性形變速率為:

式中γ是阻尼因子，阻尼因子γ越大，分辨就越低；另一方面，當(dāng)γ增大時(shí)分母就越大，從而減小了模型參數(shù)的誤差[5]，阻尼因子γ對(duì)反演結(jié)果有一定影響，本文對(duì)阻尼因子進(jìn)行實(shí)驗(yàn)分析選取了一個(gè)較為合適的數(shù)值，阻尼因子為0.25進(jìn)行阻尼最小二乘的反演計(jì)算。

將式(3)帶入式(2)得:

4 數(shù)據(jù)預(yù)處理與誤差消除

4.1 實(shí)驗(yàn)區(qū)數(shù)據(jù)準(zhǔn)備與預(yù)處理

針對(duì)天津?qū)嶒?yàn)區(qū)，選取了16景ALOSPALSAR從2007年1月17日～2009年10月25日高分辨率模式單視復(fù)影像，所獲取為升軌時(shí)期的數(shù)據(jù)。將其中任意兩幅圖進(jìn)行組合形成120個(gè)像對(duì)，使用GAMMA軟件利用常規(guī)二軌差分干涉測(cè)量方法進(jìn)行處理，對(duì)于PALSAR還需要進(jìn)一步精化基線來削弱線性相位的影響。利用美國(guó)航天飛機(jī)雷達(dá)地形測(cè)繪任務(wù)(SRTM)3″分辨率(90 m)DEM去除地形相位和地平相位，利用最小費(fèi)用流方法對(duì)干涉相位進(jìn)行相位解纏，從而得到差分干涉圖。

處理得到的差分相位由形變相位、殘余軌道誤差相位、DEM殘余誤差相位和大氣相位等組成。由In-SAR原理可知垂直基線長(zhǎng)度與高程模糊度成反比，基線越短，高程模糊度越大，地形對(duì)相位的影響也就越小，所以需要選擇時(shí)空基線較短的干涉像對(duì)組合。通過實(shí)驗(yàn)分析，最終選取了4 000 m以下的干涉質(zhì)量較好的78組干涉像對(duì)，采用短基線集方法進(jìn)行形變反演分析，對(duì)以上誤差相位進(jìn)行分離并予以削弱，得到最終的形變相位。所選取的干涉對(duì)的時(shí)空基線分布如圖1所示。

圖1 PALSAR時(shí)空基線分布圖

4.2 衛(wèi)星軌道誤差的消除

SAR衛(wèi)星軌道誤差是指用來計(jì)算的衛(wèi)星星歷中存在誤差引起的基線計(jì)算不精確而導(dǎo)致后續(xù)的干涉相位誤差。這里首先采取精化基線的方法對(duì)軌道誤差進(jìn)行一定的改正，如圖2所示。為了進(jìn)一步削弱殘余軌道誤差的影響，然后通過線性擬合的方法來對(duì)解纏后的干涉相位進(jìn)行擬合，得到最佳擬合軌道面，即殘余的軌道誤差，并將其從干涉圖中去除。按照如下方程來計(jì)算每個(gè)干涉圖的最佳的軌道面[6]:

式中，[x，y]是每個(gè)像素點(diǎn)的在像素坐標(biāo)系下的坐標(biāo)，u，v是斜率參數(shù)，w是截距，如圖3所示。

圖2 PALSAR基線精化處理

圖3 PALSAR軌道誤差改正

4.3 消除DEM殘余誤差

此外，在實(shí)際處理過程中需要考慮DEM殘余誤差的影響，則式(5)表示為:

Dv+C·ε=△φ(8)

式中cT=[(4π/λ)(B⊥1/r sinθ)，…，(4π/λ)(B⊥N/ r sinθ)]，λ為雷達(dá)波長(zhǎng)，B⊥為垂直基線，r是衛(wèi)星到地面間的距離，θ是衛(wèi)星視角。通過阻尼最小二乘算法計(jì)算出DEM殘余誤差的改正量ε，其大小在[-10，10]m范圍內(nèi)，結(jié)果如圖4所示。同時(shí)可得到沿雷達(dá)視線方向的線性形變速率v乘以時(shí)間就得到雷達(dá)視線方向的形變量再除以cosθ，從而轉(zhuǎn)化為垂直形變量。

圖4 PALSAR DEM殘余誤差改正量

4.4 大氣相位影響

短基線集方法是采用線性模型，為了進(jìn)一步消去殘余相位中的大氣相位和非線性形變相位的影響，可以在線性模型的基礎(chǔ)上進(jìn)行適當(dāng)?shù)臑V波處理[7]。根據(jù)大氣相位空域的低頻特性，本文采用在空間上實(shí)施高斯低通濾波處理，移除“形變信號(hào)”中上的高頻部分。如圖5所示，是對(duì)其中一組數(shù)據(jù)進(jìn)行濾波的前后對(duì)比圖。

圖5 PALSAR濾波前后對(duì)比圖

5 結(jié)果分析

5.1 獲取天津地面相對(duì)形變量

此次形變反演過程中相干閾值設(shè)置為0.67，得到高相干點(diǎn)個(gè)數(shù)為599 728，影像面積約為900 km2，即點(diǎn)密度為666個(gè)/km2。DInSAR測(cè)量結(jié)果是相對(duì)于基準(zhǔn)點(diǎn)而言的，該基準(zhǔn)點(diǎn)是研究區(qū)內(nèi)所有其他點(diǎn)形變量的參考基準(zhǔn)，所以需要選擇形變較為穩(wěn)定的點(diǎn)作為基準(zhǔn)點(diǎn)。同時(shí)選擇干涉質(zhì)量較好的2007年7月20日主影像為沉降基準(zhǔn)時(shí)刻，將另外的16景影像與之差分，從而得到16幅天津?qū)嶒?yàn)區(qū)地面相對(duì)形變量變化圖，如圖6所示。圖中“+”表示的參考基準(zhǔn)點(diǎn)，是研究區(qū)內(nèi)相對(duì)穩(wěn)定的點(diǎn)，負(fù)值表示沉降。通過地理編碼，將沉降結(jié)果由制圖坐標(biāo)系轉(zhuǎn)換到WGS84坐標(biāo)下，可以獲取實(shí)驗(yàn)區(qū)內(nèi)自2007年01月17日～2009年10月25日的累積相對(duì)形變量圖，如圖7所示。圖中“★”表示參考基準(zhǔn)點(diǎn)，“?”表示累積沉降量最大的點(diǎn)，我們定義為A點(diǎn)，其相對(duì)形變量不超過15 cm，將其相對(duì)于參考基準(zhǔn)點(diǎn)的沉降趨勢(shì)表示出來，如圖8所示。

圖6 天津?qū)嶒?yàn)區(qū)地面相對(duì)形變量變化圖

圖7 2007年01月17日～2009年10月25日天津地面相對(duì)形變量累積圖

圖8 A點(diǎn)相對(duì)于基準(zhǔn)點(diǎn)的沉降趨勢(shì)圖

5.2 DInSAR與水準(zhǔn)數(shù)據(jù)的融合

通過短基線集DInSAR方法獲得的是大范圍相對(duì)形變量，而通過水準(zhǔn)測(cè)量可以獲得高精度的局部范圍的絕對(duì)沉降信息。為了獲取大范圍的高精度地面沉降信息，則可以將二者進(jìn)行數(shù)據(jù)融合。由于DInSAR測(cè)量的基準(zhǔn)與水準(zhǔn)測(cè)量基準(zhǔn)不一致，所以需要進(jìn)行基準(zhǔn)的統(tǒng)一。通過數(shù)據(jù)分析，兩者監(jiān)測(cè)結(jié)果關(guān)系可采用下列函數(shù)表達(dá)[8]，令:

式中y為水準(zhǔn)測(cè)量結(jié)果，x為雷達(dá)干涉測(cè)量結(jié)果。

基于2007.01～2009.10的水準(zhǔn)測(cè)量結(jié)果和2007.01.17～2009.10.25的DInSAR結(jié)果(由于水準(zhǔn)測(cè)量的周期較長(zhǎng)，所以不能保證與DInSAR數(shù)據(jù)的時(shí)間段完全重合)，對(duì)水準(zhǔn)與DInSAR結(jié)果進(jìn)行了融合。首先按線性關(guān)系對(duì)兩者進(jìn)行了擬合分析，剔除偏差較大的點(diǎn)后，再用7個(gè)點(diǎn)分布均勻的數(shù)據(jù)對(duì)兩者重新擬合，得到新的函數(shù)關(guān)系結(jié)果如圖9所示。

圖9 DInSAR數(shù)據(jù)與水準(zhǔn)數(shù)據(jù)的線性關(guān)系圖

用于擬合的點(diǎn)的校正后的DInSAR結(jié)果與水準(zhǔn)數(shù)據(jù)的標(biāo)準(zhǔn)差為3.8 mm。則可利用水準(zhǔn)與DInSAR之間的函數(shù)關(guān)系對(duì)DInSAR形變結(jié)果進(jìn)行校正，可以得到整個(gè)區(qū)域的經(jīng)過水準(zhǔn)校正的地面沉降成果，從而獲取天津?qū)嶒?yàn)區(qū)的從2007年01月17日～2009年10月25日的絕對(duì)形變量，結(jié)果如圖10所示。5.3 實(shí)驗(yàn)區(qū)沉降結(jié)果分析

圖10 絕對(duì)累積形變量

將經(jīng)過水準(zhǔn)數(shù)據(jù)修正之后的InSAR數(shù)據(jù)疊加到GoogleEarth地圖上，各相干點(diǎn)目標(biāo)沉降累積形變量的大小以不同的顏色分級(jí)表示，結(jié)果顯示DInSAR數(shù)據(jù)點(diǎn)主要分布在天津市中心的范圍，越靠近郊區(qū)的范圍的相干點(diǎn)數(shù)據(jù)越為稀疏(主要為農(nóng)田或魚塘)，如圖11所示。將沉降結(jié)果繪制成沉降等值線疊加到GoogleE-arth地圖上，如圖12所示。

從圖5中沉降趨勢(shì)分析，自2007年1月～2009年10月天津市主城區(qū)沉降速率較為穩(wěn)定，累積沉降量約為4 cm左右。沉降速率從中心到郊區(qū)逐漸變大，說明主城區(qū)的地面沉降已經(jīng)得到了有效的遏制，而郊區(qū)的沉降速率仍然偏大，且發(fā)現(xiàn)局部地區(qū)存在明顯的沉降漏斗。

圖11 天津?qū)嶒?yàn)區(qū)DInSAR相干點(diǎn)數(shù)據(jù)疊加到Google Earth圖

圖12 DInSAR沉降等值線疊加到Google Earth圖

從圖11和圖12可以看出沉降漏斗主要分布在武清區(qū)、北辰區(qū)以及西青區(qū)。研究區(qū)域內(nèi)最大(負(fù)值代表下沉)沉降漏斗中心位于西青區(qū)南河鎮(zhèn)附近該沉降漏斗中心的累積沉降量約為26.6 cm和大寺鎮(zhèn)附近該沉降漏斗中心的累積形變量達(dá)24.6 cm。武清區(qū)的沉降漏斗位于王慶坨鎮(zhèn)附近，累積沉降量達(dá)18.5 cm。還有位于北辰區(qū)的西南方向的上河頭鎮(zhèn)附近的沉降漏斗，漏斗中心的沉降量約為21.6 cm。郊區(qū)沉降較大的原因是近幾年天津工業(yè)快速向郊區(qū)轉(zhuǎn)移和發(fā)展，郊區(qū)地下水抽取日益嚴(yán)重，導(dǎo)致地面沉降加速[9，10]。

6 結(jié) 語

本文針對(duì)天津地區(qū)地面沉降監(jiān)測(cè)試驗(yàn)主要做了以下幾項(xiàng)工作，首先是利用16景PALSAR單幅視影像進(jìn)行常規(guī)二通差分干涉測(cè)量方法處理，采用精軌數(shù)據(jù)去除較大的軌道誤差，利用SRTM DEM去除地形相位和地平相位的影響，利用最小費(fèi)用流方進(jìn)行相位解纏，從而得到差分干涉圖。然后根據(jù)短基線集方法的原理，選取78個(gè)干涉對(duì)進(jìn)行時(shí)間序列分析，采用擬合最佳軌道面的方法消除殘余軌道誤差的影響，并利用阻尼最小二乘算法求解DEM殘余誤差改正量以及所需要得到的平均線性形變速率。為了減弱大氣相位的影響，這里做了空間低通濾波的處理，得到天津地區(qū)2007年1月～2009年10月的相對(duì)累積形變。最后為了獲得絕對(duì)形變量，與7個(gè)水準(zhǔn)點(diǎn)測(cè)量數(shù)據(jù)進(jìn)行了融合，獲得外部精度為3.8 mm的沉降結(jié)果，并確定天津地區(qū)沉降漏斗的分布及其絕對(duì)沉降量。

基于上述工作得到了以下幾點(diǎn)認(rèn)識(shí):

(1)短基線集方法算法較PS方法簡(jiǎn)單，本文的特點(diǎn)在于在該算法基礎(chǔ)上對(duì)軌道誤差和大氣相位的影響做了一定的分析和削弱；

(2)在利用阻尼最小二乘算法進(jìn)行反演過程中，選擇較為合適的阻尼因子是關(guān)鍵，一般選擇在0.01～1之間；

(3)在水準(zhǔn)數(shù)據(jù)與DInSAR數(shù)據(jù)融合過程中，由于水準(zhǔn)點(diǎn)處沒有高相干點(diǎn)數(shù)據(jù)，其水準(zhǔn)點(diǎn)的DInSAR數(shù)據(jù)是通過插值得到，對(duì)最終結(jié)果有一定影響；

(4)利用水準(zhǔn)數(shù)據(jù)可以將DInSAR的相對(duì)沉降信息轉(zhuǎn)為絕對(duì)沉降信息，同時(shí)兼顧水準(zhǔn)測(cè)量高精度的優(yōu)點(diǎn)和DInSAR大面積、分辨率高的優(yōu)點(diǎn)，能夠在地面沉降監(jiān)測(cè)的實(shí)際應(yīng)用工作中發(fā)揮一定作用。

[1] 周春霞.星載SAR干涉測(cè)量技術(shù)及其在南極冰貌地形研究中的應(yīng)用[D].武漢:武漢大學(xué)，2004.

[2] 王超，張紅，劉智.星載合成孔徑雷達(dá)干涉測(cè)量[M].北京:科學(xué)出版社，2002.

[3] Berardino，P.Fornaro，G.Lanari，R.et al.A new algorithm for surface deformationmonitoring based on smallbaseline differential interferograms[J].IEEE Transactions on Geoscience and Remote Sensing，Vol.40，No.11，2375～2383，2002.

[4] 董克剛，王威，于強(qiáng)等.天津地面沉降區(qū)土水比論述[J].水文地質(zhì)工程地質(zhì)，2008，19(3):54～59.

[5] 崔希璋，於宗儔，陶本藻等.廣義測(cè)量平差(第二版)[M].武漢:武漢大學(xué)出版社，2009.

[6] 湯益先.基于相干目標(biāo)的DInSAR方法及其地表沉降應(yīng)用研究[D].北京:中國(guó)科學(xué)研究院，2006.

[7] 張紅，王超，吳濤.基于相干目標(biāo)的DInSAR方法研究[M].北京:科學(xué)出版社，2008.

[8] 張?jiān)娪?干涉圖大氣效應(yīng)校正方法研究及InSAR在天津地面沉降監(jiān)測(cè)中的應(yīng)用[D].武漢:武漢大學(xué)，2009.

[9] 范景輝，郭華東，郭小方等.基于相干目標(biāo)的干涉圖疊加方法監(jiān)測(cè)天津地區(qū)地面沉降[J].遙感學(xué)報(bào)，2008，12 (1):111～118.

[10] 張?jiān)娪?，李淘，夏?基于InSAR技術(shù)的城市地面沉降災(zāi)害監(jiān)測(cè)研究[J].武漢大學(xué)學(xué)報(bào)·信息科學(xué)版.2008，33 (8):850～853.

Study on Using PALSAR Data Extract Land Subsidence w ith Short Baseline Subset M ethod

Wang Cheng
(Chongqing Survey Insitute，Chongqing 400020，China)

The short baseline subset DInSAR technology can be very good to overcome the difficulties of poor coherence data，this paper is based on the analysis of shortbaseline subsetmethod toweak the atmospheric phase and orbital phase error，then gainmore accurate deformation phase.This paper collected 17 ALOSPALSAR single complex images since the January 17，2007 to October 2009，25which cover the Tianjin area，Intergration of DInSAR and leveling data to get external precision as the settlement results of 3.8mm.Finally，Tianjin ground subsidence funnel distribution and absolute variables are determined.

DInSAR；short baseline subset；groud subsidence；orbit error；atmospheric phase error

1672-8262(2013)05-81-06

P237

A

2013—03—18

王成(1988—)，女，工程師，主要從事航空攝影測(cè)量與遙感生產(chǎn)工作。

國(guó)家科技支撐計(jì)劃(2011BAH12B07-03)