朱 武，竇 昊，殷那政，程意清，張雙成，張 勤

1. 長安大學地質工程與測繪學院，陜西 西安 710054； 2. 西部礦產資源與地質工程教育部重點實驗室 陜西 西安 710054； 3. 生態(tài)地質與災害防控自然資源部重點實驗室，陜西 西安 710054

膨脹土是一種高塑性黏性土，土體內含有較多的蒙脫石等黏土礦物，具有很強的親水性。同時，土體的含水量變化會引起體積變化，具體表現(xiàn)為吸水劇烈膨脹、失水劇烈收縮及反復脹縮形變。除脹縮性外，膨脹土還具有裂隙性和超固結性兩個重要特性，次生裂隙的產生與超固結作用都會破壞土體的穩(wěn)定性[1]。南水北調中線工程于2014年12月正式通水，從丹江口水庫調水，為沿線的鄭州、石家莊、北京、天津等14座城市提供生產和生活用水，是一項緩解我國華北平原水資源短缺、優(yōu)化水資源配置的重大戰(zhàn)略工程[2]。南水北調中線總干渠約有1/3段穿越膨脹土(巖)區(qū)域，干渠通水后其周邊地區(qū)土壤內的含水量必然升高，而潮濕態(tài)的膨脹土在經(jīng)歷降水和蒸發(fā)時其體積和強度都會發(fā)生變化[3]。特別地，如果在裂隙發(fā)育密集的范圍內有強降水滲入，裂隙周邊的巖土體也會吸水膨脹而逐漸向周邊擴散，使渠道表層巖土體達到飽和狀態(tài)，抗剪強度降低，在表層引發(fā)膨脹土淺層滑坡。因此，膨脹土的存在增加了南水北調工程沿線邊坡失穩(wěn)的可能性，為了避免膨脹土渠段滑坡對南水北調工程通水產生影響，對膨脹土區(qū)域進行形變監(jiān)測和特征分析是非常有必要的。

GNSS和應變計等地表觀測技術是目前膨脹土形變監(jiān)測的主要方法，但是其只能獲取離散的單點觀測結果，對于大面積未知區(qū)域的形變探測和監(jiān)測則存在局限性。合成孔徑雷達干涉測量技術(synthetic aperture radar interferometry，InSAR)作為雷達遙感的一個重要分支，具有低成本、大范圍、高精度等優(yōu)勢，能夠全天時全天候對地表進行觀測，已經(jīng)成為滑坡識別[4]、地面沉降[5]、DEM重建[6]以及地震變形監(jiān)測[7]的有效手段。目前，已有一些學者利用InSAR技術開展了膨脹土形變監(jiān)測方面的研究。文獻[8]基于TerraSAR-X數(shù)據(jù)獲取了荷蘭某農田短期內的膨脹土形變，并定性分析了黏土儲水量與地表位移的線性關系。文獻[9]采用PS-InSAR方法獲取了荷蘭某膨脹土堤壩長達6 a的時間序列形變，并結合外部數(shù)據(jù)與PS點的形變序列詳細分析了降水、氣溫及土壤類型對堤壩形變的影響。文獻[10]利用多時相InSAR技術獲取了法國某試驗基地的膨脹土垂向形變，并結合SMOS衛(wèi)星的土壤表面濕度數(shù)據(jù)估計了黏土層深度和厚度的變化。國內應用InSAR技術監(jiān)測膨脹土形變的研究相對較少。文獻[11]利用PS-InSAR技術獲取了南寧市區(qū)的地表形變信息，重點關注下伏基巖以第三系膨脹土為主的地鐵修建區(qū)，研究表明地鐵沿線沉降均在安全范圍內。文獻[3]結合南水北調中線工程輝縣段膨脹土在運營期間存在的問題，分析了PS-InSAR技術在該渠段膨脹土形變監(jiān)測中存在的難點和重點。文獻[12]應用多時相InSAR技術對南水北調中線地表形變進行研究，并結合降水和土壤濕度變化分析了輝縣市段膨脹巖土的形變機理。

盡管已有學者利用InSAR技術開展了膨脹土的形變監(jiān)測研究，多數(shù)針對施工過程中回填的平坦地表，較少開展膨脹土邊坡的研究，且未對脹縮形變進行有效分離?；诖?，本文收集了2017年3月12日—2021年4月8日的112景Sentinel-1升軌SAR數(shù)據(jù)，采用多時域InSAR技術獲取了南水北調中線某一膨脹土渠道邊坡的時間序列形變，并通過奇異譜分析(singular spectrum analysis,SSA)分解得到形變序列的周期項、趨勢項和噪聲項，開展膨脹土形變與區(qū)域內氣溫和降水之間的關系研究，解析膨脹土的形變特征，為南水北調中線干渠的邊坡穩(wěn)定性評價提供科學依據(jù)。

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

1.1 研究區(qū)域

研究區(qū)域位于南水北調中線工程某一膨脹土邊坡，該區(qū)域地處暖溫帶，屬大陸性季風氣候，四季分明，降水充沛，光照充足，夏季降水最多，冬季降水最少，每年最高氣溫出現(xiàn)在7月，最低氣溫出現(xiàn)在1月。該渠段地處黃河以南，相對南水北調中線其他渠段而言，區(qū)域內降水較多、地下水位較高，且地質條件復雜，沿渠存在大量膨脹土[13]。圖1展示了研究區(qū)域的地理位置、地形變化及膨脹土分布。圖1中標識的P1—P2段為高填方和高挖方膨脹土區(qū)域。圖1(b)為美國SRTM數(shù)字高程模型與德國TanDEM-X數(shù)字高程模型之間的差值，反映了該區(qū)域2001—2015年的地形變化?？梢钥吹?，南水北調渠道大部分為挖方區(qū)域，最大挖深約40 m。同時，部分穿越河流的渠道采取填充處理，最大填充高度約為20 m，如圖1(b)中的P2點附近。圖1(c)展示了渠道沿線P1—P2的地質剖面，可以看到整個沿線分布著厚度不等的膨脹土，其中大部分為強膨脹土，這樣的條件增加了南水北調工程沿線邊坡失穩(wěn)的可能性[14]。

圖1 研究區(qū)域DEM變化及P1—P2點地質剖面Fig.1 Data coverage, DEM changes and geological profile along P1—P2 over the study area

1.2 數(shù)據(jù)收集

為了分析研究區(qū)的時空形變信息，獲取了112景哨兵1號(Sentinel-1)升軌SAR數(shù)據(jù)。其成像模式為IW單視復數(shù)(single look complex，SLC)影像，成像方式為TOPS(terrain observation by progressive scans)，時間為2017年3月12日—2021年4月8日，具體參數(shù)見表1。同時，為了分析降水、氣溫與InSAR形變之間的關系，從河南省氣象數(shù)據(jù)探測中心收集了研究區(qū)附近的月平均降水量和月平均氣溫數(shù)據(jù)。數(shù)據(jù)采集時間為每日2、8、14、20時4次氣象觀測，然后計算月平均降水量與月平均氣溫。

表1 SAR影像參數(shù)

2 研究方法

對于獲取的112景SAR影像，主要開展差分干涉、時間序列和奇異譜分析3個步驟的處理，詳細流程如圖2所示。在差分干涉環(huán)節(jié)，首先進行多視處理，本文試驗按10(距離向)×2(方位向)的比例進行多視處理，從而抑制斑點噪聲，并使得距離向和方位向具有近似的空間分辨率[15]。其次，由于Sentinel-1 TOPS模式影像在采集時天線會發(fā)生旋轉，為了避免相鄰burst之間發(fā)生相位跳變，基于SAR影像相位和強度信息估計距離向和方位向的偏移量，隨后擬合偏移量多項式模型，并進行迭代處理，直至配準精度達到千分之一個像元[16]。然后，在選取適當?shù)臅r空基線生成干涉圖后，依次去除干涉圖中的地形相位、大氣相位、噪聲相位及趨勢項，將形變相位分離出來并完成相位解纏[17]。由于解纏誤差會影響時間序列形變的估計精度，在時序分析前需要剔除含有解纏誤差的干涉基線，為此試驗中選擇利用相位閉合環(huán)方法去除含有解纏誤差的干涉圖，選取最優(yōu)的干涉對網(wǎng)絡。時間序列形變采用短基線(SBAS)干涉測量方法進行求解。最后，為了揭示研究區(qū)域膨脹土的形變特征，開展了基于奇異譜分析(SSA)的時間序列形變分解研究，并分析形變與降水、氣溫之間的關系。

圖2 方法流程Fig.2 Flowchart of used method

2.1 相位閉合環(huán)選取干涉圖

在進行InSAR相位解纏處理時，由于去相干、地形復雜等原因引起的干涉圖相位不一致，增加了相位解纏的困難，不可避免地產生相位解纏誤差，從而降低了形變參數(shù)的解算精度。因此在進行時序分析之前，需要生成最優(yōu)的干涉網(wǎng)絡以提高結果的可靠性?；诟缮婢W(wǎng)絡的冗余性，本文采用了相位閉合環(huán)的方法識別并校正解纏誤差[18-20]。假設3景SAR影像φ1、φ2和φ3，以及與之相關的3幅解纏相位φ12、φ23和φ13，則相位環(huán)閉合差被定義為[16]

φ123=φ12+φ23-φ13

(1)

式中，φ12為SAR影像φ1和φ2之間的解纏相位；φ23為φ2和φ3之間的解纏相位；φ13為φ1和φ3之間的解纏相位；φ123為相位環(huán)閉合差。理論上，若3個干涉對φ12、φ23和φ13均不存在解纏誤差，則φ123等于0[21]；如果其中干涉圖含有解纏誤差時，則φ123為2π的整數(shù)倍。

本文基于相位閉合環(huán)對解纏后的干涉圖進行篩選。由于自動校正解纏誤差可能會產生更加難以識別的誤差，導致時間序列形變精度估計誤差，因此保守起見，逐干涉對而非逐像素評估環(huán)相位閉合度[22]，計算每幅干涉圖相位環(huán)閉合差的均方根誤差(root mean square,RMS)，將RMS大于給定閾值(1.5 rad)的干涉對從干涉圖網(wǎng)絡中剔除。圖3所示為2017年3月—2021年4月Sentinel-1數(shù)據(jù)干涉生成的時空基線圖，該時間段的數(shù)據(jù)冗余最多，更適合按照給定閾值自動進行篩選。

圖3 時空基線Fig.3 Spatio-temporal baseline network

2.2 時間序列形變估計

小基線集(SBAS)技術[23]是一種時序InSAR分析方法。該方法基于多個主影像，選取較短的時空基線進行干涉以減弱InSAR處理中的去相關問題，最后合并若干個干涉對子集求解形變時間序列及形變速率。目前，SBAS-InSAR方法已經(jīng)得到多方面的擴展應用[24-26]，可以滿足大部分研究區(qū)域與數(shù)據(jù)條件的需要。對于覆蓋研究區(qū)的N+1景SAR影像，基于短基線原則獲取M幅干涉圖。假設兩景SAR影像的采集時刻分別為ta和tb，并進行差分干涉生成第i幅干涉圖，則其任一點的差分干涉相位可以表示為[27]

(2)

式中，dta和dtb分別為ta和tb時刻相對于某一參考時刻的累計位移。本文試驗中采用三次模型[28]。如果以t0作為初始時刻，以除t0時刻外任意時刻ti(i=1,2,…,N)的干涉相位φ(ti)作為未知數(shù)，以經(jīng)過相位解纏解算的差分干涉相位φ(tj)(j=1,2,…,M)作為觀測值，則有

(3)

φ(ti)=[φ(t1)φ(t2)…φ(tN)]

(4)

φ(tj)=[φ(t1)φ(t2)…φ(tM)]

(5)

將式(2)—式(5)聯(lián)立可獲得M個方程組

Aφ=φ

(6)

式中，A為M×N階系數(shù)矩陣。由于SBAS方法在生成干涉圖時按一定的時空基線進行組合，有可能出現(xiàn)不連續(xù)的干涉圖子集，因此系數(shù)矩陣A可能秩虧，導致方程組(6)產生無窮多解。這時需要采用奇異值分解(SVD)方法，對多個短基線干涉圖子集聯(lián)合求解，獲得系數(shù)矩陣A的廣義逆，求出速度矢量在最小范數(shù)意義上的最小二乘解，進而求得形變時間序列。

2.3 基于SSA的形變分解

由于環(huán)境、系統(tǒng)等的影響，由時序InSAR技術所獲取的原始形變時間序列中不可避免地含有噪聲誤差，限制了形變信息的精度與可靠性。因此，將可靠的形變信號從原始形變時間序列中提取出來并確定其變化趨勢，對地表形變的時序分析具有重要意義。SSA方法建立在Karhumen-Loeve分解理論的基礎之上，是近年來飛速發(fā)展起來的一種研究非線性時間序列的方法。該方法無須先驗信息，不受正弦波假定約束，通過對一維時間序列數(shù)據(jù)進行主成分分析，構造出原始時間序列的軌跡矩陣并在此基礎上分解、動力重構，再結合經(jīng)驗正交函數(shù)，最終提取出趨勢項和周期項等有效信息[29-30]。一維時間序列x=(x1,x2,…,xS)的SSA過程可分為分解和重構兩個步驟。首先，選擇合適的分解窗口L(1

(7)

然后，采用奇異值分解算法(SVD)對軌跡矩陣XL×N進行分解得到

(8)

(9)

最后，根據(jù)式(10)對角平均化重構時間序列(y1,y2,…,yS)

(10)

基于時間序列分解原理，本文首先采用SSA算法得到了InSAR時間序列各分量特征值，隨后計算了時間序列各分量特征值的貢獻率?；谮厔菪巫兤交?、周期形變規(guī)律明顯、殘差符合白噪聲特點等原則，依據(jù)特征值累積貢獻率，將原始監(jiān)測數(shù)據(jù)劃分為趨勢項、周期項和噪聲項。

3 結果和分析

3.1 形變空間分布特征

利用收集的112景Sentinel-1A數(shù)據(jù)獲取了觀測時段的平均形變速率，如圖4(a)所示，其中綠色代表穩(wěn)定區(qū)域，黃色及紅色代表下沉區(qū)域，藍色代表上升區(qū)域。由圖4中可以看出，InSAR觀測點位主要位于附近村莊建筑物密集區(qū)域，以及南水北調兩側人工修筑的邊坡(圖4(b))，大部分農田則表現(xiàn)為失相干?？臻g分布上，除了東北角之外，南水北調干渠沿線大部分區(qū)域表現(xiàn)為抬升形變，最大年抬升速率約為18 mm/a。而干渠沿線東北角表現(xiàn)為下沉形變，最大年下沉速率約為15 mm/a。除此之外，干渠東側區(qū)域也表現(xiàn)為較為明顯的下沉，最大年下沉速率達到22 mm/a。進一步分析圖4(a)和圖1(b)，發(fā)現(xiàn)研究區(qū)的形變空間分布與地形變化相關，即挖方區(qū)域表現(xiàn)為抬升形變，填充區(qū)域表現(xiàn)為下沉形變。同時發(fā)現(xiàn)，形變量大小與挖填方深度相關，挖方越深抬升形變越大，填方越深下沉形變越大。

圖4 研究區(qū)地表形變速率和挖方渠道現(xiàn)場Fig.4 Surface deformation rate of study area and photo of excavation channel

3.2 形變時間分布特征

為了進一步調查研究區(qū)膨脹土的時間序列分布特征，繪制了圖4(a)中A、B、C和D4個點2017年3月—2021年4月近4年的時間序列形變，如圖5所示。A點位于干渠堤岸邊坡上，挖方深度為35 m，可以看到整體表現(xiàn)為明顯的抬升形變，觀測時段內累計抬升達到48 mm，時間上經(jīng)歷先緩慢下沉后加速抬升再緩慢抬升的變化特征。B點同樣為挖方區(qū)域，但是挖方深度只有25 m，相應的累計抬升量小于A點，觀測時段內約為20 mm，時間變化特征與A點相似，經(jīng)歷緩慢下沉后加速抬升再緩慢抬升。C點位于膨脹土填方渠道上，填方厚度為15 m，總體表現(xiàn)為持續(xù)的下沉形變，累計下沉量為63 mm，時間上具有一定的波動性，每年夏季表現(xiàn)為加速下沉特征。D點位于干渠東側填方區(qū)域，填方厚度為25 m，累計下沉量達到90 mm，時間波動性與C點相似。通過A、B、C和D點的時間序列分析，進一步證實了挖方抬升、填充下沉的形變特征，且形變量大小與填挖深度正相關，時間序列上表現(xiàn)為一定的波動性。同時發(fā)現(xiàn)在相同的挖填方深度情況下，填方下沉速度大于挖方抬升速度，如圖5中的B點。

圖5 研究區(qū)地表時間序列形變Fig.5 Time series deformation of the study area

3.3 時間序列形變分解和重構

膨脹土具有吸水劇烈膨脹、失水劇烈收縮及反復脹縮形變的特性，對其時間序列形變進行分解和重構，提取趨勢項、周期性和噪聲項，對分析膨脹土形變規(guī)律具有重要意義。為此，本文利用SSA方法，基于趨勢形變平滑、周期形變規(guī)律明顯、殘差符合白噪聲特點的原則，對圖5中的A點深挖方抬升形變和D點高填方下沉形變進行了分解和重構，如圖6所示。考慮到本文時間序列觀測個數(shù)為97及填挖方形變特征，設置A、D點奇異譜分析窗口長度l為12。圖6(a)為A點趨勢形變，其主要反映了上部膨脹土卸載情況下，下部土壤自然膨脹的過程，表現(xiàn)為抬升形變，占據(jù)圖5中A點總體形變量的75%。圖6(c)為A點周期形變，為降雨、氣溫、蒸發(fā)、地下水、荷載等外部因素引起的膨脹土形變，占據(jù)圖5中A點總體形變量的12%?？梢园l(fā)現(xiàn)，除了2017年之外，A點在2018—2021年表現(xiàn)為每年4至9月抬升形變，10月至來年3月下沉形變的周期性形變規(guī)律。這一變化與大氣降雨時間強相關，可以推測挖方膨脹土在每年4至9月集中降雨階段遇水膨脹，每年10月至來年3月少雨階段失水收縮。圖6(e)為A點噪聲形變，其主要由SAR成像幾何、數(shù)據(jù)處理等引起的誤差項，具有隨機性的特點，占據(jù)總體形變量的13%，具體分析中應去除此部分貢獻。

圖6 基于SSA的時間序列分解和重構Fig.6 Decomposition and reconstruction of time series deformation based on SSA method

與A點相比較，D點呈現(xiàn)出相反的趨勢形變(圖6(b))，主要為未壓實的填充膨脹土自然固結而產生的下沉形變，占據(jù)圖5中D點總體形變量的81%。圖6(d)中D點的周期項形變也主要受降雨、溫度等外在因素的影響，發(fā)現(xiàn)2017年和2019年具有較大波動，2018年和2020年波動相對較小，占據(jù)圖5中D點總體形變量的6%。圖6(f)為A點噪聲形變，占據(jù)總體形變量的13%。

4 形變成因分析

膨脹土邊坡渠坡形變是內外因綜合作用的結果。結合圖1(c)以及相關資料分析，研究區(qū)膨脹土主要成分為伊利石，其次是蒙脫石，此外還含有少量的高嶺石、綠泥石和水云母等[31-32]。該類膨脹土的滲透系數(shù)和自由膨脹率較大，具有較強的漲縮潛勢。因此，土體物理力學性質等內因對形變起主導作用。對于挖方膨脹土邊坡，由于膨脹土的超固結性和脹縮性，渠道開挖過程中會產生卸荷回彈，在垂直方向產生抬升形變，如圖5中的A、B點。對于填方膨脹土邊坡，由于膨脹土的超固結性，填方體的自然沉降固結導致沉降產生，如圖5中的C、D點。由于干渠兩側的膨脹土本身含水量較大，土體的膨脹性不是很強，因此抬升形變量較小。

膨脹土邊坡形變的外部原因主要是降雨、蒸發(fā)、溫度、地下水等引起的膨脹土土體含水率變化，進而導致膨脹土渠道脹縮。已有試驗結果表明，每次干濕循環(huán)(降雨與蒸發(fā)交替)，膨脹土邊坡均累積了向坡下的沉降和水平位移。其中，降雨決定土體濕化程度及入滲深度，是渠坡形變的最直接氣候因素。降雨入滲有一定滯后性，且由于地表徑流、水分蒸發(fā)等原因，入滲量要小于降雨量，因此由于降雨引起的形變具有一定的滯后性。圖7(a)為深挖方A點周期形變與月平均降雨量之間關系，可以看到隨著降雨量的逐漸增大，A點抬升形變逐漸增大，且形變會滯后2～3個月，之后隨著降雨量的減小形變開始下沉。氣溫變化會引起土體溫度變化，進一步引起含水率變化，從而導致膨脹土渠坡形變。圖7(b)為深挖方A點周期形變與月平均氣溫之間關系，可以看到溫度逐漸升高形變逐漸增大，溫度逐漸降低形變也逐漸降低，且也表現(xiàn)出一定的時滯性，表明溫度與形變具有較強的相關性。地下水是影響渠坡形變的另一重要因素，地下水不僅可加快結構面軟化，使得滑面抗滑力降低，還能在底滑面上提供揚壓力、在后緣拉裂面提供靜水壓力，導致渠坡形變和滑坡啟動。由于本文中的地下水位略高于渠道水位，因此渠道形變外部因素主要是降雨和溫度。

圖7 膨脹土邊坡形變與降雨和溫度之間的關系Fig.7 The relationship between deformation and precipitation and temperature

5 結 論

本文以南水北調中線工程某一膨脹土邊坡為例，采用多時域InSAR和SSA技術提取了膨脹土邊坡趨勢項和周期項形變，并結合膨脹土物理特性以及降雨、溫度數(shù)據(jù)，分析了膨脹土形變特征。得到如下初步結論：

(1) 挖方膨脹土邊坡表現(xiàn)為抬升形變，而填方邊坡表現(xiàn)為下沉形變。利用112景Sentinel-1A數(shù)據(jù)獲取的年平均形變速率顯示南水北調干渠沿線大部分挖方區(qū)域表現(xiàn)為抬升形變，最大年抬升速率約為18 mm/a，而干渠部分填充區(qū)域表現(xiàn)為下沉形變，最大年下沉速率約為15 mm/a。

(2) 膨脹土邊坡時間序列形變呈現(xiàn)一定的周期性變化，具有明顯的脹縮特性。通過SSA對InSAR獲取的時間序列形變進行分解，發(fā)現(xiàn)挖方和填方膨脹土邊坡周期項形變具有明顯的脹縮性，具體表現(xiàn)為挖方膨脹土每年4至9月出現(xiàn)抬升形變，10月至來年3月出現(xiàn)下沉形變。同時發(fā)現(xiàn)，膨脹土形變量大小與填挖深度正相關，挖方越深抬升形變越大，填方越深下沉形變越大。在相同的挖填方深度情況下，填方下沉速度大于挖方抬升速度。

(3) 卸載回彈和土體固結是膨脹土邊坡形變的內在因素，降雨和溫度則是形變的外部因素?？紤]到膨脹土的超固結性和脹縮性，渠道開挖過程中的卸荷回彈是挖方膨脹土形變的內因，而填方體的自然固結是填方膨脹土下沉形變的內因。通過分析膨脹土形變與降雨、溫度之間的關系，發(fā)現(xiàn)膨脹土周期性形變與降雨、溫度強相關，其造成的形變會時滯2～3個月。