郭恒 劉志揚 王澤輝

摘要

合成孔徑雷達的差分干涉測量技術(shù)具有精度高、受對流層水蒸氣及氣溶膠影響較小、監(jiān)測范圍廣等特點。文章使用歐洲航天局提供的3期平?jīng)龅貐^(qū)Sentinel-lA雷達波重復(fù)降軌影像，以及相對應(yīng)的21天精密軌道參數(shù)文件，運用InSAR工作流得到的研究區(qū)DEM數(shù)據(jù)，結(jié)合DInSAR兩軌法工作模式，獲取研究區(qū)地表形變量數(shù)據(jù)。在研究期內(nèi)得出，地表抬升變化范圍在0-69.8mm之間，地表沉降變化范圍在0-22.1mm之間，地表的抬升帶位于六盤山構(gòu)造區(qū)南部，秦安一隴縣一寶雞方向，地表沉降漏斗位于甘肅省慶陽市西峰區(qū)，沉降中中心區(qū)面積136.76km²。

【關(guān)鍵詞】合成孔徑雷達 差分干涉測量Sentinel-1A 兩軌法 相位解纏

干涉測量技術(shù)（InSAR）技術(shù)綜合了合成孔徑雷達成像原理與電磁波干涉技術(shù)，運用傳感器的系統(tǒng)參數(shù)、姿態(tài)參數(shù)和軌道之間的幾何關(guān)系等精確測量地表某一點的三維空間位置。差分干涉測量（DInSAR）技術(shù)是InSAR技術(shù)的拓展與延伸，可以精確測量地表形變的微小變化量1989年，Grabriel等人[2]首次論證了DInSAR技術(shù)可運用于探測厘米級的地表形變監(jiān)測1993年，Massonnet等人[3]用ERS-1衛(wèi)星SAR干涉影像數(shù)據(jù)以及該地區(qū)的DEM數(shù)據(jù)，完整地揭示了加利福尼亞蘭德斯地震的同震位移場，得出的形變圖與野外實際監(jiān)測結(jié)果一致程度高。隨后，Zebker等人[4]提出了利用三軌法從SAR干涉圖中提取地震形變的方法。此后，隨著DInSAR技術(shù)的不斷發(fā)展，研究熱點逐漸轉(zhuǎn)移至地面沉降、山體滑坡等細微持續(xù)的地表位移。相比較常規(guī)的GPS測量和水準(zhǔn)測量手段而言，DInSAR技術(shù)可以在大面積范圍內(nèi)監(jiān)測地面的微小形變，不需要人員進入現(xiàn)場區(qū)域測量的特點，可以節(jié)省大量人力成本。其次，DInSAR技術(shù)用幾幅影像就可以監(jiān)測上萬平方公里的地表形變，監(jiān)測效率高，多幅影像連續(xù)的監(jiān)測精度可達毫米級。與傳統(tǒng)的光學(xué)傳感器相比較，合成孔徑雷達發(fā)出的電磁波信號可以很輕松的穿透大氣，受天氣影響微乎其微，更加有利于全天候連續(xù)對地監(jiān)測。

1 D-In SAR監(jiān)測地表沉降原理

DInSAR差分干涉測量技術(shù)是通過兩期或多期相干性強的合成孔徑雷達數(shù)據(jù)進行差分干涉，得到高精度的地表微小形變引起的地表形變相位。合成孔徑雷達干涉測量技術(shù)就是充分利用雷達兩次過境時刻所接收到的回波信號所攜帶的相位信息來得到地表的三維信息或地表形變信息?；驹砭褪峭ㄟ^對同一目標(biāo)的兩個回波信號之間產(chǎn)生了相位差來計算距離的，如圖1表示為通過對雷達數(shù)據(jù)干涉處理后得到形變量、相位及波長之間的關(guān)系。通常情況下，干涉相位包含地形相位、平地相位、大氣延遲相位、噪聲相位及地表形變相位5部分，計算公式為：

δini=δtopo+δdef+δflat+δatmo+δnoise（1）

式中，δtopo為地形起伏引起的地形相位;δdef表示地表形變引起相位;δflat為平地場相位;δatmo表示大氣延遲相位;δnoise為噪聲相位。

2 研究區(qū)地表沉降監(jiān)測

2.1 干涉測量InSAR獲取DEM

InSAR干涉測量獲取DEM是根據(jù)雷達傳感器位置不同，以同一區(qū)域獲取的兩幅測試?yán)走_影像作為數(shù)據(jù)基礎(chǔ)，經(jīng)配準(zhǔn)后生成干涉圖像，再經(jīng)過去平、相位解纏等一系列數(shù)據(jù)處理，從干涉條紋中獲取真實相位信息，進而反演出地形高程數(shù)據(jù)，生成DEM。具體反演技術(shù)流程如圖2。

2.1.1 基線估算

在生成干涉條紋圖前，要用基線做估算的方法來評價干涉圖像對主從影像的質(zhì)量。數(shù)據(jù)為2期成像時間相隔較近的研究區(qū)Sentinel-IA雷達波多視數(shù)據(jù)，成像時間分別為2016年10月6日與2016年10月29日，接收模式IW，極化方式VV，分辨率5*20m。研究區(qū)雷達干涉主從影像的估算結(jié)果顯示，基線長度為26.38m，小于臨界基線6438.17m。雷達波一次2π變換可以探測高程變化588.601m，可監(jiān)測地表形變的敏感度為0.028m。通過估算結(jié)果可以看出，研究所選擇的2期雷達干涉像對的數(shù)據(jù)質(zhì)量達到了干涉測量的技術(shù)要求。

2.1.2 干涉圖生成

雷達接收來自地表的反射電磁波，在反射的重疊區(qū)會有干涉條紋形成，具有周期性。文章通過使用2期Sentinel-1A的干涉寬幅模式數(shù)據(jù)，通過對距離相與方位相參數(shù)調(diào)整，結(jié)合DEM中的投影坐標(biāo)參數(shù)，可以得到配準(zhǔn)后的干涉條紋圖和影像強度圖，如圖3所示。

2.1.3 去平地效應(yīng)

雷達波影像中有平地效應(yīng)，這些區(qū)域的干涉條紋會變得密集，直接影響濾波效果，去除平地效應(yīng)的影影響后使結(jié)果更加精確。去除平地效應(yīng)的方法有優(yōu)勢條紋頻率DEM估算法、成像中心軌道參數(shù)法。文章結(jié)合NASA提供的90m分辨率的STRM4 DEM數(shù)據(jù)對干涉影像做平地效應(yīng)的去除，如圖4。

2.1.4 濾波及相干性計算

為了去除相位噪聲的干擾，需要對影像做濾波處理，如圖4所示。干涉圖通過濾波后，局部放大圖發(fā)現(xiàn)干涉條紋出現(xiàn)的噪聲顯著減少，Data相干系數(shù)在0-1之間，Data越大，表征兩期雷達波影像相干性越高，進而說明兩期雷達波數(shù)據(jù)的地表形變變化在采集期內(nèi)變化不明顯;Data越小，相干程度越低，代表兩期影像的時間周期內(nèi)地表形變發(fā)生了較為明顯的變化，如圖5。

2.1.5 相位解纏

地表的形變量信息的相位是連續(xù)變化纏繞在一起的，干涉條紋經(jīng)過去除平地效應(yīng)和濾波處理后，值的范圍在-π到π之間，加上周期2kπ，k為整數(shù)，讓上升或下降的地表趨勢得到連續(xù)變化，通過相位解纏，干涉條紋可以反映地表真實的相位變化。為保證提取精度，將解纏閾值設(shè)定為0.2，相位解纏后結(jié)果如圖6。

2.1.6 軌道精煉和重去平

軌道精煉過程需要人為選定控制點，選取有效的控制點需要經(jīng)驗與反復(fù)試驗相結(jié)合，一般多選擇地表起伏不大的區(qū)域以及Data相干系數(shù)大的區(qū)域一般會得到較好的結(jié)果。本文運用該方法選取控制點21個，較為平均分布在研究區(qū)內(nèi)，經(jīng)三次多項式校正后得出軌道均方差為30.13m，去除殘差后的平均差為0.05m，去除殘差后的標(biāo)準(zhǔn)殘差為0.35m，得到了比較理想的結(jié)果。

2.1.7 相位轉(zhuǎn)高程與驗證

通過參考相位不斷對合成相位進行校正，逐步完成了相位與高程間的轉(zhuǎn)化，最終形成了精度約為25m并附帶地理編碼的DEM數(shù)據(jù)，投影信息為UTM，橢球參數(shù)為WGS84。如圖7。在ArcGIS中對得到的DEM做渲染。并與Google Earth中對應(yīng)區(qū)域特征點做比對并驗證，發(fā)現(xiàn)兩者吻合度較好

2.2 差分干涉測量DInSAR獲取地表形變數(shù)據(jù)

使用研究區(qū)2期Sentinel-1A雷達重復(fù)降軌影像，見表1，以及歐空局提供的2期數(shù)據(jù)相對應(yīng)的21天精密軌道參數(shù)文件，結(jié)合InSAR工作流得到的研究區(qū)DEM數(shù)據(jù)，運用ENVI5.3與SARscape5.2，使用DInSAR兩軌法工作模式，流程如圖8，得到研究區(qū)地表形變量數(shù)據(jù)，并在ArcGIS中做渲染，如圖9。

3 沉降數(shù)據(jù)分析

根據(jù)差分干涉測量得到的平?jīng)鍪械乇硇巫兛臻g分布圖得出，研究區(qū)在兩期數(shù)據(jù)的時間間隔內(nèi)的地表平均形變量在-22.1mm到69.8mm之間，負號代表地表遠離雷達傳感器運動，正號代表地表朝向雷達傳感器運動。從地表垂直形變量空間分布圖可以看出，近一年來，地表的抬升帶位于六盤山構(gòu)造區(qū)南部，秦安一隴縣一寶雞方向，地表沉降中心位于甘肅省慶陽市西峰區(qū)董志源，沉降中中心區(qū)面積196.76km²。沉降區(qū)域與甘政發(fā)[2017]98號文件公布的地面沉降控制區(qū)范圍基本一致，將2016年西峰區(qū)城市建成區(qū)矢量圖與沉降中心柵格圖在ArcGIS中做空間分析得出，西峰區(qū)城市建成區(qū)的93.56%為沉降區(qū)域，平均沉降量為20.07mm。經(jīng)分析并參考相關(guān)資料，進一步推斷區(qū)內(nèi)的沉降除地質(zhì)條件因素外，大多與地下水開采有密切關(guān)系。

4 結(jié)語

利用Sentinel-1A數(shù)據(jù)，基于差分干涉技術(shù)獲取了平?jīng)龅貐^(qū)2015-2016年的地面沉降量分布狀況，初步探明了研究區(qū)內(nèi)沉降速率與沉降中心，有助于從宏觀方面掌握沉降分布的空間位置，為監(jiān)測因人類活動造成的地表變形提供有效手段。為有效掌握研究區(qū)微小形變及地下水沉降規(guī)律，常規(guī)D-InSAR技術(shù)難以滿足需求，隨著數(shù)據(jù)的進一本積累，后續(xù)將利用多時序InSAR技術(shù)進行監(jiān)測，可以獲得更高的監(jiān)測精度，配合地區(qū)降雨、徑流量變化、城市用水量等數(shù)據(jù)，可以更好的對研究區(qū)的沉降規(guī)律進行整體把握。

參考文獻

[1]張長書，朱建軍，胡俊，王興旺.基于SRTM數(shù)據(jù)的二軌D-InSAR監(jiān)測城市地表沉降研究[J]，測繪科學(xué)，2009，34（02）：45-47.

[2]Gabriel A K，Goldstein R M，Zebker HA.Mapping small elevation changesover large areas：differentialradar interferometry.J Geophys Res94（B7）：9183-9191[J].Journal ofGeophysical Research AtMOSpheres，1989，94（B7）：9183-9191.

[3]Massonnet D，Rossi M，Carmona C，etal.The displacement field of theLanders earthquake mapped by radarinterferometry[J].Nature，1993，364（6433）：138-142.

[4]Zebker H A，Rosen P A，GoldsteinR M，et a1.On the derivation ofcoseismic displacement fields usingdifferential radar interferometry：The Landers earthquake[J].Journalof Geophysical Research Solid Earth，2002，99（B10）：19617-19634.

[5]李廣宇，張瑞，劉國祥，于冰，張波，戴可人，包佳文，韋博文.Sentinel-1A TS-DInSAR京津冀地區(qū)沉降監(jiān)Ni與分析[J].遙感學(xué)報，2018，22（04）：633-646.

[6]成曉倩，馬超，康建榮，鄒友峰.聯(lián)合DInSAR和PIM技術(shù)的沉陷特征模擬和時序分析[J].中國礦業(yè)大學(xué)學(xué)報，2018（05）：1141-1148.

[7]郭恒.基于GIS與RS技術(shù)的平?jīng)鍪械刭|(zhì)災(zāi)害危險性評價[D].四川師范大學(xué)，2017.

[8]陳志謀，胡波，陳金座等.利用InSAR技術(shù)監(jiān)測石油開采引起的地表形變[J].測繪通報，2017（11）：42-46.

[9]劉勇，李培英，豐愛平，黃海軍.黃河三角洲地下水動態(tài)變化及其與地面沉降的關(guān)系[J].地球科學(xué)（中國地質(zhì)大學(xué)學(xué)報），2014，39（11）：1655-1665.

[10]雷坤超，陳蓓蓓，賈三滿，王樹芳，羅勇.基于PS-InSAR技術(shù)的北京地面沉降特征及成因初探[J].光譜學(xué)與光譜分析，2014，34（08）：2185-2189.