基于時序InSAR的東營地區地表沉降研究

劉　鵬，周志偉，汪馳升

(1.深圳大學 信息工程學院,廣東 深圳 518060;2.深圳大學 空間信息智能感知與服務深圳市重點實驗室&海岸帶地理環境監測國家測繪地理信息局重點實驗室，廣東 深圳 518060;3.武漢大學 衛星導航定位技術研究中心，湖北 武漢 430079)

?

基于時序InSAR的東營地區地表沉降研究

劉鵬1,2，周志偉3，汪馳升1,2

(1.深圳大學 信息工程學院,廣東 深圳 518060;2.深圳大學 空間信息智能感知與服務深圳市重點實驗室&海岸帶地理環境監測國家測繪地理信息局重點實驗室，廣東 深圳 518060;3.武漢大學 衛星導航定位技術研究中心，湖北 武漢 430079)

摘要：黃河攜帶大量泥沙入海，河口三角洲覆蓋20～50 m厚的沉積層。沉積層的自然壓實導致該地區的地表沉降。此外，黃河三角洲油氣資源的開發和地下水的開采也加速地表的下陷。InSAR作為一種有效的空間大地測量工具，可以提供幾十公里范圍內的高精度地表形變場。時序InSAR技術在多幅雷達影像組成的干涉網絡基礎上，識別永久散射體(PS)等雷達回波信噪比高的像素。與傳統的InSAR技術相比，時序InSAR技術削弱雷達影像去相干效應的影響，實現長期的形變序列提取。文中利用兩組雷達影像：19幅ERS衛星1992年12月至1996年1月的SAR影像，17幅ENVISAT衛星2003年5月至2008年3月的SAR影像。影像主要覆蓋山東省東營市及其周邊部分區域。結果顯示，東營地區存在每年15 mm以上的地表沉陷，該地區的地表沉降與油氣開采活動有關。

關鍵詞：InSAR；時序InSAR；地表沉降

在海平面上升的趨勢日益明顯的今天，三角洲的沉降已經成為一個世界性的問題。三角洲沉降通常由沉積層的自然壓實、油氣資源和地下水開采等因素造成，并造成一系列的社會、經濟和環境問題[1-2]。山東省東營市地處黃河三角洲，表層沉積物收縮和油氣資源的開采帶來的沉降問題非常明顯。水準測量顯示，1986～1997年間黃河三角洲幾個油田的地表不同程度地下降0.3～0.5 m，凸顯高精度形變監測在該地區的重要性。

InSAR技術通過SAR影像差分干涉測量的相位變化，探測雷達視線方向的地表形變場。該項技術已成功應用于地殼板塊或人為活動引起的地表形變監測，如地震、火山、沉降、滑坡等現象的研究[3-10]。但是InSAR技術的精度受到DEM、軌道、大氣，以及雷達影像失相干噪聲等誤差源的影響[11-12]。時序InSAR技術可以在多幅雷達影像的基礎上，估計這些誤差源的影響，提高形變測量的精度水平。這些時序InSAR技術包括永久散射體技術(PS)、短基線技術(SBAS)、分布式永久散射體技術(SqueeSAR)等不同方法。

國外學者利用時序InSAR技術，分析世界主要三角洲地區的地表沉降，探測到尼羅河三角洲、湄公河三角洲、弗雷澤河三角洲不同成因、不同幅度的地表沉降[13-15]。國內學者針對我國長江、黃河和珠江三角洲地區的沉降也做了大量工作[16-19]。本文針對黃河三角洲東營地區，通過時序InSAR的方法，分析該地區1992～2008年間地表沉降的發展趨勢。

1時序InSAR方法

時序InSAR方法由Ferretti等學者提出，并由Colesanti和Kampes等學者加以改進[20-22]。時序InSAR方法早期主要利用幅度信息識別永久散射體(PS)。基于幅度信息的時序InSAR方法在城市地區效果較好，因為城市存在很多類似于角反射器的人工結構，這些結構可增強雷達波的后向散射，有利于永久散射體的探測。但基于幅度信息的時序InSAR方法在自然地表效果一般，因為自然地表后向散射強度弱，雷達影像去相干效應明顯。相比于幅度分析方法，相位分析方法在自然地表區域提取永久散射體的效果更好。

1.1InSAR的相位模型

本文采用基于相位分析識別永久散射體的斯坦福永久散射體方法(StaMPS)，該方法對形變的時序特征不做假設，利用形變的空間相關性識別PS點[23]。相位分析采用的相位模型為

(1)

地形去除后，干涉圖的相位由形變相位φD、大氣效應φA、軌道誤差φO、視角誤差φθ和隨機噪聲φN幾部分構成。其中形變、大氣和軌道誤差具有空間相關特征，視角誤差具有基線相關特征，因此相位模型可以分解為

(2)

1.2時序InSAR相位分析

相位分析的主要目的是評估永久散射體候選像素的時序相干性。相位分析主要包括以下3個模塊：

1)相位空間相關分析：地表形變、大氣效應和軌道誤差具有空間相關的特征。空間相關相位分析的目的不是為了分離地表形變、大氣效應和軌道誤差，相反，這3種相位被視為單獨一項空間相關相位。空間相關相位是通過條帶濾波從原始相位中分離出來的。

(3)

式中：G為條帶濾波器，條帶濾波器由一個自適應相位濾波器和一個針對相位空間相關距離的低通濾波器組成，濾波后剩余的相位具有空間不相關的特征。

2)相位視角誤差分析：雷達視角的誤差是指像素幾何中心和雷達后向散射相位中心不統一造成的相位變化。

(4)

式中：Δh和Δξ分別是幾何中心與相位中心之間距離的垂直分量和水平分量;θinc是雷達入射角;ρ是衛星與像素之間的距離，θ為視角誤差。對于每個像素來說，兩個中心的距離是一定的，但干涉圖的空間垂直基線是變化的，由此衍生的空間不相關相位也是隨基線變化的。

(5)

其中，B⊥是空間垂直基線。因此，對于每個永久散射體候選像素，利用最小二乘方法，可以估計空間不相關相位與干涉圖空間垂直基線的關聯，反演視角誤差。

(6)

其中，Kmax_γ是視角誤差與基線之間的線性相關系數。

3)相位時序相干性分析：從原始相位中，減去空間相關相位和視角誤差后，殘余相位用來分析相位的時序相干性。殘余相位時序相干性高的像素受大氣、軌道和視角誤差的影響較小，這些像素被選為永久散射體像素。

(7)

其中，ψResidua l,i是殘余相位，γ是時序相干性。

(8)

1.3相位解纏

干涉圖永久散射體像素的原始相位是纏繞在一個2π的整周內，需要通過相位解纏恢復未知的整周相位。StaMPS通過三維解纏恢復永久散射體的整周相位[24]。StaMPS的三維解纏算法集成了SNAPHU二維解纏算法。SNAPHU是一種費用統計網絡流算法，在給定的纏繞相位基礎上，獲取解纏相位的最大后驗概率。

1.4殘余誤差分析

解纏后的展開相位中可能存在長波段大氣信號和軌道誤差的影響。首先，本文通過一個最優擬合的二維平面估計長波段的大氣和軌道誤差；其次，通過大氣濾波對干涉圖中殘留的大氣效應進行估計，大氣濾波的基本原理是大氣效應和形變兩者都具有空間相關性，但是大氣效應在不具有時序相關性，而形變在時序上通常具有連續性。

2SAR數據

本文將斯坦福永久散射體方法(StaMPS)應用到黃河三角洲地區的東營市地表形變監測工作中，分析東營地區在上世紀和本世紀兩個階段的地表形變時序、地表形變速率和特征。

實驗數據包括19幅ERS影像和17幅ASAR影像(見表1和表2)。ERS和Envisat衛星沿132軌道由北向南飛行，雷達側視成像，入射角度均為23°，雷達波長均為C波段5.6 cm。實驗結果包含1992～1996年和2003～2008年兩個時段東營地區的地表形變。

表1　ERS雷達影像的獲取日期和空間垂直基線

表2　ASAR雷達影像的獲取日期和空間垂直基線

3實驗結果

雷達干涉圖的相位相干性與干涉對兩幅SAR影像之間的時間去相干和空間去相干效應有關，為削弱影像的去相干效應，選取空間垂直基線和時間基線分布較為合理影像作為主影像，分別生成18幅ERS干涉圖和16幅ASAR干涉圖(見圖1)。

圖1　干涉網絡的空間與時間基線

ERS影像時序分析提取了21 475個永久散射體像素，這些PS點顯示東營地區的地表形變速率范圍從每年-18.5～12.8 mm(見圖2)。ERS地表形變速率標準差在每年2 mm和3 mm以內的PS

點數目分別達到了總數的74%和99%以上，所以ERS形變速率對時序形變具有代表性。需要注意的是，InSAR作為一種時間和空間雙差分的測量手段，提供的是相對測量結果。由于缺少水準或者GPS等實測地面數據，最初采用影像上所有PS點的平均相位作為相位參考面，然后根據初步速率圖選取無明顯形變的區域作為時序形變的參考區域，重新計算形變和速率。參考區域的選擇還考慮到相位的時序穩定性，避免在參考相位中引入季節性波動的影響。從平均速度場可以發現影像范圍內東營地區在1992～1996年期間存在南北兩個沉降區域。

注：矩形框為參考區域，PQ和ST是兩條剖面的位置，DY1、DY2、DY3和DY4是PS時序點的位置圖2　ERS影像平均速率圖

ASAR影像提取了40 252個永久散射體像素，ASAR提取的永久散射體密度比ERS影像要高。永久散射體密度的提高可能與雷達傳感器的改進和城市發展過程中人工建筑物的增多有關。PS點的速率變化范圍從每年-17.4～8.8 mm(見圖3)。ASAR地表形變速率標準差在2 mm和3 mm每年以內的PS點數目分別達到了總數的95%和99%以上，所以ASAR形變速率對時序形變也具有代表性。ASAR形變分析采用了與ERS同樣范圍的參考區域。從平均速度場可發現影像范圍內東營地區在2003～2008年期間只存在一個沉降區域。

注：矩形框為參考區域，PQ和ST是兩條剖面的位置，DY1、DY2、DY3和DY4是PS時序點的位置圖3　ASAR影像平均速率圖

雷達視線方向單點的形變時序圖顯示該地區永久散射體不同程度的沉降(見圖4)。1號點對應的散射體沉降速率明顯減弱，2號點對應散射體的沉降速率的加強，3號點和4號點顯示散射體的沉降速率略有減弱。

圖4　雷達視線方向單點形變時序

4結束語

北部沉降區在勝利油田勝采三礦范圍。北部沉降區域以黃河路和S316舊省道交叉區域為中心，北至墾利縣S316新省道，南至德州路，東至民豐大道，西至勝利村。北部沉降區包括S316舊省道以北的九戶、左王村和工農村一帶，和S316舊省道以南至鉆井公園一帶。速率剖面PQ顯示ERS影像觀測到的北部沉降在ASAR影像上不明顯(見圖5)。目前還需要確定該地區地表沉降的消失是否與石油開采活動的變化有關。

圖5　剖面PQ形變速率

圖6　剖面ST形變速率

南部沉降區域是勝利油田東辛采油廠范圍。南部沉降區以勝利油田耿井水庫為中心，北至北二路，南至南二路，東至東二路，西至郝純路。剖面PQ顯示南部沉降區耿井水庫周邊的沉降減弱，但剖面ST顯示南部沉降區西部的沉降增強(見圖6)。結合ERS和ASAR影像的形變速率，可以發現東營南部沉降區發生自東向西的移動。

綜上，本文基于不同時段的ERS和ASAR影像，通過時序InSAR分析，在東營地區探測到南北兩個沉降區域，沉降速率在每年18 mm左右。東營地區的觀測地表沉降與油氣資源的開發有關。對比兩個時段的觀測結果發現，南部沉降區域有進一步向東營西部擴展的趨勢，而北部的沉降區域在逐步的消失。

參考文獻：

[1]TORNQVIST T E..Mississippi Delta subsidence primarily caused by compaction of Holocene strata[J].Nature Geosci,2008,1(3)：173-176.

[2]MORTON R,BERNIER J,BARRAS J.Evidence of regional subsidence and associated interior wetland loss induced by hydrocarbon production,Gulf Coast region,USA[J].Environmental Geology,2006,50(2)：261-274.

[3]HU J.Spatial-temporal surface deformation of Los Angeles over 2003—2007 from weighted least squares DInSAR[J].International Journal of Applied Earth Observation and Geoinformation,2013,21(0)：484-492.

[4]SUN Q.Slope deformation prior to Zhouqu,China landslide from InSAR time series analysis[J].Remote Sensing of Environment,2015,156(0)：45-57.

[5]周紅滿,胡金星,柳想,等.基于PS技術的深圳市地面沉降監測[J].測繪工程,2013,22(2):70-74.

[6]胡波，汪漢勝.二軌法DInSAR技術監測城市地表沉降[J].測繪工程,2010，19(2):37-41.

[7]李永生，張景發，李振洪，等.時序InSAR離散相干點相位解纏誤差檢查與校正方法研究[J].武漢大學學報 (信息科學版),2014(10)：1199-1203.

[8]胡俊，李志偉，朱建軍,等.基于BFGS法融合 InSAR和GPS技術監測地表三維形變[J].地球物理學報,2013(1):117-126.

[9]丁榮榮，徐佳，林曉彬.基于PSInSAR技術的常州地表形變監測研究[J].測繪與空間地理信息，2015，38(2)：51-54.

[10] 王宏宇，張慶濤，劉杰，等.利用短基線技術監測地表形變[J].測繪科學，2015，40(10)：123-127.

[11] HANSSEN R F.Radar Interferometry,Data Interpretation and Error Analysis[M].Remote Sensing and Digital Image Processing,ed.F.v.d.Meer.Vol.2.2001,New York：Kluwer Academic Publishers.298.

[12] ZEBKER H AVILLASENOR J.Decorrelation in interferometric radar echoes[J].Geoscience and Remote Sensing,IEEE Transactions on,1992.30(5)：p.950-959.

[13] ALY M H.Permanent Scatterer investigation of land subsidence in Greater Cairo,Egypt[J].Geophysical Journal International,2009.178(3)：p.1238-1245.

[14] ERBAN L E,S.M.Gorelick,and H.A.Zebker,Groundwater extraction,land subsidence,and sea-level rise in the Mekong Delta,Vietnam[J].Environmental Research Letters,2014.9(8)：p.084010.

[15] MAZZOTTI S.Impact of anthropogenic subsidence on relative sea-level rise in the Fraser River delta[J].Geology,2009.37(9)：p.771-774.

[16] DONG S.Time-series analysis of subsidence associated with rapid urbanization in Shanghai,China measured with SBAS InSAR method[J].Environmental Earth Sciences,2014,72(3)：677-691.

[17] WANG H.InSAR reveals coastal subsidence in the Pearl River Delta,China[J].Geophysical Journal International,2012,191(3)：1119-1128.

[18] LIU Y,HUANG H J.Characterization and mechanism of regional land subsidence in the Yellow River Delta,China[J].Natural Hazards,2013,68(2)：687-709.

[19] ZHANG J Z,HUANG H J,BI H B.Land subsidence in the modern Yellow River Delta based on InSAR time series analysis[J].Natural Hazards,2014：1-13.

[20] FERRETTI A,PRATI C,ROCCA F.Permanent Scatterers in SAR interferometry[J].Geoscience and Remote Sensing,IEEE Transactions on,2001,39(1)：8-20.

[21] COLESANTI C.Permanent Scatterers：precision assessment and multi-platform analysis.in Geoscience and Remote Sensing Symposium,2003[A].IGARSS '03.Proceedings.2003 IEEE International.2003.

[22] KAMPES B M.Radar Interferometry,Persistent Scatterer Technique.Remote Sensing and Digital Image Processing[M],ed.F.D.v.d.Meer.2006,12：Springer.

[23] HOOPER A,SEGALL P,ZEBKER H A.Persistent scatterer interferometric synthetic aperture radar for crustal deformation analysis[J].Journal of Geophysical Research,2007,112(B7)：1-21.

[24] HOOPER A,ZEBKER H.Phase unwrapping in three dimensions with application to InSAR time series[J].Optical society of America,2007,24(9)：2737-2747.

[責任編輯：張德福]

Subsidence in Dongying region based on time series InSAR

LIU Peng1,2,ZHOU Zhiwei3,WANG Chisheng1,2

(1.College of Information Engineering,Shenzhen University,Shenzhen 518060,China；2.Shenzhen Key Laboratory of Spatial Smart Sensing and Services &Key Laboratory for Geo-Environment Monitoring of Coastal Zone of the National Administration of Surveying,Mapping and GeoInformation,Shenzhen 518060,China;3.Research Center of GNSS,Wuhan University,Wuhan 430079,China)

Abstract：As its flowing water contains numerous sand and clay,the Yellow River forms a delta with sediments between 20 and 50 meters thick.Subsidence occurs at the Yellow River Delta as a result of the natural compaction of sediments,petroleum and gas extraction together with underground water withdrawal.Interferometric Synthetic Aperture Radar (InSAR) is an advanced geodetic tool with sub-centimetre precision.The InSAR time series technique combines information retrieved from multiple SAR images with a common master scene based on a PS network.The advantage of the PS time series technique indicates that it is more effective in circumventing decorrelation and forming a long-term series of observations than conventional InSAR techniques.Two groups of radar images acquired at different time intervals are used in this study.The time-series results suggest two subsidence areas with a rate of up to 15 mm/yr in Dongying region,which are probably caused by oil extraction.

Key words：InSAR；time series InSAR；subsidence

DOI:10.19349/j.cnki.issn1006-7949.2016.09.001

收稿日期：2015-6-29；修回日期：2015-10-23

基金項目：深圳市科創委項目(ZDSY20121019111146499,JSGG20121026111056204);深圳市戰略性新興產業發展專項資金項目(No.JCYJ20121019111128765);中國博士后基金資助項目(2015M570667,2015M570723)

作者簡介：劉鵬(1986-)，男，博士后.

中圖分類號：TV433

文獻標識碼：A

文章編號：1006-7949(2016)09-0001-05