多年凍土區地表變形與影響因素相關性分析

趙 韜，張明義，路建國，晏忠瑞

(1.中國科學院 西北生態環境資源研究院 凍土工程國家重點實驗室，蘭州 730000；2.中國科學院大學，北京 100049)

中國多年凍土分布極為廣泛，其面積約為2.15×106km2，占中國陸地總面積的22.4%[1]。多年凍土凍融過程中冰水相變導致的土體“體縮”和“體脹”現象會使地表產生變形，無論是體積收縮引起的沉降變形還是體積膨脹引起的抬升變形，均會對多年凍土區生態環境和基礎設施的穩定性造成影響[2]。研究多年凍土區地表變形與影響因素間的相關關系，在豐富凍土學理論的同時可以為多年凍土區地表變形問題的深入研究提供參考。

多年凍土區地表變形有諸多監測方法，目前常用的有水準測量[3]、變形儀器測量[4-6]、GPS測量[7]等。上述方法可以獲得多年凍土區單點高精度的地表變形信息，但在大時空尺度變形測量中存在很多局限性。近年來，隨著遙感技術的快速發展，利用合成孔徑雷達干涉技術(interferometric synthetic aperture radar，D-InSAR)對多年凍土區地表變形的監測取得了很多研究成果。有研究表明，InSAR監測結果與野外現場實測存在一定差異，但總體而言，InSAR技術為多年凍土區地表變形大時空尺度監測提供了新的方法，其監測精度可達厘米至毫米級[8-12]。

多年凍土區地表變形是諸多因素綜合作用的復雜過程。國內外學者從不同的角度分析了地表變形的影響因素，并探討了地表變形與影響因素間的關系。張建明等[5]從多年凍土區地表變形機理來源出發，指出多年凍土區路基工程地表總變形源于路堤的壓密變形、活動層的凍融循環變形及凍土層的融沉變形，主要受活動層厚度、年平均地溫、體積含冰量、地質構造等因素的影響，且體積含冰量越大，地表融沉變形越大。馬巍等[6]結合大量的野外實測數據，對比分析得出多年凍土區地表變形演化過程與下伏多年凍土溫度變化一致。董昶宏等[13]通過對比分析66個監測點野外變形數據與年平均地溫、體積含冰量、氣溫、太陽輻射、地質條件等因素之間的關系，指出多年凍土區地表變形與年平均地溫、體積含冰量、工程地質條件密切相關，且地表沉降量與體積含冰量呈正比關系。受地表變形實測數據來源限制，上述對多年凍土區地表變形影響因素及相關性的研究主要集中在單點尺度上，且主要以定性分析為主。近年來，隨著遙感技術的快速發展，一些學者開展了面域范圍內多年凍土區地表變形與影響因素相關性的定量分析。Zhao等[10]分析得出青藏高原多年凍土區地表變形與當地氣溫和降水有明顯的負相關關系，相關系數分別為-0.80～-0.45和-0.95～-0.75。曾旭倩等[14]研究表明，東北多年凍土區地表變形隨土壤含水量的增大而增大，且二者存在明顯的正相關關系，相關系數為0.51。但總體而言，目前在面域范圍內對多年凍土區地表變形與影響因素相關性的研究還很少。在全球氣候持續升溫和人類活動日益加劇的大背景下，多年凍土區地表變形在諸多因素的影響下將如何發展變化，是亟待探究的復雜問題。

鑒于此，首先采用InSAR技術獲取了2015年6月—2019年6月間青藏工程走廊多年凍土區的地表變形信息，并利用野外實測數據對其精確度進行了驗證；然后，借助GIS平臺，根據經驗計算模型，獲取了研究區地表變形影響因素基礎數據分布情況；最后，利用簡單相關和偏相關分析法，探究了研究區地表變形與影響因素間的相關性，結果可為多年凍土區地表變形問題的深入研究提供科學參考。

1 研究區概況

青藏工程走廊始于青海省格爾木市，止于西藏自治區拉薩市，全長約1 120 km，穿越多年凍土區約550 km，其中高溫凍土約275 km，高含冰量凍土約221 km，高溫高含冰量多年凍土約134 km[15]。本文以走廊內多年凍土區段(西大灘—安多)為研究區，探究多年凍土區地表變形與影響因素的相關性。如圖1所示，研究區地理坐標位于東經91°E～95°E、北緯32°N～36°N之間，海拔介于3～7 km。研究區地形地貌復雜多變，包括中高山區、高平原、盆地、低山丘陵、河谷及融區等，復雜多變的地形地貌給地表變形信息的獲取帶來了諸多困難。近年來持續更新的遙感影像，為獲取面域范圍內精確的地表變形提供了強有力的支持。

圖1 研究區地理位置

2 多年凍土區地表變形

2.1 數據與方法

2.1.1 數據

選用歐空局(ESA)通過數據分發系統(http://scihub.copernicus.eu)提供的空間分辨率為20 m的Sentinel-1A影像數據和美國太空總署(https://nex.nasa.gov.nex/)提供的空間分辨率為30 m的高程數字模型(DEM)，獲取研究區地表變形信息。其中，Sentinel-1A數據的時間為2015年6月—2019年6月(其中2015年11月—2016年5月無數據)，平均每月一景數據，共計43景。所有數據均為降軌數據，VV極化方式，入射角約39°的干涉寬幅工作模式下的TOPS數據。

2.1.2 方法

利用SBAS-InSAR技術處理收集的Sentinel-1A數據，以獲取研究區地表變形。SBAS計算中，先根據時空基線閾值及多普勒頻率差組合生成干涉像對，然后借助DEM數據對各像對逐一進行差分干涉運算，去除總相位中包含的地形相位和其他多余相位，獲取地表變形相位(式(1))，最后利用最小二乘法或奇異分解法，根據變形相位與變形量之間的關系，獲取地表時序變形[16]。圖2為SBAS技術數據處理具體流程。

圖2 SBAS-InSAR數據處理流程

Δφ(x,y)=φdef(x,y)+φDEM(x,y)+φres(x,y)+

(1)

式中：Δφ(x,y)為差分干涉圖中像元的總相位；φdef(x,y)為地表變形引起的相位；φDEM(x,y)為地表高程誤差引起的相位；φres(x,y)為其他因素引起的相位，如平地相位、噪聲相位等；2kπ表示相位為纏繞相位;d(ta,(x,y)) 和d(tb,(x,y))分別為ta、tb時刻影像相對于初始影像在視線向的累積變形量，初始影像變形恒為0；λ為雷達波長;B⊥為垂直基線距;ΔH為高程誤差;S為視線斜距;θ為入射角。

SBAS數據處理過程中，首先根據獲取數據的時空分布規律，設定其空間和時間基線閾值分別為100 m和365 d，生成干涉像對連接圖。為提高監測結果的精度，逐一剔除了依據氣溫數據劃分的凍結期(12月—次年4月)和融化期(5—11月)數據生成的連接圖，僅篩選同時期干涉像對參與變形計算。其次，利用獲取的DEM數據逐一對干涉圖進行差分干涉處理，獲取去除地形相位后的差分干涉相位圖。再次，利用Goldstein濾波法去除噪聲誤差，提高干涉圖清晰度；利用Minimum Cost Flow方法進行相位初步解纏，解纏過程中以青藏公路沿線深入多年凍結層中的橋墩為控制點，考慮到研究區變形量級相對較小，設置解纏閾值為0.2；利用空間低通三角濾波和時間高通線性濾波法去除大氣延遲相位，以獲取精確的地表變形相位信息。最后，利用最小二乘法計算研究區地表變形量，為保證計算結果的精確性，僅對相干系數大于0.5的高相干性像元進行了計算，并利用DEM數據對結果進行地理編碼。對于沒有參與變形計算的低相干點(相干系數小于0.5)，利用克里金插值法進行差值計算[17]，從而得到整個研究區的地表變形分布圖。

2.2 地表變形結果

圖3為利用SBAS-InSAR方法，結合Sentinel-1A數據獲取的研究區地表年變形速率分布，正值表示地表呈現向上的抬升速率，負值表示地表呈現向下的沉降速率。可以看出，研究區地表變形空間分布差異大，年變形速率介于-33～15 mm/a，這可能與地質條件、局地凍土特征等因素有關。但就整個研究區而言，地表有緩慢下沉趨勢，整個研究區地表年變形速率的平均值為-13 mm/a，這與前人基于野外實測數據的分析結果一致[5,18]。

圖3 研究區地表年變形速率分布

此外，從圖3可以看出，研究區大部分地區相對穩定，地表年變形速率較小，在楚瑪爾河高平原、五道梁、沱沱河、通天河及安多等區域，地表存在較大的年變形速率，且以沉降速率為主，這可能與氣溫升高背景下研究區內多年凍土的年平均地溫升高和冰融化等因素有關。

2.3 精度驗證

為驗證InSAR監測結果的精確性，獲取了同時期野外4個監測點(見圖1)的現場實測變形數據。變形測量用沉降桿，觀測采用水準儀進行人工定期觀測。由圖4可以看出，4個監測點現場測量值與InSAR測量值之間的絕對誤差分別為1.4～12.6、1.6～15.3、1.2～9.4、1.3～38.2 mm，平均絕對誤差分別為9.8、7.2、4.3、7.9 mm。可以看出，對于監測點1、2和3，兩種監測結果間的絕對誤差范圍均較小，而監測點4，大部分時間點的絕對誤差較小，個別時間點的絕對誤差較大，可達38.2 mm，這可能與局地降水等因素有關。但總體而言，InSAR監測的多年凍土區地表變形與野外現場監測數據較吻合，4個監測點的平均絕對誤差均小于10 mm，結果可信度高，可用于后續地表變形與影響因素相關性的分析。

圖4 InSAR監測與現場監測地表變形結果對比

3 多年凍土區地表變形影響因素

研究結果顯示，多年凍土區地表總變形主要由活動層的凍融循環變形、多年凍土上限處的融沉變形以及多年凍土層的蠕變變形組成[5]。其中凍融循環變形是多年凍土凍融過程中土結構受冷生作用影響產生的，與活動層的厚度關系最為密切[5]。融沉變形主要是土中冰融化后多年凍土上限下移造成的，與體積含冰量密切相關[18]。蠕變變形是多年凍土升溫導致土物理力學性質改變而產生的，主要與年平均地溫有關[6]。由此可見，活動層厚度、體積含冰量及年平均地溫是影響多年凍土區地表變形的主要因素。鑒于此，主要探討地表變形與這3個影響因素間的相關關系。

3.1 活動層厚度

活動層厚度是指地表至多年凍土上限處的深度，其計算方法眾多[19-21]。本文研究區的活動層厚度分布采用龐強強等根據野外實測數據校正得到的模型[21]計算：

(2)

式中：h為活動層厚度，m；λf為土體導熱系數，W/(m·℃)；L為冰融化潛熱，L=3.3×105J/kg；γck為土體干容重，kg/m3；W為土體中總的體積含水量，%；Wu為土體中的未凍水體積分數，%。各參數具體取值見文獻[21]。I為多年凍土區地表的融化指數(℃·d)，計算過程如下：

(3)

式中：E和N分別為十進制表示的經度(°)和緯度(°)，H為地表高程(m)，tm為地表的年平均溫度(℃)，tc為最冷月份地表的月平均溫度(℃)，tw為最暖月份地表的月平均溫度(℃)。圖5(a)為計算得到的研究區活動層厚度分布。

3.2 年平均地溫

年平均地溫指地溫年變化深度(即多年凍土年較差為零)處的地溫，是反應多年凍土變化動態和劃分凍土帶的主要指標之一[22]。諸多學者根據年平均地溫與緯度、經度、高程等因素間的相關關系，建立了青藏高原年平均地溫計算模型[22-24]。基于GIS平臺，利用經驗模型[24]獲取研究區年平均地溫的分布:

t=50.633-0.830N-0.005H

(4)

式中：t表示年平均地溫(℃)，N和H分別表示緯度(°)和地表高程(m)。計算結果如圖5(b)所示。

圖5 青藏工程走廊地表變形影響因素分布

3.3 體積含冰量

體積含冰量主要是指多年凍土上限附近處的體積含冰量。學者們基于野外實測數據建立了體積含冰量計算模型[25-26]。利用下述模型計算研究區的體積含冰量[26]：

Iv=0.34×ST+0.29×INDV+0.24×SD+0.13×t

(5)

式中：Iv表示體積含冰量(%)；ST表示土質類型，由中國科學院南京土壤研究所提供，并根據土體導熱系數對土質進行賦值[21]；SD表示地表坡度信息，通過美國太空總署提供的30 m精度的數字高程模型(DEM)提取得到；INDV表示地表植被覆蓋歸一化指數，利用Landsat數據借助ENVI平臺提取；t表示年平均地溫，通過式(4)計算。式中ST、SD、INDV和t均為利用GIS平臺得到的歸一化值(無量綱)。根據式(5)對各變量進行疊加，計算的研究區體積含冰量分布圖如5(c)所示。

4 多年凍土區地表變形與影響因素相關性

4.1 相關性分析方法

目前，大時空尺度范圍內多變量相關性研究最常用的方法有簡單相關分析和偏相關分析。簡單相關分析可以確定兩個變量間的線性相關性，但該結果隨其他變量的變化而變化[27]。偏相關分析是在控制或消除其他變量影響的條件下，衡量兩個變量間的凈相關關系[27]。各相關系數具體計算如下。

4.1.1 簡單相關系數計算

(6)

4.1.2 偏相關系數計算

(7)

式(7)也稱一階偏相關系數，用于分析3個變量間的偏相關性。其中rxy.z為將變量z消除后變量x與y之間的偏相關系數，rxy、rxz、ryz分別為變量x與y、x與z及y與z間的簡單相關系數。4個及以上變量間的高階偏相關系數計算公式如下：

(8)

4.1.3 顯著性檢驗

顯著性檢驗用于反映樣本統計量和假設總體參數間的顯著性差異[27]。檢驗中運用小概率原理，事先假定判斷界限，即顯著性水平(α=0.05)。當計算的顯著性概率(P)小于α時，兩個變量間存在顯著相關性，反之兩個變量間無顯著相關性。

4.2 地表變形與影響因素相關性

圖6為逐區域統計計算的各像元地表變形隨影響因素的變化關系。可以看出，研究區地表年變形速率隨體積含冰量增大而增大，在體積含冰量小于10%的區域，年變形速率均值為-8 mm/a，在體積含冰量大于50%的區域，年變形速率均值增至-17 mm/a。同時可以看出，地表年變形速率隨年平均地溫的升高增幅明顯，在年平均地溫低于-2 ℃的區域，地表年變形速率均值為-11 mm/a；在年平均地溫介于-0.5 ～0 ℃時，地表年變形速率均值為-19 mm/a。此外，地表年變形速率隨活動層厚度的增大而增大，但增幅不明顯，在活動層厚度小于2 m的區域，年變形速率均值為-14 mm/a；在活動層厚度大于4 m的區域，年變形速率均值為-18 mm/a。

圖6 研究區地表年變形速率隨影響因素變化

由圖6可以看出，多年凍土區地表年變形速率隨3個影響因素的增大而增大，且隨體積含冰量和年平均地溫的增大幅度較活動層厚度明顯，這與前人得出的年平均地溫和體積含冰量是導致多年凍土區地表融沉的主要原因的結論一致[5-6]。上述分析顯示了地表變形與影響因素間的變化關系，但無法體現地表變形與影響因素的相關程度。為進一步定量分析地表變形與影響因素間的相關性，逐區域逐像元統計計算地表變形與影響因素間的簡單相關系數及偏相關系數。

4.2.1 地表變形與影響因素簡單相關分析

表1～3為研究區地表年變形速率與影響因素在不同區域的簡單相關系數(R)。由表1可以看出，當含冰量高于30%時，多年凍土區地表變形與體積含冰量之間存在密切關系，簡單相關系數大于0.8，相應的P<0.05。而在其他體積含冰量區域相關性較差，簡單相關系數介于0.2～0.6，P均大于0.05。總體而言，多年凍土區地表變形與體積含冰量具有正強相關關系，平均簡單相關系數為0.61。

表1 地表變形與不同體積含冰量簡單相關系數

如表2所示，在年平均地溫高于-1 ℃的區域，地表變形與年平均地溫關系密切，簡單相關系數大于0.8，相應的P<0.05。在年平均地溫低于-1 ℃的區域，二者間的簡單相關系數介于0.2～0.6，P均大于0.05，說明相關性不顯著。與體積含冰量相似，多年凍土區地表變形與年平均地溫具有正強相關關系，平均簡單相關系數為0.64。

表2 地表變形與不同年平均地溫簡單相關系數

由表3可知，多年凍土區地表變形與活動層厚度呈現一定的相關性，在小于3.0 m區間范圍內，二者存在弱相關關系，簡單相關系數介于0.2～0.4。在3.0 m以上區間，二者簡單相關系數介于0.4～0.6，有中等相關性。除4.0 m以上活動層厚度，其他區域的P均大于0.05，表明相關性不顯著。總體而言，多年凍土區地表變形與活動層厚度具有正弱相關關系，平均簡單相關系數為0.38。

表3 地表變形與不同活動層厚度簡單相關系數

4.2.2 地表變形與影響因素偏相關分析

表4～6為地表變形與影響因素在不同區域的偏相關系數(r)。由表4可知，去除活動層厚度和年平均地溫因素的影響后，當體積含冰量高于30%時，多年凍土區地表變形與體積含冰量依然存在極強正相關關系，偏相關系數大于0.8，相應的P<0.05。但與簡單相關分析結果不同的是，在體積含冰量為10%～20%和20%～30%的區域，地表變形與體積含冰量的相關關系顯著，偏相關系數介于0.6～0.8，相比簡單相關系數有明顯的增幅。此外，在體積含冰量低于10%的區域，雖偏相關系數相比簡單相關系數有所增大，但地表變形與體積含冰量的關系仍不顯著，P>0.05。總體而言，去除活動層厚度和年平均地溫的影響后，多年凍土區地表變形與體積含冰量之間存在正強相關關系，平均偏相關系數為0.75。

表4 地表變形與不同體積含冰量偏相關系數

由表5可以看出，去除活動層厚度和體積含冰量因素的影響后，在年平均地溫高于-1 ℃區域，多年凍土區地表變形與年平均地溫仍具有顯著相關關系，偏相關系數大于0.8，相應的P<0.05。在年平均地溫低于-2 ℃的區域，二者間的相關性仍不顯著，偏相關系數介于0.2～0.4，P>0.05。與簡單相關分析結果不同的是，在年平均地溫介于-2～-1 ℃時，地表變形與年平均地溫的相關性明顯增強，二者存在正強相關關系，相比簡單相關系數，偏相關系數增幅明顯，約為0.62。總體而言，去除活動層厚度和體積含冰量影響后，多年凍土區地表變形與年平均地溫間存在正強相關關系，平均偏相關系數為0.70。

表5 地表變形與不同年平均地溫偏相關系數

如表6所示，去除體積含冰量和年平均地溫因素的影響后，多年凍土區地表變形與活動層厚度的相關關系變化不大，在3.0 m以上活動層厚度區間，二者有正中等相關關系，偏相關系數介于0.4～0.6。在3.0 m以下活動層厚度區間，二者相關性不顯著，P>0.05。與簡單相關分析結果不同的是，在活動層厚度為3.0～4.0 m時，多年凍土區地表變形與活動層厚度相關性顯著，P<0.05。總體而言，去除體積含冰量和年平均地溫影響后，多年凍土區地表變形與活動層厚度間存在正中等相關關系，平均偏相關系數為0.42。

表6 地表變形與不同活動層厚度偏相關系數

綜上，在體積含冰量大于30%、年平均地溫高于-1 ℃的多年凍土區，地表變形與影響因素存在強相關性，偏相關系數均大于0.8，P<0.05。在體積含冰量為10%～30%、年平均地溫為-2～-1 ℃、活動層厚度大于3 m的多年凍土區，地表變形與影響因素存在較強的相關性，偏相關系數介于0.4～0.8，P<0.05。而在其他區域，地表變形與影響因素間的相關關系較弱，偏相關系數介于0～0.4，P>0.05。但總體而言，多年凍土區地表變形與體積含冰量和年平均地溫均存在正強相關關系，平均偏相關系數分別為0.75和0.70。而多年凍土區地表變形與活動層厚度間的相關性較差，二者存在正中等相關關系，偏相關系數為0.42。由此可以看出，多年凍土區地表變形與體積含冰量和年平均地溫的相關性較強，與活動層厚度的相關性較弱，且前兩者的相關性明顯強于后者。

對比多年凍土區地表變形與影響因素的簡單相關分析和偏相關分析結果，可以看出二者間存在較大的差異。主要表現為去除其他兩個因素的影響后，地表變形與影響因素凈相關系數(偏相關系數)相比綜合相關系數(簡單相關系數)明顯增大。這也說明了由于不同變量間的復雜相關關系，簡單相關分析在一定程度上增強或削弱了地表變形與影響因素間的相關性，而偏相關分析通過去除多余因素的影響，更好地揭示了多年凍土區地表變形與不同因素間的相關關系。

5 結 論

1)研究區地表變形空間差異大，年變形速率介于-33～15 mm/a，但整個研究區的地表年變形速率均值為-13 mm/a，地表有緩慢下沉趨勢。InSAR監測結果有較高的可信度，與野外4個監測點現場實測數據的平均絕對誤差分別為9.8、7.2、4.3、7.9 mm，均小于10 mm。

2)在體積含冰量大于30%、年平均地溫高于-1 ℃的多年凍土區，地表變形與影響因素存在強相關性。在體積含冰量介于10%～30%、年平均地溫介于-2～-1 ℃、活動層厚度大于3 m的多年凍土區，地表變形與影響因素存在較強的相關性。在其他區域，地表變形與影響因素間的相關性較弱。

3)多年凍土區地表變形與體積含冰量存在正強相關關系，簡單相關系數和偏相關系數均值分別為0.61和0.75。多年凍土區地表變形與年平均地溫也存在正強相關關系，簡單相關系數和偏相關系數均值分別為0.64和0.70。多年凍土區地表變形與活動層厚度相關性較弱，簡單相關系數和偏相關系數均值為0.38和0.42。

4)總體而言，多年凍土區地表變形與體積含冰量和年平均地溫的相關性較強，與活動層厚度的相關性較弱，且前兩者的相關性明顯強于后者。