Sentinel-1A京津冀平原區(qū)2016—2018年地面沉降InSAR監(jiān)測

史 珉， 宮輝力， 陳蓓蓓， 高明亮， 張舜康

(1.首都師范大學(xué)地面沉降機(jī)理與防控教育部重點(diǎn)實(shí)驗(yàn)室，北京 100048； 2.首都師范大學(xué)水資源安全北京實(shí)驗(yàn)室，北京 100048； 3.首都師范大學(xué)北京市城市環(huán)境過程與數(shù)字模擬國家重點(diǎn)實(shí)驗(yàn)室培育基地，北京 100048； 4.首都師范大學(xué)三維信息獲取與應(yīng)用教育部重點(diǎn)實(shí)驗(yàn)室，北京 100048)

0 引言

地面沉降是一種地表高程下降的環(huán)境地質(zhì)現(xiàn)象，其發(fā)生發(fā)展增加了城市內(nèi)澇、海水倒灌等風(fēng)險(xiǎn)，并可以誘發(fā)地面塌陷、地裂縫等系列環(huán)境災(zāi)害，形成災(zāi)害鏈。區(qū)域地面沉降不僅對城市基礎(chǔ)設(shè)施、高速鐵路等重大工程產(chǎn)生潛在威脅，還影響著城市地下空間資源利用效率，成為制約區(qū)域可持續(xù)發(fā)展的重大問題。

京津冀地區(qū)所處的華北平原是我國沉降速率最大、影響面積最大的地區(qū)[1]。自1923年首次在天津發(fā)現(xiàn)地面沉降以來[2]，華北平原地面沉降按期發(fā)展過程可分為4個(gè)階段： 局部孕育階段(1954年之前)、加速發(fā)展階段(1955—1974年)、急速擴(kuò)張階段(1975—1984年)和控制減緩階段(1985年至今)[3]。截至2012年，華北平原沉降速率大于20 mm/a的區(qū)域面積約占其總面積的12.5%，其中90%的沉降區(qū)位于北京、天津、河北三地，局部地區(qū)累計(jì)沉降量超過3.4 m[4]。京津冀協(xié)同發(fā)展地質(zhì)調(diào)查中將地面沉降列為京津冀一體化進(jìn)程中需要關(guān)注的重大地質(zhì)問題之一。對其地面沉降進(jìn)行監(jiān)測，認(rèn)識其分布規(guī)律、演化特征，對京津冀城市和重大工程規(guī)劃建設(shè)，為實(shí)現(xiàn)京津冀協(xié)同發(fā)展的重大國家戰(zhàn)略提供地質(zhì)環(huán)境保障。

合成孔徑雷達(dá)干涉測量技術(shù)(interferometric synthetic aperture Radar，InSAR)已被廣泛應(yīng)用于地表形變監(jiān)測中。與水準(zhǔn)測量、分層標(biāo)、全球定位系統(tǒng)(global positioning system，GPS)測量技術(shù)相比，InSAR技術(shù)具備大范圍、全天時(shí)全天候獲取高時(shí)空分辨率地表形變的能力。多時(shí)相InSAR(multi-temporal InSAR，MT-InSAR)技術(shù)的發(fā)展(如永久散射體干涉(persistent scatterer InSAR，PS-InSAR)、小基線 (small baseline subsets，SBAS)等)更是克服了大氣延遲、時(shí)空失相干對傳統(tǒng)InSAR方法的影響[5-6]。目前，MT-InSAR技術(shù)已被廣泛應(yīng)用于京津冀地面沉降研究中，但這類研究多關(guān)注北京[7-11]、天津[12-15]、滄州[16-18]和廊坊[19-20]等城市區(qū)域，通過多源SAR數(shù)據(jù)獲取其地面沉降信息，而對京津冀平原地面沉降的分布、范圍及沉降程度缺乏整體性認(rèn)識。針對單軌SAR數(shù)據(jù)源幅寬對區(qū)域地面沉降監(jiān)測范圍的限制問題，本研究基于多軌SAR數(shù)據(jù)融合技術(shù)，獲取了2016—2018年京津冀平原地面沉降信息； 為保證InSAR監(jiān)測結(jié)果的可靠性，選用同時(shí)期水準(zhǔn)測量數(shù)據(jù)對InSAR結(jié)果進(jìn)行精度驗(yàn)證，并對相鄰軌道結(jié)果進(jìn)行一致性檢驗(yàn)； 并結(jié)合GIS空間分析、剖面分析等方法，獲取InSAR監(jiān)測時(shí)段內(nèi)京津冀平原區(qū)地面沉降演化特征。

1 研究區(qū)概況及其SAR數(shù)據(jù)源

京津冀地區(qū)(N36°05′～42°37′，E113°11′～119°45′)位于我國華北地區(qū)，包含北京、天津和河北省的11個(gè)地級市，屬暖溫帶大陸性季風(fēng)型氣候，降水多集中在7—8月。該地區(qū)地形條件復(fù)雜，總體地勢為西北高、東南低，自西北向東南依次分布壩上高原、燕山—太行山山地和東南部平原[21]。京津冀地區(qū)存在多種環(huán)境地質(zhì)問題，其中以地面沉降等為例的緩變性地質(zhì)主要發(fā)生在平原區(qū)。京津冀平原區(qū)是華北平原的重要組成部分，自山麓到海岸，按成因和形態(tài)特征可分為第四紀(jì)山前沖積、洪積傾斜平原(山前平原)，中部沖積、湖積多層疊加平原(中部平原)和東部沖積、湖積夾海積的濱海平原(濱海平原)[22]。截至2015年，京津冀平原區(qū)累積沉降量大于200 mm的面積約為6.4萬km2[21]，約占平原區(qū)面積的73%，沉降速率大于50 mm/a的嚴(yán)重沉降區(qū)面積占全國嚴(yán)重沉降區(qū)面積的92%(http: //www.tianjin.cgs.gov.cn)。

Sentinel-1A衛(wèi)星是歐洲空間局于2014年發(fā)射的重訪周期為12 d的C波段對地觀測衛(wèi)星。該衛(wèi)星使用近極地太陽同步軌道，軌道高度、傾角和周期分別為693 km，98.18°和99 min，具多種成像模式和極化方式(https: //sentinel.esa.int/)。單軌SAR影像受其幅寬限制，監(jiān)測范圍有限，為獲取京津冀平原區(qū)的地面沉降信息，本次研究共選用了3軌時(shí)間跨度為2016—2018年的Sentinel-1A數(shù)據(jù)(Track 40，Track 142和Track 69)。所用3軌Sentinel-1A數(shù)據(jù)均為寬幅干涉模式(interferometric wide，IW)獲取的單視復(fù)數(shù)(single look complex，SLC)升軌存檔影像，其極化方式均為垂直極化(vertical-vertical，VV)，空間分辨率為5 m×20 m，單幅幅寬達(dá)250 km。影像覆蓋范圍及主要參數(shù)分別如圖1和表1所示。

圖1 研究區(qū)概況示意圖Fig.1 Location of study area

表1 Sentinel-1A影像信息Tab.1 Satellite information for the data used in this study

2 研究方法

2.1 MT-InSAR處理過程

通過SARPROZ軟件(https: //www.sarproz.com/)對3軌Sentinel-1A影像分別進(jìn)行PS-InSAR處理，以獲取每軌SAR影像范圍內(nèi)時(shí)序地面沉降信息。處理過程包括以下幾個(gè)主要步驟： 主影像選取，數(shù)據(jù)集配準(zhǔn)及重采樣，基線建立，數(shù)字高程模型(digital elevation model,DEM)模擬，生成差分干涉圖，選取PS候選點(diǎn)，相位解纏，大氣相位估計(jì)和去除，重新選取PS點(diǎn)和估計(jì)PS點(diǎn)形變速率[23-24]。由前人研究成果已知，京津冀地區(qū)水平向形變值約為1～3 mm/a[25-28]，因此在本次研究中忽略水平向形變結(jié)果，通過式(1)將視線向(line-of-sight，LOS)InSAR形變結(jié)果轉(zhuǎn)為垂向，并將經(jīng)過處理后的垂向結(jié)果用于與水準(zhǔn)結(jié)果的精度驗(yàn)證和后續(xù)的SAR影像結(jié)果拼接中，即

dv=dLOS/cosθ

(1)

式中：dv為垂向形變結(jié)果，mm；dLOS為LOS向形變結(jié)果，mm；θ為雷達(dá)LOS向與垂向的夾角。

2.2 多軌SAR數(shù)據(jù)拼接

為了獲取大區(qū)域的地面形變信息，需對上文所得單軌SAR影像結(jié)果進(jìn)行拼接。首先對各軌SAR影像垂向形變結(jié)果進(jìn)行了一致性檢驗(yàn)和精度驗(yàn)證，在保證數(shù)據(jù)精度可靠的情況下，對相鄰軌道數(shù)據(jù)進(jìn)行融合。在相鄰軌道沉降速率融合中，Ge等[29]通過引入與參考點(diǎn)距離權(quán)重的關(guān)系進(jìn)行整體平差; 孫赫[30]通過影像重疊區(qū)內(nèi)PS點(diǎn)形變量差值眾數(shù)對相鄰軌道數(shù)據(jù)結(jié)果進(jìn)行校正; 熊思婷等[31]比較了區(qū)塊法和插值法在求解異軌重疊區(qū)的形變差中的應(yīng)用; 張永紅等[32]通過大量的水準(zhǔn)觀測數(shù)據(jù)對SAR數(shù)據(jù)結(jié)果進(jìn)行了平差控制。本次實(shí)驗(yàn)中，依據(jù)重疊區(qū)內(nèi)同名高相干點(diǎn)的形變速率一致準(zhǔn)則[32-33]，計(jì)算平差權(quán)重。

由于受不同軌道衛(wèi)星成像幾何差異，不同軌道SAR數(shù)據(jù)所得LOS向地表形變結(jié)果間存在系統(tǒng)性偏差。為消除不同視角投影的影響，首先依據(jù)式(1)將InSAR監(jiān)測所得各軌LOS向形變值轉(zhuǎn)為垂向形變值。如圖1所示，Track 142與Track 69和40間都具有重疊區(qū)，因此在多軌SAR數(shù)據(jù)融合時(shí)，將Track 142形變結(jié)果作為基準(zhǔn)。Track 142與Track 40影像重疊區(qū)域內(nèi)，雷達(dá)LOS向與垂向的夾角分別為33.72° 和43.75°； Track 142與Track 69影像重疊區(qū)內(nèi)LOS向與垂向夾角分別為43.77°和33.67°。根據(jù)上述θ值，將重疊區(qū)域內(nèi)每軌影像的各PS點(diǎn)LOS向結(jié)果依據(jù)式(1)分別投影到垂直向上，選取的PS點(diǎn)的相干性均在0.9以上。之后，選用最鄰近法對重疊區(qū)內(nèi)各軌PS點(diǎn)進(jìn)行同名點(diǎn)匹配。最后，采用最小二乘法(least square method, LSM)構(gòu)建相鄰軌道重疊區(qū)域的平差方程，并分別對Track 40和Track 69垂向形變結(jié)果進(jìn)行偏移量補(bǔ)償。該方法是通過實(shí)際觀測值與模型計(jì)算值的誤差的平方和最小來尋找最優(yōu)的匹配函數(shù)。假設(shè)相鄰軌道分別為x和y，其中x為參考基準(zhǔn)，y為需要進(jìn)行平差的數(shù)據(jù)，相鄰軌道之間存在n個(gè)重合點(diǎn)，基于LSM的原理，則

(2)

此時(shí)，通過對δ求極限，即可求得a和b的值。最終，獲得多軌SAR數(shù)據(jù)融合結(jié)果。

3 結(jié)果及精度驗(yàn)證

3.1 京津冀2016—2018年地面沉降監(jiān)測結(jié)果分析

為獲取京津冀平原區(qū)地面沉降信息，基于MT-InSAR技術(shù)分別獲取了3軌Sentinel-1A數(shù)據(jù)源地面沉降信息，后對相鄰軌道數(shù)據(jù)進(jìn)行融合獲取了京津冀平原區(qū)內(nèi)2016—2018年地面沉降速率(圖2)。為了顯示清晰，對結(jié)果進(jìn)行了柵格化處理，像元大小為1 km。InSAR監(jiān)測結(jié)果顯示， 2016—2018年間京津冀平原區(qū)內(nèi)最大沉降速率達(dá)164 mm/a，分布在北京朝陽—通州沉降區(qū)和廊坊—天津沉降區(qū)內(nèi)。研究區(qū)內(nèi)沉降分布范圍廣，且空間分布不均，在北京、天津及河北省各地級市均有不同程度的地面沉降問題，且三地沉降區(qū)有連片趨勢。其中北京平原區(qū)沉降較為嚴(yán)重的地區(qū)主要分布在北京東南部朝陽、通州等地，且北京東南部沉降區(qū)已擴(kuò)張到廊坊三河市界內(nèi)。天津、廊坊兩地沉降較為嚴(yán)重的區(qū)域分布在兩地交界處。河北平原沉降較為嚴(yán)重的區(qū)域主要分布在其中東部地區(qū)，且該地區(qū)沉降區(qū)呈南北貫通趨勢。研究區(qū)內(nèi)共獲得了1 426 967個(gè)PS點(diǎn)，其中2016—2018年沉降速率大于20, 50, 100 mm/a的點(diǎn)數(shù)分別占總點(diǎn)數(shù)的53.9%, 23.9%和3.6%。沉降速率大于20, 50, 100 mm/a的面積分別為3.3萬, 1.4萬和0.14萬km2。據(jù)統(tǒng)計(jì)，2016年沉降速率大于20, 50, 100 mm/a的點(diǎn)數(shù)分別占總點(diǎn)數(shù)的53.8%, 23.7%和3.6%； 到2017年沉降速率大于20和50 mm/a的點(diǎn)數(shù)略有增長，分別占總點(diǎn)數(shù)的54.4%和23.8%，沉降速率大于100 mm/a的點(diǎn)數(shù)占比則降至3.5%。2018年沉降速率大于20和50 mm/a的PS點(diǎn)數(shù)與2017年比有所增長，占比達(dá)57.2%和24.0%，沉降速率大于100 mm/a的點(diǎn)數(shù)與2017年相比保持穩(wěn)定。

圖2 2016—2018年京津冀平原區(qū)PS-InSAR地表形變速率Fig.2 Mean displacement velocities throughout the Beijing-Tianjin-Hebei derived from the Sentinel-1A data by using PS-InSAR from 2016 to 2018

研究區(qū)內(nèi)地面沉降大于50 mm的區(qū)域主要分布在北京東南部，廊坊北部，廊坊東部到天津西部，天津南部，保定東部到廊坊西部，保定西部，唐山—秦皇島沿海區(qū)域，以及河北平原東南部區(qū)域。對監(jiān)測時(shí)段內(nèi)研究區(qū)各子平原內(nèi)沉降速率大于50 mm/a的嚴(yán)重沉降區(qū)面積進(jìn)行統(tǒng)計(jì)發(fā)現(xiàn)，山前沖積洪積平原沉降量大于50 mm的面積有輕微增長(2016—2018年)，由2 937 km2變?yōu)? 982 km2，漲幅為1.5%。2016—2018年中部平原沉降量大于50 mm的面積的漲幅為0.3%，基本維持穩(wěn)定。2016—2018年，嚴(yán)重沉降區(qū)的面積分別為10 547 km2，10 574 km2和10 577 km2。與2016年相比，濱海平原嚴(yán)重沉降區(qū)的面積在2018年增長較為明顯，由382 km2變?yōu)?77 km2，漲幅達(dá)33.8%。其中唐山—秦皇島沿海沉降區(qū)范圍明顯擴(kuò)大(如圖3所示)。為進(jìn)一步分析研究時(shí)段內(nèi)地面沉降發(fā)展過程，對沉降較嚴(yán)重區(qū)

(a) 2016年 (b) 2017年 (c) 2018年

圖3 2016—2018年各年度InSAR平均形變速率Fig.3 Mean displacement velocities derived by InSAR between 2016 and 2018

域沿東西、南北方向分別進(jìn)行剖面分析，剖面線位置如圖3(c)中虛線所示。2016—2018年，不同地面沉降區(qū)的沉降空間演化存在顯著差異(如圖4所示)。其中位于唐山沉降區(qū)的TS剖面線處(圖4(e)和(f)所示)的沉降速率在監(jiān)測時(shí)段呈增長趨勢。而位于北京、廊坊—天津、保定—滄州—衡水、衡水—邢臺、邯鄲沉降區(qū)的BJ，LT，BCH，HX和HD剖面處2016—2018年沉降速率基本保持一致，沉降中心發(fā)育較為穩(wěn)定。綜上，InSAR監(jiān)測時(shí)段內(nèi)山前平原和中部平原地面沉降發(fā)展相對穩(wěn)定，濱海平原地面沉降范圍及沉降速率在監(jiān)測時(shí)段內(nèi)增長明顯，尤其是唐山沿海地區(qū)，沉降速率漲幅最為明顯。

(a) BJ西—東 (b) BJ北—南

(c) LT西—東 (d) LT北—南

(e) TS西—東 (f) TS北—南

(g) BCH西—東 (h) BCH 北—南

(i) HX西—東 (j) HX北—南

(k) HD西—東 (l) HD北—南

圖4 剖面時(shí)間序列InSAR形變速率結(jié)果Fig.4 Profiles of the InSAR displacement velocities

為了獲取2016—2018年研究區(qū)內(nèi)地面沉降時(shí)序特征，選取了8個(gè)不同位置上的沉降特征點(diǎn)，各沉降點(diǎn)的位置如圖3(b)所示，其對應(yīng)時(shí)序特征如圖5所示。北京朝陽金盞、天津王慶坨、保定大崔營、保定高陽、邢臺南宮和邯鄲臨漳等沉降區(qū)，形變特征以線性下沉為主，最大累積沉降量超過300 mm。 在唐山豐南、樂亭沉降區(qū)，區(qū)域時(shí)序沉降具有明顯的季節(jié)性特征，在2016—2018年沉降呈下降趨勢，如圖5(g)和(h)所示，但在2016年9月—2017年1月和2017年9月—2018年1月沉降減緩甚至出現(xiàn)回彈現(xiàn)象，之后地表繼續(xù)下沉，至2018年10月地面累計(jì)沉降量達(dá)180～254 mm。

(a) 北京朝陽 (b) 天津王慶坨 (c) 保定大崔營 (d) 保定高陽

(e) 邢臺南宮 (f) 邯鄲臨漳 (g) 唐山豐南 (h) 唐山樂亭圖5 研究區(qū)內(nèi)沉降特征點(diǎn)對應(yīng)時(shí)序特征Fig.5 Accumulative time-series deformation revealed by the InSAR technique

3.2 精度驗(yàn)證

為保證InSAR數(shù)據(jù)結(jié)果的準(zhǔn)確性，選用了26個(gè)2017年9月—2018年9月的水準(zhǔn)觀測值對Sentinel-1A數(shù)據(jù)所得垂向形變速率進(jìn)行精度驗(yàn)證。水準(zhǔn)數(shù)據(jù)位置如圖1所示，位于Track 142和Track 69影像覆蓋范圍內(nèi)，因此通過水準(zhǔn)測量數(shù)據(jù)對這2軌影像所得形變結(jié)果進(jìn)行驗(yàn)證。在驗(yàn)證過程中，以水準(zhǔn)點(diǎn)為圓心，100 m為半徑建立緩沖區(qū)，將緩沖區(qū)內(nèi)所有PS點(diǎn)的平均值視為該水準(zhǔn)點(diǎn)對應(yīng)的InSAR結(jié)果。驗(yàn)證結(jié)果如圖6和表2所示。

(a) Track 142 (b) Track 69圖6 InSAR結(jié)果與水準(zhǔn)測量結(jié)果精度驗(yàn)證Fig.6 Comparison between the InSAR measurements and levelling data

表2 InSAR所得形變信息與水準(zhǔn)測量結(jié)果比對Tab.2 Comparison of the mean subsidence rate between the Sentinel-1A PS and leveling data

其中Track 142所得研究區(qū)內(nèi)垂向形變結(jié)果與水準(zhǔn)測量結(jié)果相比，誤差最大值為11.6 mm/a，最小值為0.2 mm/a，均方根誤差為5.5 mm/a。對2組數(shù)據(jù)進(jìn)行線性回歸分析，R達(dá)0.97。Track 69所得研究區(qū)內(nèi)垂向形變結(jié)果與水準(zhǔn)測量結(jié)果相比，誤差最大值為13.5 mm/a，最小值為1.0 mm/a，均方根誤差為7.2 mm/a，R達(dá)0.98。驗(yàn)證結(jié)果表明InSAR監(jiān)測結(jié)果與水準(zhǔn)測量結(jié)果在同時(shí)間內(nèi)具有較好的吻合度，表明InSAR結(jié)果的可靠性。

為進(jìn)一步驗(yàn)證不同軌SAR影像結(jié)果的準(zhǔn)確性，對其進(jìn)行了交叉驗(yàn)證。由于Track 142的影像與其他2軌影像之間都具有重疊區(qū)，將其他2軌數(shù)據(jù)所得地表形變結(jié)果分別于Track 142所得結(jié)果進(jìn)行比對。結(jié)果如圖7所示，由Track 142與Track 69所得形變結(jié)果之間的相關(guān)性達(dá)0.96，由Track 142與Track 40所得形變結(jié)果之間的相關(guān)性為0.97。多軌SAR自相關(guān)驗(yàn)證結(jié)果表明了本次研究所用多軌InSAR數(shù)據(jù)結(jié)果的可靠性。

(a) Track 69與Track 142 (b) Track 40與Track 142圖7 多軌SAR影像形變速率交叉驗(yàn)證結(jié)果Fig.7 Consistency between the vertical displacement rates derived from different datasets

4 結(jié)論

本研究基于多軌SAR數(shù)據(jù)融合技術(shù)，獲取了2016—2018年京津冀平原區(qū)地表形變信息。不同軌SAR數(shù)據(jù)集交叉驗(yàn)證和水準(zhǔn)驗(yàn)證結(jié)果均表明本次研究InSAR測量所得地面沉降信息的可靠性，滿足地面沉降監(jiān)測及其演化規(guī)律的研究。

從沉降場演變過程來看，2016—2018年京津冀平原區(qū)內(nèi)沉降速率大于20 mm/a和50 mm/a的PS點(diǎn)數(shù)占總點(diǎn)數(shù)的比例均呈增長趨勢，沉降速率大于100 mm/a的點(diǎn)數(shù)占總點(diǎn)數(shù)的比例基本維持不變(3.6%，3.5%，3.5%)。InSAR監(jiān)測時(shí)段內(nèi)，濱海平原范圍內(nèi)地面沉降發(fā)展較為迅速，其中2016—2018年，該區(qū)域地面沉降嚴(yán)重區(qū)(年沉降量>50 mm)面積增長195 km2，漲幅達(dá)33.8%。沉降區(qū)時(shí)序剖面結(jié)果顯示，唐山沉降區(qū)呈加速發(fā)展趨勢，其余沉降區(qū)發(fā)展較為穩(wěn)定。InSAR所得地面沉降時(shí)序結(jié)果表明，研究時(shí)段內(nèi)，京津冀平原區(qū)內(nèi)地面沉降時(shí)序特征以線性下降為主，且唐山等沿海地區(qū)沉降區(qū)內(nèi)地面沉降時(shí)序具有明顯的季節(jié)性波動(dòng)特征。

研究成果一定程度上填補(bǔ)了京津冀平原區(qū)地面沉降整體演化特征研究較少的空缺，為京津冀地區(qū)地面沉降防治、調(diào)控，城市地下空間安全等提供科學(xué)依據(jù)。