基于PS
--InSAR技術的哈大高鐵扶余至德惠段地表形變監(jiān)測

王莉，王志明，楊瑩

1.西藏自治區(qū)重點公路建設項目管理中心, 拉薩 850001；2.交通運輸部天津水運工程科學研究所, 天津 300456

0 引言

隨著中國國民經濟的快速發(fā)展，道路交通網絡迅猛擴張，但是區(qū)域性的地表形變給道路交通的安全運行帶來了許多潛在危害。因此，對道路交通區(qū)域地表形變的監(jiān)測日益受到關注[1--2]。

傳統(tǒng)地表形變監(jiān)測主要通過精密水準測量、GPS測量等，采用人工方式展開，精度很高，但存在效率低、時效性差和覆蓋范圍有限等問題，很難實現大范圍的動態(tài)連續(xù)監(jiān)測。而合成孔徑雷達干涉測量(InSAR)技術，因具有全天時、范圍廣、受天氣影響小等優(yōu)點，可以獲取大范圍、高精度(mm級)的地表形變信息，在地表形變監(jiān)測方面得到了越來越廣泛的應用。Heleno et al.利用永久散射體干涉測量(PS--InSAR)技術對Lisbon地區(qū)的主要公路和鐵路線路進行沉降監(jiān)測，在道路沿線發(fā)現了兩個新形成的沉降區(qū)[3]。李英會基于穩(wěn)定點目標分析方法，利用PS--InSAR和小基線集干涉測量(SBAS--InSAR)方法開發(fā)了一種新的InSAR技術，對京津高鐵天津段2009—2010年的地面沉降進行監(jiān)測，發(fā)現區(qū)域最大沉降量達400 mm[4]。Perissin et al.使用InSAR技術監(jiān)測上海的地表形變情況，發(fā)現地鐵隧道和公路是地形沉降的主要發(fā)生區(qū)域[5]。楊艷等分析了京津城際高鐵線路運營五年多來的差異沉降與坡度變化特征，認為區(qū)域沉降導致的差異沉降量較大，但沉降變化仍符合高鐵軌道的設計要求[6]。張學東等融合PS--InSAR技術和SABS--InSAR技術，分析了京滬高速公路和上海市地鐵沿線的地表形變，得到了較好的結果[7]。劉歡歡等監(jiān)測了京津高鐵北京段地面沉降情況，并對其成因及影響因素進行了分析[1]。與傳統(tǒng)方法相比，PS--InSAR技術在道路區(qū)地表形變監(jiān)測方面，表現出精度較高、連續(xù)性好的優(yōu)勢。

但是，當前尚缺乏對中國東北高鐵運行區(qū)地表形變的監(jiān)測。哈爾濱到大連高速鐵路(哈大高鐵)是中國“四縱四橫”快速鐵路網京哈高鐵的重要組成部分，對東北地區(qū)的經濟社會發(fā)展具有重要戰(zhàn)略價值；也是國內開通的首條高寒地帶客運專線，運行路線經過松花江、拉林河、飲馬河及遼河等多條大河，惡劣的氣候、不穩(wěn)定的地下特征和復雜的地質結構都會帶來不可避免的地表形變，從而給高鐵的安全運行帶來潛在風險[8--9]。因此，筆者以跨越松花江、飲馬河的哈大高鐵扶余至德惠段為例，擬采用Sentinel--1B衛(wèi)星獲取的合成孔徑雷達(SAR)數據，嘗試基于PS--InSAR方法，進行哈大高鐵運行以來的地表形變分析。

1 研究區(qū)域和數據

1.1 研究區(qū)概況

哈大高鐵北起哈爾濱西站、南至大連北站，線路全長921 km。高鐵的運行極大縮短了中國東三省主要城市間的時空距離，為東北區(qū)域的經濟一體化創(chuàng)造條件，是中國東北地區(qū)的干線鐵路之一。哈大高鐵也是世界上首條在高寒地區(qū)運行的高鐵，路基凍脹等地表形變也是高鐵安全運行面臨的主要問題之一。

為了驗證基于PS--InSAR技術監(jiān)測高鐵行進區(qū)間地表形變的可行性，選擇哈大高鐵扶余至德惠段作為研究區(qū)(圖1)。研究區(qū)內高鐵全長約60 km，大部分位于中國東北季節(jié)性凍土區(qū)；高鐵運行區(qū)間經過了松花江及其重要支流飲馬河，高鐵沿線的土地利用類型以旱田、水田為主。耕地的季節(jié)性特征、城市的擴張，都增加了鐵路沿線地表形變監(jiān)測的難度。

圖1 研究區(qū)域示意圖(底圖為Sentinel--2B影像)Fig.1 Sketch map of study area

1.2 Sentinel--1B SAR數據

Sentinel--1B衛(wèi)星是歐空局(ESA)實施的涉及環(huán)境和安全信息服務的哥白尼計劃的一部分，搭載有C波段(3.8～7.5 cm)合成孔徑雷達(SAR)，主要用于陸地、海洋和海冰監(jiān)測服務。衛(wèi)星于2016年發(fā)射，軌道高度為693 km，重訪周期為12 d；中心入射角為38.89°，距離向分辨率為 5 m，方位向分辨率為20 m。合成孔徑雷達采用了4種工作模型，本研究選擇了寬幅干涉模式成像，VV極化方式。所有影像均為單視復數類型的數據產品。

研究中，選擇了2016年12月至2019年11月的32景SAR數據，覆蓋研究區(qū)近3年時間，采樣的平均時間間隔約為34 d。

另外，在PS--InSAR數據處理方面，衛(wèi)星星歷軌道信息在不同時相SAR的影像配準和形變信息生成等環(huán)節(jié)都起著重要作用。原始的衛(wèi)星軌道信息存在諸多誤差，會造成干涉圖中存在因基線誤差而出現的殘余條紋誤差。因此，基于精密的衛(wèi)星星歷軌道數據對衛(wèi)星原始軌道信息進行校正，可有效去除因軌道誤差引起的殘余條紋誤差[1,5]。為了確保Sentinel--1B數據的軌道精度，本實驗采用的是歐空局公布的POD精密定軌星歷數據，定位精度優(yōu)于5 cm。

從美國太空總署官方網站下載了由美國太空總署和國防部國家測繪局聯(lián)合測量的航天飛機雷達地形測繪使命(SRTM--DEM)數據，空間分辨率為30 m，為Sentinel--1B SAR數據的處理提供參考DEM，用于SAR影像集差分干涉處理時的地形相位去除[1,5]。

實地調查發(fā)現，鐵路附近90%以上的區(qū)域為耕地，平均氣溫、地表特征隨季節(jié)的變化很大；而PS--InSAR技術可以較好地抑制由于季節(jié)性溫度變化和植被覆蓋等因素引起的失相干現象。因此，本研究采用PS--InSAR技術進行研究區(qū)的地表形變分析。

2 數據處理

2.1 PS--InSAR方法

PS--InSAR技術是Ferretti et al.提出的，利用同一區(qū)域的多時相SAR影像，通過適當的數據處理和分析程序，將形變相位分量與其他相位分量分離[10]。PS--InSAR技術的關鍵是確定PS點，用以估計形變時間序列和形變速率。PS--InSAR方法的另一個結果是地形殘差，即給定PS點的散射相位中心的真實高度和該點的DEM高度之差，這也是獲得精確PS地理編碼的關鍵參數之一[10]。PS--InSAR技術處理流程包括以下幾個步驟：

SAR影像定標和配準進行SAR影像配準和輻射定標，將時間序列數據的幾何位置和輻射強度統(tǒng)一起來，用以消除成像過程中受衛(wèi)星視角、軌道高度、入射角、大氣效應以及SAR系統(tǒng)熱噪聲等造成的振幅信息變化。

生成差分干涉數據集選取其中一個時相的SAR數據作為主影像，其他時相數據作為副影像，分別進行干涉處理，生成干涉圖集；利用DEM數據去除地形相位，生成時間序列差分干涉圖集。

PS候選點的識別利用振幅離差指數法與相關系數法雙閾值法進行PS候選點的識別，其中，相關系數雙閾值法能夠快速縮小PS點探測范圍，振幅離差指數法能夠逐一進行PS點的精確篩選。

基于多圖像稀疏網格的相位解纏利用Delaunay不規(guī)則三角網建立PS候選點之間的連接關系，采用加權最小二乘法進行相位解纏。

大氣相位的估計和去除在估計每個PS點的線性畸變和DEM誤差后，從初始差分干涉相位減去殘余相位，包括非線性畸變相位和大氣相位，采用時域濾波和空間濾波進行分離。在估計時間序列SAR影像上PS候選點的大氣相位后，通過Kriging插值方法，得到SAR影像任意像素的大氣相位。

PS點的時序分析與形變量估計利用相干系數法，確定最終的PS點，重復第(4)步，重新構建Delaunay不規(guī)則三角網，進行多影像稀疏格網相位解纏，獲得所有PS點的時序形變值與形變速率。

2.2 Sentinel--1B SAR數據處理

本實驗使用Esri公司SARscape軟件進行PS--InSAR處理，選取高鐵沿線左右各延伸12 km作為研究區(qū)域。根據所選擇的32景Sentinel--1B SAR數據，時空基線的最優(yōu)值(圖2)，選擇2018年8月10日的SAR數據作為主影像，其余31景數據為從影像，分別與主影像進行配準和干涉處理。

經過對比分析，選取時成像時間為2018年8月10日的SAR數據作為主影像，共生成31個干涉對，其中最大時間基線600 d，最小時間基線為24 d(圖2)。

圖2 Sentinel--1B SAR數據的時間位置(a)和Sentinel--1B SAR干涉對的時間基線(b)Fig.2 Time position of Sentinel--1B SAR data (a) and time baseline of Sentinel--1B SAR interference pairs (b)

首先進行Sentinel--1B SAR影像集預處理，包括SAR數據導入、研究區(qū)裁剪和SAR影像配準、輻射定標處理。然后，根據所選擇的32景SAR數據，制作了空間基線和時間基線分布圖(圖2)；根據空間基線和時間基線的最優(yōu)值和多普勒質心差閾值，選擇2018年8月10日影像為主影像，由此將數據集劃分為不同的子集，生成31個干涉測量對；這樣，對研究區(qū)地表形變的平均采樣間隔約為34 d，最大時間基線600 d，最小時間基線為24 d，平均絕對基線和最大絕對基線分別為38.01 m和106.2 m。

利用SRTM DEM數據去除地形相位，利用POD精密軌道數據消除因軌道誤差導致的基線誤差和衛(wèi)星系統(tǒng)誤差導致的平地效應相位，生成31幅時間序列差分干涉圖集。將振幅離差指數法與相干系數雙閾值法相結合用以識別和確定候選PS點；基于Delaunay三角網法，確定相鄰PS點間最優(yōu)解關系，并進行空間上的相位解纏處理，利用已知的地面控制點解算任一點的線性形變速率與高程殘差；基于時空濾波處理，去除大氣相位誤差，進而計算各PS點整體相干性；進行二次解纏處理，解算最終PS點的稀疏格網，得到雷達視線(line of sight, LOS)方向年平均形變速率分布圖(圖3)。

3 結果與討論

3.1 研究區(qū)年平均地表形變特征

圖3是高鐵扶余至德惠段LOS向的年平均形變速率圖。圖3表明，從2016年12月到2019年11月，高鐵沿線的年平均形變速率為-16.8～9.1 mm/a，形變速率集中在-2.4～2.0 mm/a的變化區(qū)間，表明研究區(qū)總體的地表形變特征是比較穩(wěn)定的。

圖3 雷達LOS向的年平均地表形變速率分布圖(底圖為Sentinel--2B影像)Fig.3 Annual average surface deformation rate map in radar LOS direction

研究中共得到51 169個PS點，其中發(fā)生上升形變的PS點有19 626個，占總數的38.4%；發(fā)生沉降形變的PS點31 543個，占總數的61.6%；表明研究區(qū)的地表形變還是以沉降為主。發(fā)生上升形變的PS點集中在德惠市區(qū)的南部，結合實地調查和Google Earth影像，表明上升形變與城市區(qū)的擴張有關。沉降形變集中在扶余市區(qū)北部和南部、德惠市區(qū)西北部和松花江兩岸。其中，松花江附近的地表形變均為負值，最大的年平均沉降速率可達到-12.2 mm/a。而在扶余市區(qū)，地表形變特征表現復雜，西部和南部以沉降為主，最大沉降量達16 mm/a；而中部和東部則是上升形變和沉降形變同時出現。德惠市區(qū)以上升形變?yōu)橹?，僅在中部沉降形變以線性形式呈北西向展布。

在扶余至德惠段的高鐵沿線，地表形變特征表現特殊，以上升形變?yōu)橹鳌＝y(tǒng)計直方圖(圖4)表明，高鐵沿線的年平均形變速率集中在-1.8～2.0 mm/a，整體的年平均形變速率均值為0.44 mm/a。共采集PS點1 295個，發(fā)生上升形變的點883個，占總數的68.2%。發(fā)生沉降的PS點集中在松花江附近和扶余市區(qū)附近。

圖4 高鐵沿線PS點的年平均地表形變速率直方圖Fig.4 Annual average surface deformation histogram of PS points along high-speed railway

3.2 典型位置累積地表形變特征

為了更清楚地掌握高鐵沿線的地表形變特征，選擇4個典型PS點位置(圖3)，分別是扶余北站(A點)、松花江大橋(B點)、耕地(C點)和德惠西站(D點)，根據所生成的31個干涉影像對，制作了典型位置累積地表形變量的時序變化曲線(圖5)。圖5表明，A點和D點表現出了總體的沉降特征，B點的沉降量保持不變，而C點則出現了緩慢上升的地表形變特征。

圖5 典型PS點(圖3)地表累積形變量的時序變化曲線Fig.5 Variation of accumulated surface deformation with time in typical PS points as shown in Fig.3

但是，累積地表形變量隨時間的變化幅度非常大。A點在扶余北站附近，從2016年12月至2017年3月，累積地表形變呈明顯增加的趨勢；而在2017年3至4月，出現大幅度下降；從2017年4月至8月，累積地表沉降呈上升趨勢；之后，至2017年12月，累積地表形變呈下降趨勢，僅在11月份有小幅波動；從2017年12月至2018年3月，累積地表形變呈快速增大趨勢；至2018年5月，再次呈直線下降趨勢；隨后，至2018年9月，累積地表形變表現為緩慢上升趨勢，之后，至2018年12月，再次呈現快速沉降特征；2018年12月至2019年4月，總體表現為上升趨勢；隨后至2019年7月，再次表現為直線下降趨勢；2019年7月至11月，總體上表現為上升趨勢。

B點在松花江大橋附近，累積地表形變隨時間的變化表現復雜。從2016年12月至2017年4月，表現為增大的沉降形變；2017年4月至10月，累積地表形變波動很大，而2017年11月至2018年2月，表現為上升形變；2018年2月至11月，累積地表形變波動再次增大；2018年11月至2019年3月，再次表現為上升形變；之后，累積地表形變再次呈大幅波動形式。

C點在接近德惠的典型耕地區(qū)，2016年12月至2017年4月，累積地表形變量呈明顯增大趨勢；隨后，至2017年8月，呈沉降趨勢；2017年8月至2018年3月，累積地表形變一直呈緩慢增大趨勢，僅在2017年11月份有小幅波動；2018年3月至10月，地表形變波動很大；2018年10月至2018年12月、2019年2月至3月，累積地表形變均表現為增大特征，中間出現小幅減小特征；隨后，至2019年7月，呈現快速下降趨勢；7月至11月，累積地表形變特征表現復雜。

D點在德惠西站附近，2016年12月至2017年1月，地表形變有小幅增大；隨后，至2017年6月，呈緩慢下降趨勢；8月有小幅上升，之后再次沉降；至9月到達最大累積沉降量；2017年9月至2018年2月，累積地表形變再次表現為上升趨勢；之后，至2018年4月，再次呈現沉降特征；2018年4月至10月，累積地表形變量隨時間的變化波動很大；2018年10月至2019年2月，累積地表形變表現為快速上升趨勢；2019年2月至6月，表現為沉降特征；隨后至2019年11月，再次出現較大波動的累積地表形變量。

總體上，4個典型位置的累積地表形變都表現出了隨時間的明顯變化特征。

3.3 高鐵安全運行的潛在風險分析

綜合研究區(qū)總體的地表形變特征，認為高鐵扶余至德惠段安全運行的潛在風險有以下兩個：

(1)地表的季節(jié)性形變特征

對4個典型區(qū)的3年累積地表形變分析結果表明，無論是在高鐵站點(圖5，A點、D點)，還是河流(圖5，B點)、耕地區(qū)域(圖5，C點)，地表在每年12月至第二年2月(或3月)表現為上升形變，3月至6月表現為沉降形變，6月至10月的地表形變量變化很大。結合東北區(qū)的氣候特征，每年12月至第二年2月(或3月)的上升地表形變，明顯是寒冷氣候條件下由于土壤的凍脹生成的；而3月至6月的沉降形變應該與凍土層的融化有關；出現較大波動地表形變的6月到10月，正好處于東北地區(qū)的雨季。

因此，東北地區(qū)高鐵行進路線的地表形變，主要是季節(jié)性的路基凍脹和降雨造成的。

(2)高鐵線路附近的潛在風險

從年平均沉降速率圖(圖3)可以看出，在德惠西站和扶余北站附近都出現了大面積的地表沉降；從累積地表形變曲線(圖5)來看，這兩個點的時序變化曲線與其他兩個點的區(qū)別，在于每年8月或9月有一個明顯的上升形變。

結合東北區(qū)的氣候特征，初步分析很可能是地下水開采造成的地下漏斗，總體的沉降與沉降漏斗的加劇有關；而在8月或9月的上升形變，則是與該時期較豐富的地表降水有關，減少了地下水的開采量，并對地下水產生一定的補充。目前，尚缺乏扶余地區(qū)的地下水勘探資料。徐世明和李虎等分別在2005和2002年的地下水資源調查中，已明確指出，扶余地區(qū)地下水位呈連年下降趨勢[11--12]。另外，圖3表明，扶余市區(qū)的地表形變特征表現復雜，市區(qū)中部和東部的上升形變和沉降形變同時出現，進一步證實了該結論的可靠性；這里，上升形變很可能與城市東擴有關，而沉降形變則很可能來自新居住區(qū)的地下水開采。

總體地形沉降量在德惠市<4 mm，對高鐵的安全運行影響不大；扶余地區(qū)的累積地表形變量在2019年7月之前較為穩(wěn)定，之后變化幅度加劇，變化量接近或超過5 mm，這已經對高鐵的安全運行具有潛在的風險。

4 結論

(1)總體上，高鐵行進區(qū)間的地表形變表現穩(wěn)定，LOS向年平均形變速率在-2.4 ～2.0 mm/a之間；尤其是在松花江大橋附近，累積地表形變量隨時間的變化表現穩(wěn)定,表明研究區(qū)的地表形變特征符合高鐵安全運行的要求。

(2)季節(jié)性的路基凍脹和降雨是東北高鐵區(qū)的地表形變的主要影響因素。

(3)在德惠西站和扶余北站附近，都發(fā)現了可能與地下水開采有關的明顯地表沉降；其中扶余北站的地表沉降有加劇趨勢，已對高鐵的安全運行構成潛在風險。

由于缺少高鐵線路和地下水位監(jiān)測數據，該研究未進行地表形變監(jiān)測精度的評價。但從雷達LOS向的年平均地表形變速率分布圖(圖3)和高鐵沿線PS點的年平均地表形變速率直方圖(圖4)來看，研究區(qū)的地表形變速率在合理的變化范圍內，表明了該研究成果的可靠性。后續(xù)將結合精確的地表形變測量結果，進行哈大高鐵沿線的長時間序列地表形變監(jiān)測，為高鐵的安全運行提供科學的數據支撐。

致謝感謝美國地質調查局STRM DEM數據和歐空局Sentinel--1 B SAR數據的免費使用。