聯合主/被動微波遙感技術的活動層季節性凍融變化監測研究

寇 曉 康，劉 信 達，李 旭 東

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

0 引言

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

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

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

1 研究區及數據

1.1 研究區概況

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

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

1.2 研究數據

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

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

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

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

2 研究方法

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

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

2.1 SBAS-InSAR方法

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

2.2 去趨勢方法

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

(1)

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

2.3 重采樣方法

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

2.4 地表凍融判別式算法

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

(2)

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

3 結果分析

3.1 季節形變量提取

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

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

圖5 研究區內季節形變量分布Fig.5 Distribution of seasonal deformation in the study area

圖6 QT08站點對應的3.125 km×3.125 km像元范圍內季節形變量平均值時序變化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監測結果對比分析

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

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

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

3.3 PMRS監測結果與活動層土壤溫度對比分析

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

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

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

3.4 AMRS季節形變監測結果與活動層土壤濕度對比分析

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

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

4 結論

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