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

王成

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

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

王成?

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

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

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

1 引 言

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

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

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

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

2 實驗區(qū)概況

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

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

3 短基線集DInSAR方法與原理

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

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

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

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

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

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

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

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

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

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

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

圖1 PALSAR時空基線分布圖

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

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

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

圖2 PALSAR基線精化處理

圖3 PALSAR軌道誤差改正

4.3 消除DEM殘余誤差

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

Dv+C·ε=△φ(8)

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

圖4 PALSAR DEM殘余誤差改正量

4.4 大氣相位影響

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

圖5 PALSAR濾波前后對比圖

5 結(jié)果分析

5.1 獲取天津地面相對形變量

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

圖6 天津?qū)嶒瀰^(qū)地面相對形變量變化圖

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

圖8 A點相對于基準(zhǔn)點的沉降趨勢圖

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

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

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

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

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

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

圖10 絕對累積形變量

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

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

圖11 天津?qū)嶒瀰^(qū)DInSAR相干點數(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)附近該沉降漏斗中心的累積形變量達24.6 cm。武清區(qū)的沉降漏斗位于王慶坨鎮(zhèn)附近，累積沉降量達18.5 cm。還有位于北辰區(qū)的西南方向的上河頭鎮(zhèn)附近的沉降漏斗，漏斗中心的沉降量約為21.6 cm。郊區(qū)沉降較大的原因是近幾年天津工業(yè)快速向郊區(qū)轉(zhuǎn)移和發(fā)展，郊區(qū)地下水抽取日益嚴(yán)重，導(dǎo)致地面沉降加速[9，10]。

6 結(jié) 語

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

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

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

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

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

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

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

[2] 王超，張紅，劉智.星載合成孔徑雷達干涉測量[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] 董克剛，王威，于強等.天津地面沉降區(qū)土水比論述[J].水文地質(zhì)工程地質(zhì)，2008，19(3):54～59.

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

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

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

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

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

[10] 張詩玉，李淘，夏耶.基于InSAR技術(shù)的城市地面沉降災(zāi)害監(jiān)測研究[J].武漢大學(xué)學(xué)報·信息科學(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—)，女，工程師，主要從事航空攝影測量與遙感生產(chǎn)工作。

國家科技支撐計劃(2011BAH12B07-03)