劉　鵬，周志偉，汪馳升

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

?

基于時序InSAR的東營地區(qū)地表沉降研究

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

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

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

關(guān)鍵詞：InSAR；時序InSAR；地表沉降

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

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

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

1時序InSAR方法

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

1.1InSAR的相位模型

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

(1)

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

(2)

1.2時序InSAR相位分析

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

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

(3)

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

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

(4)

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

(5)

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

(6)

其中，Kmax_γ是視角誤差與基線之間的線性相關(guān)系數(shù)。

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

(7)

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

(8)

1.3相位解纏

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

1.4殘余誤差分析

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

2SAR數(shù)據(jù)

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

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

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

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

3實驗結(jié)果

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

圖1　干涉網(wǎng)絡(luò)的空間與時間基線

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

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

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

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

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

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

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

4結(jié)束語

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

圖5　剖面PQ形變速率

圖6　剖面ST形變速率

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

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

參考文獻：

[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技術(shù)的深圳市地面沉降監(jiān)測[J].測繪工程,2013,22(2):70-74.

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

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

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

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

[10] 王宏宇，張慶濤，劉杰，等.利用短基線技術(shù)監(jiān)測地表形變[J].測繪科學(xué)，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.

[責(zé)任編輯：張德福]

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

基金項目：深圳市科創(chuàng)委項目(ZDSY20121019111146499,JSGG20121026111056204);深圳市戰(zhàn)略性新興產(chǎn)業(yè)發(fā)展專項資金項目(No.JCYJ20121019111128765);中國博士后基金資助項目(2015M570667,2015M570723)

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

中圖分類號：TV433

文獻標識碼：A

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