多時相Sentinel-1A InSAR的連鹽高鐵沉降監測分析

何秀鳳，高 壯，肖儒雅，羅海濱，馮 燦

1. 河海大學地球科學與工程學院,江蘇 南京 211100； 2. 南京信息工程大學遙感與測繪工程學院，江蘇 南京 210044； 3. 北方信息控制研究院集團有限公司，江蘇 南京 211153

近年來，隨著我國經濟發展和城市化進程的不斷加快，我國高鐵事業得到了快速發展。2013—2019年是我國鐵路投資最集中、強度最大的時期，截至2017年年末，我國高鐵運營里程為25 000 km，占世界高鐵里程總量的66%[1]。雖然高鐵為人口產業集聚、要素流動、時空壓縮和地方高鐵新城的涌現提供了一系列有利條件，但是在建造長距離軌道線路的過程中，脆弱的工程地質環境、縱橫交錯的地表構筑物載荷和地下水的過度開采等自然和人為因素，極易引發區域地表沉降和不均勻的縱向變形，給高鐵安全運營帶來了極大的隱患，甚至造成嚴重的經濟損失和不良社會影響。高鐵線路通常呈線性走向分布，具有距離長、跨度大、分布不均等特點[2]。在設計、施工和運營監測過程中，線路需要進行形變監測。而常規單一的形變測量手段，如水準測量、GNSS測量[3]和三維激光掃描，主要依賴監測站點的重復觀測，對觀測環境條件要求高，需要耗費大量的人力物力，難以監測高鐵沿線大跨度、長時間的沉降趨勢情況，從而給高鐵的運營和維護帶來極大的困難。

合成孔徑雷達干涉測量(interferometry synthetic aperture radar，InSAR)作為近些年發展迅猛的大地測量手段，具有監測范圍大、時空分辨率高，可利用數據源多，可以全天時全天候觀測的優點，在形變監測方面具有獨特的優勢。文獻[4—5]利用D-InSAR(differential InSAR)手段進行地表沉降監測，并融合GPS技術獲取三維形變場信息。但是常規D-InSAR方法易受時空去相關和大氣延遲誤差的影響，形變監測的精度和可靠性受到限制。為此，多時相InSAR(multiple temporal InSAR，MT-InSAR)技術在D-InSAR的基礎上逐漸發展起來，如永久散射體干涉測量(permanent scatterer interferometry，PSI)技術[6-7]和小基線集(small baseline subset，SBAS)技術[8-9]。它克服了常規D-InSAR中的時空去相關和大氣效應問題，在監測人工線狀地物形變上具有獨特的優勢，展現出極具潛力的應用前景。文獻[10]利用26景TerraSAR-X數據，通過PS-InSAR技術獲得了上海軌道交通的形變結果；文獻[11]利用相同的數據源和處理方法對上海寶山區1、3號地鐵線進行監測，揭示了地鐵沿線沉降的時空格局；文獻[12]利用天津地區某條鐵路沿線的水準數據對時序InSAR監測結果進行驗證，二者結果取得很好的一致性；文獻[13]利用相同的方法對京津線進行監測，獲取了京津線在運營過程中的沉降情況；文獻[14]對青海省內G214高速公路進行PS-InSAR監測，并探討了凍土季節性解凍現象對公路沉降的影響；文獻[15]利用X波段的COSMO-SkyMed數據，聯合PS-InSAR與DS-InSAR技術對北京至石家莊的高鐵展開監測，驗證了COSMO數據能夠獲取亞毫米級線性形變速率。

應用多時相InSAR技術監測人工線狀地物(如鐵路、高速公路)雖取得了一些研究成果，但大多都是宏觀的對監測結果進行解釋，加以定性的分析，對InSAR監測結果進行定量評價的較少，同時就應用的SAR衛星數據源來說大都是X波段的數據，應用C波段數據的案例鮮少。因此，本文以剛開始運營的連云港至鹽城的高鐵作為研究對象，獲取研究區2018年1月—2019年9月期間共47景C波段Sentinel-1A數據，利用MT-InSAR技術進行處理分析。此外，采用同時期內的BDS觀測數據進行對比驗證，探究應用C波段SAR數據的多時相InSAR技術在高鐵形變監測的應用能力。由結果可知，利用C波段的Sentinel-1A數據能夠獲取毫米級的線性形變速率和時序位移序列，與BDS監測結果相比取得很好的一致性，且連鹽高鐵線路整體表現穩定。但在與BDS監測站驗證過程中也發現，利用C波段SAR數據探測出的高相干性點在線路空間分布上具有不均勻性，在一些橫向形變明顯區域很難獲得高相干性點。

1 研究區與數據集

1.1 SAR數據集

圖1 研究區地理位置Fig.1 Geographic location of the studied HSR in Jiangsu province

連鹽高鐵位于我國蘇北平原，沿線分布眾多城鎮，人口分布密集，平均每平方千米600～700人。在過去的幾十年中，蘇北平原地下水面臨著過量開采的問題，地面沉降現象也因此變得嚴重起來。地下水開采和重大工程建設一直是引起該區域地表沉降的主要因素，特別是在過去15年，沿海城鎮化的加速與擴張加劇著部分地區地表沉降，對人工線狀構筑物帶來了極大的影響。而高鐵對軌道的穩定程度又要求極高，地表輕微的形變就會給高鐵運營帶來較高的潛在風險。為此，本文收集了高鐵建成后檢測階段共47景TOPS成像模式的S1數據，具體時間跨度為2018年1月10日至2019年9月2日。S1數據可以免費從歐空局下載，入射角為33.6°，影像覆蓋區域約為110×250 km2，像元分辨率在距離向是5 m，方位向是20 m。由于研究區在影像中相對較小，因而本文只對裁剪后的影像進行處理，裁剪區域長約30 km，寬約36 km，面積約為1080 km2。表1列出了影像的獲取日期以及時空基線分布情況。選取2018年8月26日獲取的影像作為主影像，其他從影像均配準并重采樣到主影像上。

表1 Sentinel-1A數據獲取時間以及干涉對參數

* 為參考主影像。

1.2 BDS實時形變監測系統

中國兵器工業集團下的北方信息控制研究院作為國內首家將北斗高精度衛星定位技術應用于監測鐵路基準點坐標準確性及路基沉降的單位，負責連鹽線北斗基礎測量基準控制網的布設。布設的控制網分為基礎控制網和加密控制網，其中基礎控制網由15個基準點組成，間距2 km，全長29.60 km；加密基準點布設在灌河特大橋的CPⅡ控制點處，共布設12個基準點，間距500 m，全長10.396 km，形成北斗鐵路加密測量基準控制網，為目標路段范圍內的變形監測及軌道控制提供測量基準。此外，在路基段四標段一工區選擇3個監測斷面，共布設9個路基沉降監測點，各斷面之間的距離約為80 m，進行基于北斗的鐵路路基沉降監測應用研究，并基于此研究確定路基沉降監測的方式及預警機制。各個BDS基準點、加密點和斷面位置如圖2所示。由于各個監測站與基準站距離相對較近，故BDS監測數據處理主要是采用靜態相對定位[16]。對于短基線的相對定位，大氣延遲(對流層延遲和電離層延遲)能較好地被削弱或消除。其中采用LAMBDA方法固定整周模糊度[17]，采用高次差法和Turbo Edit組合法進行周跳探測與修復。

圖2 連鹽線BDS形變監測系統裝置Fig.2 Photos of BDS deformation monitoring system alone Lian-Yan HSR

2 算法與結果分析

相比于城市地面沉降監測等應用領域，高鐵等這種線狀人工地物具有距離長、跨度大、分布不均等特點，給時序InSAR數據處理造成了一定困難，主要表現在：①高鐵在影像上分布不均，如走向上占很多像元，而鐵路的寬度可能只占幾個像元；②大部分線路位于非城市區域，成像范圍內后向散射強度低，經典PS選取方法得到的空間點分布較為稀疏，無法滿足沉降監測的要求，而簡單放寬選取條件可能引起大量錯選；③PSI算法無法像SBAS算法一樣可以使用相位環閉合差進行解纏相位誤差檢驗，當解纏相位含有明顯誤差時，一般最小二乘估計的結果會發生偏差。針對上述現實問題，本文對常規時序InSAR數據處理流程進行改進，通過應用適當的策略，具體包括采用全分辨率干涉、相位穩定性分析探測PS點和基于L1范數的抗差M估計解算時序形變等。

2.1 MT-InSAR技術及數據預處理

考慮到本文監測對象是線狀人工地物，影像時間間隔小且影像數量也足夠，故采用最具有代表性的PS-InSAR處理流程，并采用全分辨率法提取小尺度形變。PS-InSAR處理流程具體包括主影像的選擇、影像配準、差分干涉圖的生成、PS點的探測、PS點網絡構建、相位解纏、大氣效應去除、形變速率的估算。本文數據預處理利用POD(precise orbit ephemerides)精密定軌星歷數據去除平地相位，定位精度可達5 cm。利用30 m分辨率的SRTM DEM數據模擬地形相位。由于雷達是采用側視成像方式，波束斜向照射地表時會導致雷達圖像中的疊掩和陰影現象[18]，借助DEM(digital elevation model)的局部入射角影像和疊掩、陰影影像補償由傾斜地形引起的SAR后向散射。在將DEM配準、重采樣到SAR坐標系過程中，需要利用基線信息模擬干涉條紋。經驗顯示，即使是利用精密軌道數據計算得到的基線信息也不能保證完全正確，生成的差分干涉圖中會存在以斜坡形式、與地形或系統有關的殘余相位。因此，需要對生成的差分干涉圖進行快速傅里葉變換(FFT)確定殘余的基線分量，對基線模型進行改正，以達到最佳的干涉效果。

2.2 主圖像的選取

PS-InSAR算法首先需要根據總體相干性最大的原則選擇主影像，然后將其他從影像與之配準，生成差分干涉圖集。公共主影像的選取需要顧及時空基線和多普勒質心頻率差，其函數模型[19]為

（3）完善國家創新體制機制，提高政府管理效率。從美國創新戰略可以看出，美國政府將民眾、學生和企業加入創新體制中，共同關注世界熱點問題，解決現實難題。在全球科技創新主體多元化、創新形式多樣化的大背景下，為滿足我國科技創新發展需求，提升科技成果轉化效率，我國的科技創新體制機制也需要進一步的完善和更新。政府在科技創新中的地位和作用有待進一步明確，在堅持以市場為導向、企業為主體、政策為引導、產學研協同創新的體制機制下，確立企業在產業導向的科技計劃中決策者、組織者、投資者的地位。與此同時，無論是從法律、法規層面的約束，還是普惠性政策的實施，都要確保實施過程得到了有效的監督、評估，提高政府的管理效率。

(1)

式中，γtotal為綜合選取系數；T為時間基線；B為垂直基線；FDC為多普勒質心頻率差；C表示各自的臨界參數值。根據式(1)計算每一幅影像的選取系數，系數最大的作為主影像。最大綜合選取系數的影像為2018年8月26日獲取的影像，將其作為主影像，共生成46個干涉對，平均相干性為0.82。連鹽高鐵區域生成的時空基線圖如圖3所示。

圖3 46對差分干涉對時空基線分布Fig.3 Baseline-time plots for the 46 differential interferograms

2.3 高相干點目標的選取

像素子視相干系數閾值法、相位離差法和振幅離差指數法是進行PS點探測最常用的幾種方法[20]。這幾種方法均是利用高相干點的穩定相位特征或者低相位離散特性，但是各有優缺點，實際數據處理中需要結合其中兩種或多種方法進行高相干點的選取。由于連鹽鐵路沿線大部分位于非城市區域，而振幅離差方法適合城區高相干點的探測，故采用常用的振幅離差方法探測出的PS點密度可能不高，不利于分析線路的沉降細節特征。為了增加沿線區域PS點的密度，本文采用相位穩定性分析方法。該方法判定點在時空范圍內的相位特征和基于點的空間距離越小相位相關性越高的假設。因此，首先進行空間低通濾波，然后計算相位離散度，進而選取具有穩定相位的點。這種選取方案，既保證了點位密度，也為三維相位解纏提供足夠多可靠的點。首先利用振幅離差指數閾值(閾值設為0.4)方法，選擇在長時間間隔內保持穩定的候選點集(PSCs)；然后運用相位穩定性分析估計每個候選點的噪聲相位，將PS點概率函數值作為權重迭代進行精化。基于計算出每個候選點的時序相位相干因子γx值并結合振幅離差選取PS像元。γx計算公式為

(2)

圖4 研究區PS點分布Fig.4 Distribution maps of PS in the study area

2.4 沉降速率和沉降量分析

對于監測高鐵這種線狀人工地物，C波段S1數據的分辨率略低于TerraSAR-X和COSMO等分辨率達1～3 m的衛星影像數據。為了盡可能地提高可監測形變梯度，文中采用全分辨率法提取小尺度形變。所謂的全分辨率干涉就是不對影像進行多視操作，直接使用SLC數據進行干涉處理。這樣既可以保證對強相干目標如鐵路、橋等構筑物準確識別，也提高像元中強散射體占主導地位的概率，獲得更多的高相干點目標，不會造成孤立強散射體。

基于選取的高相干點目標，采用3D相位解纏方法[22]：①在時間域內，基于這些PS點建立Delaunay三角網，建立相鄰兩點之間的相位模型，經過去除大氣效應和殘余的軌道誤差后，對時間序列相位差進行低通濾波和相位解纏；②運用解纏的差分相位時間序列，對每一個解纏相位值構建一個費用函數，最終得到解纏后的相位如下

Δφx,i=W{Δφx,i}+2kx,iπ=Δφh,x,i+φdef,x,i+φorb,x,i+φatm,x,i+φnoise,x,i

(3)

式中，Δφx,i為第i對干涉圖像元x的解纏相位；kx,i為相位整周數，若解纏精確，則kx,i對干涉圖上大部分像元是相同的；W{·}代表纏繞；Δφh,x,i為DEM不精確導致的地形殘差相位；φdef,x,i為視線向上形變相位；φatm,x,i為傳播過程中大氣引起的路徑延遲相位；φorb,x,i為殘余軌道誤差相位；φnoise,x,i為熱噪聲、配準誤差等引起的噪聲相位。

為了獲取最終的形變相位，還需要扣除各噪聲分量。對于軌道誤差相位，文中采用二階多項式對方位向和距離向進行擬合去除。對于殘余高程和形變相位的估計，一般采用最小二乘(LS)估計。然而LS估計沒有考慮到殘余軌道誤差和大氣相位的影響，它只是一個理想的全局最優解。只有當觀測相位值是相互獨立并且服從正態分布，解算的結果才是無偏和有效的，當解纏相位觀測值含有明顯殘余誤差時，LS估計的結果會發生明顯歪曲。因此，本文在解算參數時，采用抗差M估計[23]為

(4)

由于試驗區地形平坦，大氣中垂直分層效應不明顯，湍流效應占據主導地位。而湍流大氣一般被認為在時間和空間上是隨機分布的，為了避免引入大氣誤差，數據處理中并沒有直接利用外部數據進行大氣相位的去除，而是采用時空濾波進行大氣相位估計，最終得到平均沉降速率圖和監測周期內時序位移序列。圖5為分別利用抗差M估計和最小二乘方法解算的形變速率標準差分布情況，可以看出利用抗差M估計解算的形變速率標準差比LS估計的結果普遍要小，形變速率標準差降低的PS點約占89.2%。圖6為利用抗差M估計得到的線性沉降速率結果。

圖5 形變速率標準差分布圖Fig.5 Distribution map of the standard deviation of deformation rate

圖6顯示了利用時序InSAR獲取的整個區域在2018年1月至2019年9月期間的平均線性形變速率，整個區域最大形變速率達到了15 mm/a，而高鐵沿線區域形變更是微小，平均形變速率約為5 mm/a。在將近兩年的監測時間內，人工構筑物保持了很好的相干性。就高鐵沿線區域而言(圖6(b)所示)，盡管相較X波段，C波段的S1數據地面分辨率中等，但本文采用的是全分辨率干涉，保證了最高分辨率，共獲取10 199個PS點，密度遠遠超過BDS監測點。結果表明利用C波段S1數據監測像高鐵這類人造線狀地物目標，具有很大的應用潛力。

圖6 線性形變速率Fig.6 Linear deformation velocity map

3 沉降監測結果對比驗證及分析

MT-InSAR技術監測的是視線向形變，而BDS監測的是垂直方向和水平方向形變，為了驗證MT-InSAR的結果，一般假設沒有水平變形，將InSAR視線向形變轉換到垂直方向上。根據雷達成像幾何，將InSAR監測結果除以cosθ(θ為入射角，每個PS點入射角不同)即可完成轉換。BDS監測數據后處理采用短基線相對定位，精度在亞太地區優于GPS，在垂直向優于2 cm，水平向上優于1 cm[24-25]。關于InSAR結果的驗證，本文將從兩個方面展開，即依據BDS監測結果驗證比較InSAR的平均形變速率和時序位移量結果。

3.1 平均形變速率結果驗證

根據BDS測站的位置選取相應PS點，將計算得到的BDS形變速率與PS結果進行比較。值得注意的是，BDS監測時間段與InSAR數據時間范圍并不完全一致，僅是其一部分，因此本文僅比較兩者在重疊時間內的形變速率，并引入均方根誤差(RMSE)評價精度。RMSE計算公式為

(5)

式中，Δvk為第k個BDS監測站形變速率與InSAR結果形變速率之差；n為測站點的個數。

其比較結果見表2和圖7所示。

表2 BDS監測點與PS點形變速率結果比較

由表2可以看出，BDS監測點與PS點形變速率均方根誤差為6.2 mm/a，最大偏差為13.9 mm/a，且偏差落在5 mm/a內占70.6%。偏差大于1 cm/a僅有2個測站(圖7(b))，對應圖7(a)中的序號2和17，回歸分析結果負相關系數R為0.79，表現出很好的一致性。關于圖7(a)中序號2和17對應的測站出現較大偏差的原因，4.2節將詳細分析。圖8顯示了MT-InSAR在時間跨度為21個月內，高鐵沿線各個BDS監測點位置的平均形變速率變化情況，從圖8中可以發現整條線路平均形變速率在1 cm/a內，線路比較穩定。為了對比驗證MT-InSAR監測的精度，又分別獲取了二者在重疊時間內(2018年11月～2019年9月)的沿線平均形變速率變化情況，可以發現InSAR結果與BDS監測的結果較為吻合，趨勢相近。

圖7 BDS監測與MT-InSAR監測結果回歸分析與誤差直方圖Fig.7 The regression analysis and error histogram of velocity and displacement comparison between PS results and BDS measurements along a railway

圖8 高鐵沿線各個BDS監測站與PS點平均形變速率結果Fig.8 Average subsidence velocity comparison of BDS monitoring stations along the railway

3.2 時序位移量結果對比驗證

以BDS監測站天線位置為中心，100 m作為半徑搜選所有PS點，并求取該范圍內所有PS點的平均時序位移序列。為了匹配兩種技術獲取的位移序列，需要將位移時間序列轉化為重疊時間內共同的時間參考。由于篇幅所限，文中只顯示了其中10個監測站與MT-InSAR獲取的時序位移量的比較結果。

為了比較MT-InSAR監測結果，本文將2018年11月6日設為時間參考基準點，求取BDS與MT-InSAR相對基準點在垂直方向上的位移變化量。從圖9中可以看出，除了JM0006Q和JC0006Q測站，其他測站的結果十分吻合。JM0006Q和JC0006Q測站偏差相對較大是因其周圍的PS點很少，且距離周圍最近的PS點相對較遠。選取的PS點并不能完全反映這兩個測站的真實形變情況，而地表形變具有空間相關性，故雖然整體略有偏差，但整體形變趨勢較為相近。從圖9中還可以發現，2018年11月之前，形變位移序列起伏較大，例如JC0007Q、JC0015Q站，除了有一部分原因是殘余大氣的影響，還可能是因為前期施工過程中地基需要一個穩定的過程，即下沉-抬升過程。

為了揭示前期沉降過程，利用同樣的算法，但只處理了前期即2018年間獲取的影像，結果如圖10所示。對比圖6形變速率，可以發現前期地面沉降更為明顯。尤其是第2段線路(圖10(c))，形變速率最大達到2 cm/a。從早期整個線路的形變來看，可以看出前期施工后一段時間內線路出現地基回升與緩慢下沉，這在鐵路建造過程中是常出現的現象。再結合圖9，可以發現線路在進入2019年后，地基逐漸變得穩定，無論抬升還是下沉均在數毫米范圍內波動。

圖9 InSAR與BDS在垂直方向上時序位移對比結果Fig.9 The displacement comparison along the vertical direction between InSAR and BDS results along Lian-Yan HSR

表3為各個BDS監測站與其相應的PS時序位移統計結果，其中最大均方根誤差為6.1 mm，最小為1.6 mm，平均均方根誤差為3.8 mm。選取的所有樣本點中最大偏差為14.1 mm，偏差位于2 mm內的樣本占樣本總數的47%，最小偏差為0，二者偏差大于5 mm僅占19%，以上可以看出二者結果是十分吻合的。

表3 BDS監測點與PS點均方根誤差

3.3 連鹽高鐵沿線地面沉降分析

為了更加全面的分析連鹽線的形變特征，圖11可視化了3個區域(A、B和C區域)中的一些探測點的形變序列。在A和C區域中，探測到的PS點主要分布在線路上以及相鄰的人工建筑物。其中A區域位于第1段線路(圖10)，選取的兩個PS點，形變速率分別為10.7 mm/a和12.3 mm/a。從其形變序列上可以看出，二者有著相近的形變趨勢，形變主要發生在2018年5月至8月，具體表現為抬升，幅度大約在1 cm左右，8月之后，地基逐漸變得穩定。C區域主要位于第2段線路灌河特大橋北側，可以發現區域內的形變特征雖然和A區域相近，但是在2018年間形變更為劇烈，這也驗證了圖10得到的第2段線路較第1段線路形變速率更為明顯的結論。推測這可能是由于該區間建成時間比第1段線路晚，且位于軟土分布區，主要表現為高含水量、高壓縮性、低承載力，因此地表固結穩定的時間相對延后。相比A區域，C區域在2018年5月至2018年8月地基回彈更為劇烈，時間更長，直到12月中旬線路開通后才逐漸穩定。無論A區域還是C區域在線路開通運營后，線路均變得較為穩定。而B區域是連鹽線路基段四標段一工區，區域內有3個監測斷面，共布設9個路基沉降監測點，各斷面之間的距離約為80 m。可以發現在該區域并沒有探測到相干點，顯現了高相干性點在線路空間分布上具有不均勻性的特征。但是從9個監測站的監測結果來看，3個斷面在2018年11月至次年4月的監測周期內，沒有明顯的形變發生，只在1～2 mm之間變化。以上可以得出結論，連鹽線總體表現穩定，尤其是12月中旬線路開通運營后。

圖11 連鹽線重點監測區域時序形變Fig.11 Time series deformation of the points in the areas A, B, C along the Lian-Yan Line

4 結 語

本文將C波段Sentinel-1A數據應用到連鹽高鐵沉降監測，采用多時相InSAR技術進行處理并與同時期獲取的BDS監測結果進行對比驗證分析。結果發現，在連鹽高鐵建造完成至試運營階段，平均形變速率在1 cm/a之內，整體線路路基在2018年11月之前經歷一段時間的沉降與回彈后，11月以后，尤其在進入2019年后，沿線逐漸平穩。與BDS監測數據驗證結果來看，二者線性平均速率均方根誤差為6.2 mm/a，最大偏差為13.9 mm/a，最小偏差為1.1 mm/a，相關系數達到了0.79；10個BDS監測站平均時序位移均方根誤差為3.8 mm，且與BDS監測結果呈現微小的系統誤差影響，偏差在2 mm內約占47%。無論是獲取的平均形變速率還是時序形變序列，均取得很好的一致性，顯現出利用C波段S1數據對人工線狀地物進行時序InSAR監測的巨大潛力。同時在監測過程中，雖然利用多時相InSAR技術獲得了相當高密度的PS點，但是試驗發現在3個重點監測的斷面處，并沒有獲取足夠的PS點，導致無法與布設的BDS監測點進行比較，這也暴露了在利用C波段SAR數據進行時序InSAR監測過程中，雖可以獲取可觀密度的PS點，但是在空間分布上具有不對稱、不均勻的特點。因此，綜合考慮高鐵線路的走向、坡度以及雷達成像視角，探測沿線區域分布更加均勻、穩定可靠的相干性點，并獲取高鐵線路的三維形變場，有待進一步深入研究。

致謝：感謝歐空局提供的2018年1月至2019年9月的Sentinel-1A數據。