周超凡， 宮輝力， 陳蓓蓓， 雷坤超， 施轢原， 趙 宇

(1.首都師范大學水資源安全北京實驗室，北京 100048； 2.首都師范大學地面沉降機理與防控教育部重點實驗室，北京 100048； 3.三維信息獲取與應用教育部重點實驗室，北京 100048； 4.首都師范大學資源環(huán)境與旅游學院，北京 100048； 5.北京市水文地質工程地質大隊，北京 100195 )

0 引言

地面沉降是一種由多種因素引發(fā)的緩變性地質災害現(xiàn)象，影響范圍廣、成因機制復雜、防治難度大。不均勻地面沉降發(fā)育嚴重，已對地下管道、線性交通等重要公共設施建設和運營安全產生較大影響[1]。

京津冀地區(qū)位于華北平原北部地區(qū)，區(qū)域地質構造復雜，水資源緊缺。為滿足城市建設和人口急劇增長的用水需求，長期超量開采地下水，使得區(qū)域內形成了大面積的地下水超采降落漏斗[2-3]。加之高密度城市集群快速發(fā)展、高速立體化交通網絡急劇拓展，導致區(qū)域內形成了多個地面沉降區(qū)，并呈現(xiàn)出跨區(qū)域連片分布的特征，沉降速率較大。區(qū)域內部分高鐵，如京滬、京津、津保和石濟高速鐵路已穿過沉降區(qū)。隨著區(qū)域地面沉降尤其是不均勻地面沉降的發(fā)展，會改變高速鐵路沿線的坡度，降低高速鐵路的平順性，影響高速鐵路的運營安全[4-5]。國內眾多學者對高速鐵路開展了監(jiān)測和預測等方面研究。在監(jiān)測方面，眾多專家早期使用水準、全球定位系統(tǒng)(global positioning system,GPS)、分層標等傳統(tǒng)技術對高速鐵路沿線沉降進行監(jiān)測，隨著合成孔徑雷達干涉測量(interferometric synthetic aperture Radar,InSAR)技術的發(fā)展，部分學者使用InSAR技術監(jiān)測高鐵沿線地面沉降[6-8]，研究發(fā)現(xiàn)，除了傳統(tǒng)的測量形變手段外，InSAR技術可用于線性地物的沉降信息監(jiān)測，監(jiān)測精度較高。在預測方面，主要有基于經典物理模型以及基于時序沉降演化特征的數(shù)學模型預測2大類。經典物理模型主要為基于太沙基固結理論和滲流理論建立的水-土耦合模型[9-10]，該類模型涉及參數(shù)較多，預測精度高; 另一類數(shù)學模型主要探查時序地面沉降演化過程特征[11-13]，在時間序列上建立地面沉降模型，對于物理參數(shù)限制較小，目前已從基于線性模型[14]發(fā)展到以機器學習為主的非線性模型[15]進行沉降量的預測。

津保高速鐵路是連接天津、河北及中西部地區(qū)的便捷通道，于2015年12月28日正式運營。高鐵沿線先后穿越王慶坨、廊坊和雄縣沉降區(qū)域，區(qū)域內最大沉降速率超過100 mm/a[16]，開展高鐵沿線地面沉降預測及坡度分析，對高鐵安全運營具有重要意義。本文利用小基線干涉測量-相干目標點分析(SBAS-IPTA)技術獲取京津冀平原典型區(qū)地面沉降時空分布信息，解決傳統(tǒng)隨機森林模型預測時未考慮時序數(shù)據(jù)內部復雜規(guī)律的問題，構建了小波變換-隨機森林(wavelet transform-random forest,WT-RF)預測模型，對小波變換分解后地面沉降趨勢分量和隨機分量分別學習獲取預測分量，再通過小波重構得到沉降完整預測信息并進行高鐵沿線坡度變化分析。

1 研究區(qū)概況

京津冀位于華北平原北部，平原區(qū)屬于華北平原的一部分，地理坐標為E114°16′～118° 04′，N36° 01′～40° 32′，主要包括北京平原、天津平原以及河北省京津之間和以南的平原區(qū)，總面積約為80 059 km2，占華北平原面積的26%，占全國陸地總面積的0.83%[17]。京津冀地區(qū)地質條件復雜、水資源緊缺，長期超量開采地下水，導致區(qū)域內形成了全球最大的地下水降落漏斗、范圍最廣的地面沉降區(qū)。

津保高速鐵路(圖1)在2010年3月21日開始建設，全長157.8 km，于2015年12月28日正式運營，最高時速不超過250 km/h。津保高速鐵路是連接天津、河北及中西部地區(qū)的便捷通道，在橫向上更是連通了京廣、京九、京滬鐵路三大繁忙干線，對加快京津冀一體化具有重要意義。然而，津保高速鐵路沿線穿過王慶坨、勝芳和雄縣沉降區(qū)[18]，開展高鐵沿線地面沉降坡度變化及預測研究對于高鐵運營安全具有重要意義。

圖1 津保高速鐵路位置及Sentinel-1A數(shù)據(jù)空間分布示意圖Fig.1 Location of the Tianjin-Baoding high-speed railway and the spatial distribution of Sentinel-1A data

2 數(shù)據(jù)源及研究方法

2.1 研究數(shù)據(jù)源

本研究數(shù)據(jù)包括140景Sentinel-1A(S1A)升軌存檔數(shù)據(jù)，時間跨度為2016年1月14日—2019年12月24日。數(shù)據(jù)覆蓋范圍和數(shù)據(jù)基本情況如圖1和表1所示。去除平地相位的外部數(shù)字高程模型(digital elevation model,DEM)數(shù)據(jù)選取美國國家航空航天局發(fā)布的SRTM DEM，空間分辨率為90 m。沿津保高鐵線間隔300 m取一個樣本點，利用樣本點提取小基線集干涉測量(SBAS-InSAR)獲取的沉降信息，采樣時間序列與影像獲取時間保持一致。將樣本點的沉降時間序列，批量重采樣為每個月的1號和15號，作為小波分解的輸入數(shù)據(jù)。

表1 雷達影像信息Tab.1 Radar image information

2.2 研究方法

本次研究的整體流程如圖2所示，首先利用SBAS-IPTA技術，獲取研究區(qū)2016—2019年地面沉降監(jiān)測信息(2019年形變結果僅作為預測模型的驗證數(shù)據(jù))，針對2016—2018年地面沉降監(jiān)測信息，在分析高鐵沿線坡度變化的基礎上，聯(lián)合小波變換和機器學習技術方法，精準預測高鐵沿線地面沉降。

圖2 整體技術路線Fig.2 Overall technology roadmap

1)SBAS-IPTA。SBAS-InSAR最早由Berardino等[19]在2002年提出，該方法為了提高干涉的相干性，將所有獲取到的雷達影像組合成若干個集合，保證在一定程度上解決了時間空間失相關等限制因素，最終獲取高精度時序地表形變信息。相干點目標分析技術(interferometric point target analysis，IPTA)是時序InSAR技術的另一種應用，相比于其他時序InSAR方法，IPTA技術通過矢量格式存儲數(shù)據(jù)結果，極大得減小了數(shù)據(jù)存儲容量并且提高了計算效率。另外，IPTA通過對相干目標點進行多次迭代處理，可獲取高精度地表形變信息。由于研究區(qū)覆蓋范圍和數(shù)據(jù)量較大，為保證數(shù)據(jù)結果精度的前提下提高處理效率，減小數(shù)據(jù)存儲量，故本研究結合SBAS-InSAR和IPTA方法，獲取京津冀典型區(qū)地面沉降信息[20]。

2)鐵路沿線坡度變化。本文選取高鐵坡度變化作為評價高鐵平穩(wěn)性的指標，《城際鐵路設計規(guī)范》[21]指出城際鐵路正線最大坡度不宜大于20‰，由于高速鐵路為百年工程，可以假定區(qū)域地面沉降對高鐵沿線坡度的年度改變安全范圍為小于0.2‰。在一定年限內，高鐵沿線2點產生的地面沉降坡度變化計算方法為[22]：

(1)

式中：i為高鐵沿線2點產生的地面沉降坡度變化； Δh為高鐵沿線2點地面沉降差值；n為沉降年限；a1和a2為高鐵2點上的沉降量；L為高鐵沿線2點間的距離。

3)小波變換。小波變換方法是Moelet等在1974年提出來的一種新的時頻分析方法。該方法通過對時間序列中的頻率信號進行細節(jié)分析，識別數(shù)據(jù)內部復雜規(guī)律信息[23-24]。時間序列數(shù)據(jù)是由多個序列通過交錯時間和疊加頻率組合而成，利用小波變換，分離出時間序列中的高頻信息和低頻信息，其中高頻信息為隨機分量、低頻信息為趨勢分量。設L2(R)為希爾伯特(Hillbert)空間內平方可積一維函數(shù)，若φ(t)∈L2(R)，其傅里葉變換ψ(ω)滿足條件(式(2)和式(3))，ω表示頻率，則φ(t)為母小波函數(shù)，公式為：

(2)

(3)

式中：a為尺度因子；u為平移因子。

小波分解和重構計算公式分別為：

(4)

(5)

式中：Wf(a,u)為小波分解；f(t)為時間序列數(shù)據(jù)。由于地面沉降受地下水開采、地層結構、動靜載荷等多種因素影響，時序地面沉降量為非平穩(wěn)的時間序列數(shù)據(jù)，但時序沉降仍然包含了趨勢分量，該趨勢分量與地下水的季節(jié)性開采有關，具有一定的周期性，而測量誤差等因素使得時序沉降中包含隨機分量。本次研究選取5級小波分解，利用上述公式可得到趨勢分量s5，隨機分量d1,d2,d3,d4,d5，基于隨機森林模型分別對趨勢分量和隨機分量進行預測。在此基礎上，利用小波重構將預測的趨勢分量和隨機分量結果進行疊加，得到最終地面沉降量，該方法可以在一定程度上減弱時間序列數(shù)據(jù)的擾動誤差。

4)隨機森林。隨機森林是Breiman在2001年提出的一種集成機器學習算法。隨機森林利用可重復取樣的方法對一系列的決策樹進行訓練，但同時根據(jù)各個決策樹的特點進行一定的改進，使得最終構建的決策樹間相關性盡量小，從而顯著提高了最終模型的性能[25-26]。隨機森林的具體算法為，首先對樣本數(shù)據(jù)進行有放回抽樣，獲得N個采樣集； 然后對每個采樣集隨機抽取m個特征構成N個決策樹； 接著對N個決策樹的輸出結果進行簡單平均法得到決策結果，完成訓練過程； 最后利用訓練得到的模型進行預測。模型中涉及到的參數(shù)分別有： 子模型數(shù)為決策樹的個數(shù)，不純度指標為決策樹做劃分時候的評價指標，葉節(jié)點最小樣本數(shù)為決策樹劃分時包含的樣本數(shù)量。地面沉降由趨勢和隨機等分量組成，對于傳統(tǒng)隨機森林直接對地面沉降數(shù)據(jù)進行模型訓練并未對沉降特性進行有效學習，本文中WT-RF將地面沉降趨勢分量和隨機分量作為多特征輸入給隨機森林進行學習，獲取預測分量，再通過小波重構得到沉降完整的預測信息。

本次研究的數(shù)據(jù)樣本時間序列為2016年1月—2018年11月，對每個月份取2個時間樣本，共有70個時間點。將2016年1月—2017年12月(48個)作為隨機森林模型的訓練數(shù)據(jù)集，將2018年1—11月(22個)作為模型的驗證數(shù)據(jù)集，用于評價模型性能。評價指標選取決定系數(shù)(R2)、均方根誤差(root mean square error,RMSE)和平均絕對誤差(mean absolute error,MAE)。R2越接近1、RMSE和MAE值越小，表明模型的預測效果越好。

3 研究結果與分析

3.1 津保高速鐵路時序地面沉降現(xiàn)狀

基于SBAS-IPTA方法，獲取2016—2018年京津冀平原典型區(qū)地面沉降分布特征。從圖3中可以看出，京津冀平原典型區(qū)地面沉降空間差異明顯，區(qū)域內形成了多個沉降漏斗，已構成復合沉降漏斗。主要分布在平原北部(北京市)、平原東北部(廊坊市和天津市)和平原南部(保定市、滄州市、衡水市和邢臺市)。津保高鐵穿過王慶坨—勝芳、雄縣和保定沉降漏斗。2016—2018年間，數(shù)據(jù)覆蓋范圍區(qū)域內最大沉降速率達到123 mm/a，有部分抬升值，分布在平原東北部山區(qū)。

圖3 2016—2018年京津冀平原典型區(qū)地面沉降空間分布特征Fig.3 Spatial distribution characteristics of land subsidence in typical areas of the Beijing-Tianjin-Hebei Plain from 2016 to 2018

進一步提取高速鐵路沿線2016—2018年平均沉降速率及高鐵沿線區(qū)段內沉降嚴重點時序沉降量(圖4)，進行鐵路沿線地面沉降時空差異特征研究。從圖4(a)中可以看出，2016—2018年，高鐵沿線地面沉降差異明顯，最大沉降速率為89 mm/a(P2)，發(fā)生在距離起始站91 km左右，位于雄縣沉降漏斗附近； 其次為P1點，位于王慶坨—勝芳沉降漏斗附近，距離起始站27 km左右，平均沉降速率達到46 mm/a； 距離終點站11 km附近的P3點，平均沉降速率達到44 mm/a。從圖4(b)中可以看出，2016—2018年，3個重點觀測區(qū)域地面沉降呈現(xiàn)減弱趨勢，P2點和P3點的沉降變幅分別達到45 mm和27 mm，雖然沉降處于減緩狀態(tài)，但仍需要重點關注。為進一步分析地面沉降對高鐵安全運營的影響，獲取2016—2018年高鐵沿線坡度變化信息。

(a) 鐵路沿線地面形變速率 (b) 3個重點觀測區(qū)圖4 2016—2018年津保高速鐵路沿線地面形變速率Fig.4 Rate of land subsidence along the Tianjin-Baoding high-speed railway from 2016 to 2018

3.2 津保高速鐵路沿線坡度演化特征

按照公式(1)，分別計算2016年度、2017年度、2018年度及2016—2018年度津保高鐵沿線沉降坡度變化情況(圖5)。圖5(d)結果表明，2016—2018年，津保高速鐵路沿線沉降坡度變化范圍為0～0.16‰。坡度變化小于0.05‰的區(qū)域范圍占總長度的89.41%，坡度變化在0.05‰～0.1‰，0.1‰～0.15‰和0.15‰～0.16‰的區(qū)域分別占比9.26%，1.13%和0.19%，線路長度分別為14.6 km，1.8 km和0.3 km，其中坡度變化大于0.1‰的區(qū)域長度約為2.1 km。進一步分析圖5(a)—(c)中2016—2018年各年度沉降坡度變化信息發(fā)現(xiàn)，2017年較2016年，坡度變化在0.04‰～0.06‰的區(qū)域呈現(xiàn)減小趨勢，2018年較2017年和2016年，坡度變化在0.02‰～0.04‰和0.04‰～0.06‰的區(qū)域均呈現(xiàn)減小趨勢。從以上分析結果可以發(fā)現(xiàn)，津保線路坡度變化范圍在工程技術標準之內，地面沉降暫時未威脅津保高鐵的運行安全。為進一步分析地面沉降對高鐵運行安全的影響，集成WT-RF模型，預測高鐵沿線地面沉降信息，獲取坡度變化情況。

(a) 2016年 (b) 2017年

(c) 2018年 (d) 2016—2018年圖5 2016—2018年津保高鐵沿線坡度變化Fig.5 Changes in slope along the Tianjin-Baoding high-speed rail from 2016 to 2018

3.3 津保高速鐵路沿線地面沉降預測

3.3.1 基于小波基的高鐵沿線地面沉降信息分解

為進一步評價地面沉降對高鐵運營的影響，需對地面沉降進行預測。選取小波分解方法，分解出現(xiàn)有地面沉降時間序列信息中的趨勢分量和隨機分量。隨機選取db5小波基對地面沉降信息進行分解，分解級數(shù)為5級，分解層數(shù)為6層[27]，由于表征地面沉降信息點數(shù)據(jù)量較大，本次僅展示P2點的小波分解結果，如圖6所示。

(a) d1分量 (b) d2分量 (c) d3分量

(d) d4分量 (e) d5分量 (f) s5分量圖6 地面沉降時間序列信息5級小波分解結果Fig.6 5-level wavelet decomposition results of land subsidence time series information

從圖6中可以看出，s5為地面沉降時間序列小波分解后得到的趨勢分量，d1—d5分別為每級分解得到的隨機分量，其中d1—d4隨機性明顯。在得到小波分解的結果后，選取隨機森林模型，將小波分解得到的趨勢分量和隨機分量分別作為模型的學習特征，獲得地面沉降趨勢分量和隨機分量的預測信息，再對預測結果進行小波重構，得到完整的地面沉降時間序列信息。

3.3.2 高鐵沿線地面沉降預測

將小波分解后的5個隨機分量和1個趨勢分量對應的時間序列沉降信息作為隨機森林模型預測的樣本集，其中70%(前48個時間點)的時序沉降量作為訓練數(shù)據(jù)，30%(后22個時間點)的時序沉降量作為驗證數(shù)據(jù)。隨機森林模型參數(shù)如表2所示。

表2 隨機森林模型參數(shù)Tab.2 List of random forest model parameters

利用測試集對模型進行驗證，模型的R2均達到0.8以上，表明模型性能準確。為進一步驗證聯(lián)合WT-RF模型方法預測沉降的效果，選取基礎模型ARIMA和隨機森林進行對比，對比結果如表3所示。

表3 3種模型精度對比Tab.3 Accuracy comparison of the three models

可以看出，基于WT-RF模型對地面沉降預測精度較高，其R2達到0.97，RMSE和MAE分別為5.87 mm和2.03 mm，在3種模型中為最小值。利用WT-RF模型對5個隨機分量和1個趨勢分量分別進行2 a(到2020年)的預測，并進行小波重構，得到地面沉降量的最終預測結果，如圖7所示。為進一步驗證預測結果的準確性，將預測得到的2019年地表形變信息與基于InSAR技術監(jiān)測的形變信息進行對比分析，統(tǒng)計津保沿線各研究點預測值與InSAR實測值差值情況，結果如表4所示?？梢园l(fā)現(xiàn)，基于WT-RF模型的預測結果誤差主要集中在1～5 mm范圍內，占總研究點數(shù)的36%，小于10 mm的研究點占總研究點數(shù)的76%。

圖7 基于WT-RF模型的津保高速鐵路沿線地面沉降預測結果Fig.7 Forecast results of land subsidence along the Tianjin-Baoding high-speed railway based on the WT-RF model

表4 基于WT-RF模型預測的2019年形變結果與InSAR監(jiān)測信息對比Tab.4 Comparison of subsidence results predicted by WT-RF model and InSAR monitoring in 2019

對預測的高鐵沿線地面沉降量，按照式(1)，分別計算2016—2019年和2016—2020年津保高速鐵路沿線沉降坡度變化情況(圖8)，結果表明，2016—2020年鐵路沿線沉降坡度變化范圍(圖8(b))較2016—2018年(圖5(d))均呈現(xiàn)增大趨勢，坡度變化范圍在0.05‰～0.1‰的線路增大了5.61個百分點，長度增加了9 km； 坡度變化范圍在0.1‰～0.15‰的線路增大了2.77個百分點，長度增加了4.5 km； 坡度變化范圍在0.15‰～0.2‰的線路增大了0.74個百分點，長度增加了1.2 km，使得坡度變化范圍在0.15‰～0.2‰的線路長度達到1.5 km； 綜上結果說明，區(qū)域地面沉降改變著津保高速鐵路坡度，為保證高速鐵路的運營安全，應采取措施控制地面沉降的發(fā)展。

(a) 2016—2019年 (b) 2016—2020年圖8 2016—2019年和2016—2020年津保高鐵沿線坡度變化Fig.8 Changes in gradients along the Tianjin-Baoding high-speed railway from 2016 to 2019 and from 2016 to 2020

4 結論與展望

本文采用SBAS-IPTA技術獲取京津冀平原典型區(qū)地面沉降時空分布，分析地面沉降對高速鐵路坡度變化的影響。構建WT-RF模型預測高速鐵路沿線地面沉降，評價地面沉降對高速鐵路坡度的影響。研究結論如下：

1)2016—2018年，京津冀平原典型區(qū)地面沉降差異分布明顯，區(qū)域內形成了多個沉降漏斗，并且構成了復合型沉降漏斗，研究區(qū)內最大沉降速率達到123 mm/a。

2)2016—2018年，地面沉降對高速鐵路沿線坡度影響的變化范圍在0～0.16‰，沿線坡度變化超過0.1‰的長度約為2.1 km。

3)基于WT-RF模型對地面沉降預測可以達到較好的效果。2016—2020年，地面沉降仍持續(xù)發(fā)展。高速鐵路沿線坡度影響的變化范圍在0～0.2‰，其中坡度變化范圍在0.15‰～0.2‰的線路長度達到1.5 km，應對地面沉降采取防控措施，保證高鐵的安全運營。

本文選取地面沉降時序結果作為機器學習模型的樣本集，側重使用時間序列地面沉降信息，基于時序地面沉降演化特征進行沉降的預測，而地面沉降受到多種因素影響，未來可將地面沉降的影響因素作為特征集，結合機器學習模型進一步預測時序地面沉降。

志謝：感謝歐洲航天局提供的Sentinel-1A系列數(shù)據(jù)。本論文中涉及的代碼已上傳至網址： https: //github.com/chaofan-cnu/paper-subsidence。