寇 曉 康，劉 信 達，李 旭 東

(1.石家莊鐵道大學土木工程學院，河北 石家莊 050043；2.道路與鐵道工程安全保障省部共建教育部重點實驗室(石家莊鐵道大學)，河北 石家莊 050043；3.中核三維地理信息工程技術研究中心，河北 石家莊 050043)

0 引言

多年凍土[1]一般分布在高緯度或高海拔地區(qū)，我國青藏高原地區(qū)多年凍土覆蓋區(qū)域占比近70%[2]。天然地面到多年凍土上限(多年凍土的頂板)之間是活動層(又稱季節(jié)融化層)，活動層季節(jié)性變化引起的地表凍脹融沉是凍土區(qū)鐵路、公路、管線、建筑物等產生形變的主要原因之一。近年來，隨著全球不斷升溫[3]，多年凍土上限下移，活動層厚度及凍脹融沉幅度均隨之發(fā)生改變，深刻影響著寒區(qū)基礎設施的運營與安全。因此，開展活動層的凍融變化監(jiān)測研究至關重要。

傳統(tǒng)活動層凍融變化監(jiān)測方法通過埋設傳感器獲取溫濕度數(shù)據，輔助利用精密水準儀、GNSS等傳統(tǒng)測量手段獲取凍融引起的形變信息，得到的點位信息精確，數(shù)據準確度高，但點位布設密度低會產生較大誤差，有些地方因很難抵達而無法監(jiān)測，同時存在觀測周期長、成本高等問題。相比而言，微波遙感穿透力強，且不受云、雨和太陽照度影響[4]，可用于大范圍活動層的凍融變化形變監(jiān)測，相關研究常采用被動微波遙感(Passive Microwave Remote Sensing,PMRS)技術，其基本探測原理為：水和冰的介電特性在微波波段差異巨大，當活動層內的土壤水發(fā)生冰/水相變之后，根據微波傳感器接收的微波輻射亮溫差異，可反推出土壤的凍結/融化狀態(tài)。常用的基于PMRS技術的地表凍融監(jiān)測方法有單/雙指標算法[5-8]、決策樹算法[9]、季節(jié)閾值算法[8,10]、判別式算法[11-13]等，由此生成了多種地表凍融狀態(tài)數(shù)據集[14-16]。上述算法中用到的亮溫數(shù)據多來源于SMMR、SSM/I、SSMIS、AMSR-E、AMSR-2、FY3B-MWRI、SMOS、SMAP等衛(wèi)星傳感器，時間分辨率較高(1～2 d)，但空間分辨率低，單個像元對應的地面單元邊長約為25～40 km。因被動微波穿透能力與活動層內土壤水分含量及探測波段有關，而被動微波一般只能探測到土壤最表層的凍融狀態(tài)信息[17]，對青藏高原地區(qū)而言，利用PMRS探測到的主要是活動層表層的凍融狀態(tài)。在活動層季節(jié)性變化引起的凍脹融沉形變監(jiān)測方面，一般采用主動微波遙感(Active Microwave Remote Sensing,AMRS)中的雷達差分干涉測量(D-InSAR)技術，利用合成孔徑雷達(SAR)兩次觀測所獲得的雷達復數(shù)影像數(shù)據中的相位信息，通過過濾觀測區(qū)域的地形及平地相位反演地表形變，但該技術易受時(空)間失相干和大氣延遲影響，精度難達要求，為此，學者們提出了以永久散射體法(PS)和小基線集(SBAS)方法為代表的時序差分雷達干涉測量 (MT-InSAR) 技術[18-20]，其克服了傳統(tǒng)D-InSAR技術的缺點，理論精度能達到毫米級[18]。目前已將MT-InSAR技術應用于青藏高原的凍土形變監(jiān)測中[21-29]，研究所需數(shù)據多來源于ERS、RADARSAT、ALOS、TerraSAR-X、ENVISAT、Sentinel-1、GF-3等衛(wèi)星，空間分辨率在1 m～1 km之間，時間分辨率在11～46 d之間[30]。

活動層表層的凍融變化與整個活動層的凍融周期緊密相關，伴隨活動層凍脹融沉的變化表層會產生形變,可見二者間有天然聯(lián)系，但受制于PMRS的低空間分辨率及AMRS的低時間分辨率，現(xiàn)有研究鮮有從遙感視角分析二者間的關聯(lián)。鑒于此，本文以青藏高原五道梁地區(qū)為研究區(qū)，擬聯(lián)合主/被動微波遙感技術獲取地表凍融與形變信息，并結合地面實測數(shù)據深入剖析活動層季節(jié)性凍融變化中凍融狀態(tài)、土壤溫度、地表形變、土壤濕度間的內在關聯(lián)，從而明確主/被動微波遙感監(jiān)測信息與活動層內部溫濕度參量的協(xié)同變化機制，以期從遙感視角深入揭示活動層的季節(jié)性變化過程。

1 研究區(qū)及數(shù)據

1.1 研究區(qū)概況

五道梁地區(qū)位于青藏高原腹地青海省西南部曲麻萊縣境內，海拔4 600～4 800 m，屬高原山地氣候，四季皆冬,1月平均氣溫為-16.2 ℃，7月平均氣溫為6.0 ℃，全年平均氣溫為-5.1 ℃，是全國最低值。研究區(qū)南部地表覆蓋類型主要為高寒草甸，北部主要為裸地，該地區(qū)還分布著山體、熱融湖塘等，青藏鐵路貫穿其中(圖1)，研究該地區(qū)地表凍融變化及形變對于凍土區(qū)鐵路的維護具有重要意義。

圖1 青藏高原五道梁地區(qū)地理位置及地表覆蓋類型Fig.1 Geographical location and land cover types of Wudaoliang area in the Qinghai-Tibet Plateau

1.2 研究數(shù)據

(1)AMRS數(shù)據，來源于Sentinel-1衛(wèi)星系統(tǒng)，下載網址https://search.asf.alaska.edu/#/，該系統(tǒng)由Sentinel-1A和Sentinel-1B兩顆SAR衛(wèi)星構成，單顆星重訪周期為12 d，雙星協(xié)同工作的重訪周期為6 d。2017年10月至2018年12月研究區(qū)Sentinel-1B衛(wèi)星數(shù)據缺失較嚴重，考慮到數(shù)據的一致性，本研究選取Sentinel-1A 衛(wèi)星在干涉寬幅模式(IW)下的C波段VV極化方式升軌影像，共計36景，影像距離向分辨率為5 m，方位向分辨率為20 m，空間范圍如圖1所示。

(2)PMRS數(shù)據，來源于美國冰雪數(shù)據中心發(fā)布的增強型空間分辨率的亮度溫度(CETB)數(shù)據集，該數(shù)據集包含SSM/I、SSMIS、AMSR-E、AMSR-2等傳感器的長時序亮溫數(shù)據，且各數(shù)據源之間進行了一致性校準，下載網址https://nsidc.org/data/NSIDC-0630/versions/1。選取DMSP F18衛(wèi)星SSMIS傳感器2017年10月至2018年12月的升軌數(shù)據，時間間隔為1 d,數(shù)據采用網格化存儲，投影方式為EASE-Grid 2.0。不同通道的空間分辨率有區(qū)別，其中19 GHz亮溫的分辨率為6.25 km，37 GHz亮溫的分辨率為3.125 km。本文首先對19 GHz亮溫進行降尺度，使之與37 GHz相同，再進行凍融判別。

(3)地面觀測數(shù)據，來源于青藏高原多年凍土綜合監(jiān)測數(shù)據集[31,32]，可從國家青藏高原科學數(shù)據中心下載，選取位于研究區(qū)編號為QT08 2017年10月至2018年12月的站點數(shù)據，該站點對表層以下10 cm、40 cm、120 cm、200 cm、240 cm深度處的土壤溫度和土壤濕度進行了長期觀測。

(4)其他輔助數(shù)據，包括30 m空間分辨率的地表分類數(shù)據集(http://www.globallandcover.com/)以及30 m空間分辨率的SRTM DEM數(shù)據集(https://lpdaac.usgs.gov/)。

2 研究方法

本文技術流程(圖2)為：首先利用Sentinel-1 SLC雷達影像數(shù)據，基于SBAS-InSAR監(jiān)測方法獲取研究區(qū)內的地表時序形變信息，并將其作為線性形變分量(年際形變)與周期性形變分量(季節(jié)形變)的組合，通過逐像元線性趨勢擬合及趨勢剔除，獲取研究區(qū)地表季節(jié)形變信息，在此基礎上通過空間聚合提取站點所對應微波像元(約3 km)的平均時序形變信息；其次利用SSMIS亮溫數(shù)據，基于凍融判別式算法，通過對亮溫重采樣獲取研究區(qū)內約3 km×3 km像元尺度的每日地表凍融狀態(tài)信息；最后結合地面站點觀測數(shù)據，分析活動層季節(jié)性凍融變化過程中的內部物理參量變化機制及其與外在遙感特征之間的關聯(lián)。

圖2 技術流程Fig.2 Technical flow chart

2.1 SBAS-InSAR方法

SBAS-InSAR方法是對相干目標進行相位分析以獲取時序形變，通過選擇合適的空間基線和時間基線閾值組成差分干涉對，然后選取相干目標點利用線性相位變化模型進行建模和解算，通過時空濾波去除大氣延遲，在減少D-InSAR處理過程中的失相干影響及高程、大氣誤差的同時獲取地表形變時間序列[20]。本文SBAS-InSAR方法實現(xiàn)流程為：根據研究區(qū)范圍對衛(wèi)星影像進行剪裁，將生成的數(shù)據進行干涉像對配對，加入外部DEM數(shù)據去除地形相位；依據小基線集技術原理，設置最大臨界基線百分比為2%、時間基線為60 d，選擇20～30個相干性較好的點作為相位解纏的參考點；將干涉圖進行自適應濾波，使用最小費用流(MCF)方法進行相位解纏和基線精化估計；對殘余的相位進行自適應濾波和第二次解纏，對輸入的干涉圖進行優(yōu)化，在第一步得到的形變速率基礎上進行大氣濾波，從而估算和去除大氣相位，最后將結果編碼到WGS-84坐標系下，得到時序形變結果。

2.2 去趨勢方法

多年凍土區(qū)地表形變分為年際尺度形變與季節(jié)尺度形變：在年際尺度上，夏季多年凍土區(qū)地表差異主要由多年凍土層上限的上下移動導致；在季節(jié)尺度上，地表形變主要由活動層的周期性凍脹融沉導致，即夏季活動層中的冰融化成水，地表發(fā)生下沉，冬季活動層水分重新凝結成冰，地表發(fā)生隆升，呈現(xiàn)周期性的形變特征[33]。本研究假設年際形變量在1年內為線性形變，可通過方程(Y=K×X+B，Y為年際形變量，X為時間(年)，K和B為線性擬合系數(shù))進行擬合，基于該方程可獲取各像元的年際變化趨勢；假定季節(jié)形變量為原始形變量在年際變化趨勢線上的投影轉換分量D，計算公式為：

(1)

式中：X0、Y0為原始形變量在以時間和形變量組成的平面坐標系中的坐標。

2.3 重采樣方法

本文采用雙三次插值方法將SSMIS 19 GHz亮溫數(shù)據的空間分辨率從6.25 km降至3.125 km，該方法用三階多項式函數(shù)通過相鄰的網格點擬合曲面，最后根據曲面擬合函數(shù)得到查詢點的插值結果。具體計算參考文獻[34]。

2.4 地表凍融判別式算法

凍融判別式算法比其他國際上主流算法的判對率高[14]，對裸地、植被或積雪覆蓋下的土壤凍融情況均能進行有效區(qū)分，因此，本研究采用該算法進行活動層表層凍融狀態(tài)的識別。該算法包括36.5 GHz V極化下的亮溫(TB36.5V)和18.7 GHz H極化與36.5 GHz V極化下的亮溫比(Qe18.7H/36.5V)兩個主要指標，分別用于指示地表溫度和土壤發(fā)射率的變化，具體判別公式為：

(2)

式中：DF、DT分別為凍土和融土的判別方程函數(shù)值,若DF-DT>0，則判別為凍土，否則判別為融土。

3 結果分析

3.1 季節(jié)形變量提取

3.1.1 假定合理性分析 本研究假定基于SBAS-InSAR方法獲取的形變結果由線性形變分量(年際形變量)與周期性形變分量(季節(jié)形變量)兩部分組成。由于研究期較短，為驗證年際形變趨勢是否存在，首先對時序形變結果進行逐像元線性擬合，統(tǒng)計各像元時序形變線性擬合方程的顯著性檢驗結果(表1)并進行可視化(圖3)，進而判定年際形變趨勢的穩(wěn)定性。由表1可知，顯著相關區(qū)域面積之和占統(tǒng)計總面積的77.41%，表明研究區(qū)內絕大部分區(qū)域存在較明顯的年際變化趨勢，證明了本文假設的合理性。進一步對QT08站點3.125 km范圍內的像元形變平均值進行統(tǒng)計(圖4)，可以看出，其地表形變既表現(xiàn)出周期性特征，又表現(xiàn)出一定的年際沉降線性趨勢。由于常年積雪等因素影響，部分區(qū)域失相干嚴重，造成形變結果存在少量的空白區(qū)(圖3中的白色區(qū)域)，這部分并未參與計算。

圖3 研究區(qū)逐像元時序形變線性擬合方程的顯著性檢驗結果空間分布Fig.3 Spatial distribution of significance test results of linear fitting equation for time series deformation of pixels in the study area

圖4 QT08站點對應的 3.125 km×3.125 km像元范圍內形變平均值的時序變化Fig.4 Time series variation of average deformation in one 3.125 km×3.125 km pixel corresponding to QT08 station

表1 顯著性檢驗標準及結果Table 1 Statistics of significance test criteria and results

3.1.2 季節(jié)形變量分布 利用2.2節(jié)的去趨勢方法獲取研究區(qū)2017年10月至2018年12月之間36景影像所對應的季節(jié)形變量的空間分布(圖5)，可以看出，研究區(qū)大部分區(qū)域表現(xiàn)出較明顯的季節(jié)性凍脹融沉趨勢，其中北部尤為明顯，南部稍弱，可能由于研究區(qū)北部以裸地為主，而南部以草地為主(圖1)，受植被干擾產生失相干現(xiàn)象，進而對反演結果產生一定影響。進一步對QT08站點周圍3.125 km的像元網格的季節(jié)形變量平均值進行可視化(圖6)，可以看出，經去趨勢處理后，該區(qū)域內的地表形變具有明顯的凍脹融沉周期性變化趨勢。

圖5 研究區(qū)內季節(jié)形變量分布Fig.5 Distribution of seasonal deformation in the study area

圖6 QT08站點對應的3.125 km×3.125 km像元范圍內季節(jié)形變量平均值時序變化Fig.6 Time series variation of average seasonal deformation in one 3.125 km×3.125 km pixel corresponding to QT08 station

3.2 AMSR和PMSR監(jiān)測結果對比分析

本研究利用SSMIS亮溫數(shù)據，在對19 GHz亮溫進行降尺度處理的基礎上，基于凍融判別式算法最終獲取研究區(qū)范圍內3.125 km分辨率的每日地表凍融狀態(tài)信息，然后提取QT08站點所對應的被動微波像元凍融狀態(tài)信息，并與該地區(qū)的地表形變信息進行對比(圖7)。可以看出，凍融狀態(tài)及形變均表現(xiàn)出明顯的周期特征，為更清晰地對比二者間的聯(lián)系與差異，進行分時段對比分析。1)凍結期：自2017年10月18日開始，PMRS監(jiān)測結果顯示該區(qū)域地表經歷一次凍融交替后持續(xù)處于凍結狀態(tài)，說明自此進入凍結期；AMRS監(jiān)測結果顯示，隨著地表凍結狀態(tài)的持續(xù)，地表不斷抬升，且有慢—快—慢—平穩(wěn)的變化特征，而該特征在PMRS凍融狀態(tài)監(jiān)測信息中并未體現(xiàn)。2)凍融交替期：自2018年4月初開始，PMRS監(jiān)測結果顯示該區(qū)域地表進入凍融交替期；AMRS監(jiān)測結果對應出現(xiàn)起伏特征，但變化量很小。3)融化期：自2018年5月中旬開始，PMRS監(jiān)測結果顯示除了偶爾的凍融交替現(xiàn)象之外(考慮可能為誤判)，該區(qū)域地表始終保持融化狀態(tài)，說明進入穩(wěn)定的融化期；AMRS監(jiān)測結果顯示，隨著地表融化狀態(tài)的持續(xù)，地表不斷沉降，沉降趨勢同樣有慢—快—慢—平穩(wěn)的特征，但該特征在PMRS凍融狀態(tài)監(jiān)測信息中并未體現(xiàn)。4)再次凍融交替期：自2018年9月初開始，PMRS監(jiān)測結果顯示該區(qū)域再次進入凍融交替期；AMRS監(jiān)測結果顯示形變量再次出現(xiàn)波動，且波動幅度較小。5)再次凍結期：自2018年10月中旬之后，該區(qū)域地表持續(xù)處于凍結狀態(tài)，地表不斷抬升，說明再次進入凍結期。

圖7 地表凍融狀態(tài)與季節(jié)形變量的時序變化對比Fig.7 Temporal variation comparison of surface freezing-thawing state and seasonal deformation

綜上，在活動層季節(jié)性變化的各個周期內，PMRS地表凍融狀態(tài)監(jiān)測結果與AMRS季節(jié)形變量提取結果在時間上均具有較強的關聯(lián)性與一致性，但值得注意的是，在凍結期和融化期內，AMRS所顯示的地表形變速率特征并未在PMRS監(jiān)測結果中體現(xiàn)。

3.3 PMRS監(jiān)測結果與活動層土壤溫度對比分析

有學者依據活動層凍融過程不同階段的水熱特征差異，將青藏高原活動層的年際變化過程劃分為夏季融化、秋季凍結、冬季降溫和春季升溫4個階段[35]。為明晰主/被動微波遙感技術在活動層季節(jié)性凍融變化監(jiān)測中的作用，本研究將年際凍融變化過程重新劃分為穩(wěn)定凍結期、穩(wěn)定融化期和兩次凍融交替期。由于被動微波的穿透能力有限，PMRS的監(jiān)測結果主要反映近地表約0～5 cm的土壤凍融狀態(tài)，現(xiàn)結合QT08站點在10 cm、40 cm、120 cm、200 cm和240 cm深度處的實測溫度進行對比分析(圖8)。從圖8可以看出，基于PMRS獲取的凍融狀態(tài)變化情況與各層溫度之間既存在明顯的相關性，又存在一定的差異性，具體到各時期總結如下：1)凍結期：PMRS監(jiān)測結果顯示，兩次地表凍結期分別自2017年10月18日和2018年10月18日左右開始，而從QT08站點各層土壤溫度看，淺層土壤(40 cm以內)凍結開始時間與其大體相當，伴隨外界低溫向下傳導，凍結起始時間隨活動層深度的增加有不斷向后延遲趨勢。2)兩次凍融交替期：PMRS監(jiān)測結果顯示，2018年4月初及9月下旬的凍融交替現(xiàn)象主要出現(xiàn)在活動層淺層，120 cm及以下溫度振幅較小，幾乎無明顯的凍融交替現(xiàn)象。由于活動層淺層厚度較小，推斷其為凍融交替期季節(jié)形變幅度較小的原因。3)融化期：PMRS監(jiān)測結果顯示，自2018年5月初表層融化狀態(tài)基本保持穩(wěn)定，活動層淺層土壤的融化期基本與PMRS監(jiān)測結果一致，隨著外界熱量不斷向下傳導，雖然各層土壤溫度均有抬升，但中層(120 cm)及以下土壤溫度仍維持在0℃以下，其凍結期并未結束，隨后經過約1個月的熱量向下傳導，這些層次才相繼進入融化期。

圖8 地表凍融狀態(tài)與各層溫度時序變化對比Fig.8 Temporal variation comparison of surface freezing-thawing state and soil temperature for various layers

綜上可知，PMRS監(jiān)測結果顯示的凍結期/融化期的開始和結束時間與活動層淺層凍結期/融化期的開始和結束時間基本一致，活動層中層、深層凍結期/融化期的開始和結束時間相對滯后，其原因在于淺層土壤與外界環(huán)境交互較劇烈，而熱量在土壤中的傳遞需要時間；此外，PMRS監(jiān)測結果中的凍融交替現(xiàn)象主要出現(xiàn)在活動層的淺層，中層和深層土壤溫度振幅較小，無明顯凍融交替現(xiàn)象，因此凍融交替期內AMRS獲取的季節(jié)形變幅度一般也較小。

3.4 AMRS季節(jié)形變監(jiān)測結果與活動層土壤濕度對比分析

AMRS監(jiān)測到的形變實際為活動層各層冰/水相變的積累，本文結合QT08站點在不同層次的實測土壤濕度進行對比分析。從圖9可以看出，基于AMRS獲取的季節(jié)形變量與各層濕度之間有明顯的相關性，具體到各時期總結如下：1)凍結期：AMRS監(jiān)測結果顯示，自2017年10月中旬至2018年2月中旬期間，地面整體出現(xiàn)抬升趨勢，對比各層土壤濕度變化可知，地面抬升主要與各層水分持續(xù)下降有關，冰/水相變(相變水量)速率決定了抬升趨勢的變化速率，自2018年2月下旬開始，當水分全部凍結后，地表保持穩(wěn)定狀態(tài)。2)兩次凍融交替期：AMRS監(jiān)測結果顯示，在2018年4月初及9月中旬，地表形態(tài)出現(xiàn)微小波動，這主要與淺層土壤的冰/水相變有關，在3.3節(jié)中已做出解釋。3)融化期：AMRS監(jiān)測結果顯示,在2018年4月中旬至9月初期間，地面整體出現(xiàn)沉降趨勢，對比各層土壤濕度變化可知，地面沉降主要與各層水分的回升有關，水分的回升速率與地表形變速率呈正相關，個別層次土壤濕度突然增高可能是降雨所致。綜上可知：AMRS監(jiān)測得到的季節(jié)形變量與活動層各層冰/水相變總量有關，季節(jié)形變速率與冰/水相變速率呈正相關，當活動層內部冰/水不再發(fā)生變化時，季節(jié)形變趨于穩(wěn)定。

圖9 季節(jié)形變反演結果與各層土壤水分時序變化情況對比Fig.9 Comparison between inversion results of seasonal deformation and temporal variation of soil moisture in each layer

4 結論

本研究以青藏高原五道梁地區(qū)為研究區(qū)，分別基于主/被動微波遙感技術獲取了3 km×3 km分辨率的每日活動層表層凍融狀態(tài)信息及5 m×20 m分辨率、時間間隔約12 d的活動層季節(jié)形變信息；在此基礎上，通過空間重采樣及時序擬合技術，結合地面實測土壤溫度、濕度數(shù)據，對比分析了凍結期、融化期和凍融交替期內活動層各凍融參量之間的內在聯(lián)系，得出以下結論：1)PMRS凍融狀態(tài)監(jiān)測結果可用于提取活動層表層凍融周期，該周期與活動層淺層凍融周期基本一致，活動層中層及深層的凍融起始時間相對滯后；2)AMRS季節(jié)形變監(jiān)測結果與活動層內部總的相變水量緊密相關，形變速率與冰/水相變速率呈正相關關系；3)凍融交替現(xiàn)象主要發(fā)生在活動層的淺層，深層不明顯，此為AMRS監(jiān)測結果波動幅度較小的原因；4)PMRS監(jiān)測結果與AMRS監(jiān)測結果在時間上有較強的關聯(lián)性和一致性，在凍融周期上的統(tǒng)一是季節(jié)形變監(jiān)測中主/被動微波遙感協(xié)同應用的基礎。本文關于AMRS和PMRS的聯(lián)合分析對于大尺度范圍內的活動層季節(jié)性變化監(jiān)測研究具有重要意義。由于研究區(qū)內可用的站點數(shù)據較少，通過綜合分析本文只得出一些定性結論，未來可結合更多實測數(shù)據，將AMRS與PMRS技術深入結合，應用于活動層凍融參量的定量監(jiān)測中。