楊 通，張智杰，王效科※，郭旭東，于 瀟

（1.中國科學(xué)院生態(tài)環(huán)境研究中心，城市與區(qū)域生態(tài)國家重點實驗室，北京 100085；2.中國科學(xué)院空天信息研究院，北京 100094；3.中國國土勘測規(guī)劃院，自然資源部土地利用重點實驗室，北京 100035）

0 引 言

土壤侵蝕是中國干旱半干旱地區(qū)面臨的主要生態(tài)問題之一，嚴重的土壤侵蝕一方面迅速惡化生態(tài)環(huán)境，另一方面加劇當(dāng)?shù)厝嗣袢罕姷呢毨С潭萚1]?，F(xiàn)階段土壤侵蝕監(jiān)測方法包括建立觀測站、樣地調(diào)查、立體攝影、高精度GPS、三維激光測量、元素示蹤、模型模擬以及遙感監(jiān)測等[2-4]。遙感方法具有宏觀性、綜合性、可視性、實用性、經(jīng)濟性等特點，突破了傳統(tǒng)的樣地調(diào)查微觀機理研究范疇，將研究層面拓展到面上研究和動態(tài)研究[5]，對土壤侵蝕研究發(fā)揮重要作用[6]?，F(xiàn)行遙感監(jiān)測土壤侵蝕方法屬于光學(xué)遙感范疇，監(jiān)測地表正射水平面的二維特征，以像元流失作為土壤侵蝕的主要依據(jù)[7]。然而，土壤侵蝕的二維變化特征并不顯著，光學(xué)遙感傳感器的分辨率也不足以捕捉以毫米計量的土層剝蝕特征，從而低估土壤侵蝕的實際情況，誤以為土壤侵蝕很少或侵蝕程度較輕。

土壤侵蝕的主要特征是土層厚度的剝蝕，是高程向的累計形變。水利部發(fā)布的《土壤侵蝕分類分級標準SL190—2007》規(guī)定：水力侵蝕的土壤平均流失厚度介于5.9～11.1 mm/a，其侵蝕強度為極強烈；風(fēng)力侵蝕的風(fēng)蝕厚度介于10～25 mm/a，其侵蝕強度為中度。理論上，若高程向監(jiān)測精度達到厘米級，即可監(jiān)測到中度以上風(fēng)蝕和極強烈水蝕。

雷達遙感對地物的高程變化信息敏感，監(jiān)測尺度由二維空間轉(zhuǎn)換為三維空間，高程向監(jiān)測精度可達厘米甚至毫米級[8]，理論上具備識別土壤侵蝕的能力。其中，永久散射體合成孔徑雷達干涉測量 PS-InSAR（Persistent Scatterer Synthetic Aperture Radar Interferometry）是雷達微波成像和電磁波干涉的集成技術(shù)[9]，是一種新的D-InSAR（Differential SAR Interferometry）技術(shù)。該技術(shù)通過合成孔徑雷達在不同位置對觀測區(qū)域的相位差等信息分析處理，獲取觀測區(qū)域的三維空間信息。該技術(shù)以其覆蓋范圍大、不受云霧影響、厘米/毫米級監(jiān)測精度的特點，成為監(jiān)測地表形變最有前景的空間大地測量技術(shù)之一，并廣泛應(yīng)用于火山活動[10]、地震[11]、滑坡[12]、不均勻沉降[13]、濕地變化[14]、地形測繪等研究領(lǐng)域[15-16]。PS-InSAR通過分析時間序列影像，提取時序內(nèi)相位和幅度保持穩(wěn)定的離散點，即PS點（Persistent Scatterer），反演地表緩慢形變[17]。與傳統(tǒng)的 InSAR技術(shù)相比，PS-InSAR克服了空間去相干、時間去相干、大氣延遲，顯著提高了測量精度和信噪比，其高程向測量精度優(yōu)于厘米級甚至達到毫米級[8,12]。

目前，國內(nèi)外少見PS-InSAR監(jiān)測土壤侵蝕的研究，有學(xué)者指出土壤侵蝕造成了上海市的地面沉降[2]，也有學(xué)者注意到浙江四明山自然地表在雨季發(fā)生沉降，推測為土壤侵蝕所致[18]；國外在伊朗西部[19]、意大利中部[20]有應(yīng)用InSAR技術(shù)監(jiān)測土壤侵蝕的研究案例。本研究擬利用哨兵一號衛(wèi)星平臺數(shù)據(jù)（Sentinel-1A），基于改進后的PS-InSAR方法，對內(nèi)蒙古和林格爾縣進行土壤侵蝕的面上觀測，結(jié)合多源數(shù)據(jù)對監(jiān)測結(jié)果進行評估和驗證，檢驗PS-InSAR技術(shù)監(jiān)測土壤侵蝕的可行性，為該方法在土壤侵蝕監(jiān)測領(lǐng)域的推廣提供助益。

1 研究方法與數(shù)據(jù)

1.1 研究區(qū)概況

研究區(qū)位于內(nèi)蒙古自治區(qū)中南部（39°58′～40°41′N，111°26′～112°18′E），總面積 3 436 km2，海拔 1 400～2 028 m，年平均降雨量為 376 mm，年均蒸發(fā)量為2 073 mm。和林格爾縣屬內(nèi)蒙古高原和黃土高原的過渡地帶，東南部山區(qū)、黃土丘陵區(qū)的面積占縣域面積的77.7%；西北部屬土默川平原，面積占縣域面積的22.3%。研究區(qū)超過一半?yún)^(qū)域?qū)儆邳S土高原。黃土高原是世界上土壤侵蝕最嚴重區(qū)域，每年入黃泥沙量達 16億 t，且每年由于人類活動新增水土流失面積900～1 100 km2。自大規(guī)模農(nóng)耕活動開始，人類活動就已經(jīng)與水力、風(fēng)力并列，成為土壤侵蝕的主要營力[1]。根據(jù)《土壤侵蝕分類分級標準SL190—2007》分類體系，研究區(qū)既屬于西北黃土高原水力侵蝕區(qū)，又屬于內(nèi)蒙古高原草原中度風(fēng)蝕水蝕區(qū)，水力侵蝕與風(fēng)力侵蝕并存。

1.2 數(shù)據(jù)概況

監(jiān)測數(shù)據(jù)為2017至2018年的17期Sentinel-1A雷達數(shù)據(jù)，每期1景，用于監(jiān)測土壤侵蝕。Sentinel-1A的軌道類型為近圓形的太陽同步軌道，其軌道高度為693 km，傾角為98.18°。衛(wèi)星重訪周期為12 d，內(nèi)有175個軌道。為定向監(jiān)測土層厚度變化，選用干涉寬模式（Interferometric Wide，IW）下的S1 TOPS-modes單視復(fù)數(shù)SLC（Single Look Complex）數(shù)據(jù)，幅寬250 km，入射角29.1°～46°。該數(shù)據(jù)為C波段的升軌數(shù)據(jù)，C波段屬于長波，可以穿透稀疏植被直接作用于地表。SLC數(shù)據(jù)處理后的距離向分辨率為2.3 m，方位向分辨率為13.9 m。

輔助數(shù)據(jù)為同步觀測的哨兵2號（Sentinel-2A）數(shù)據(jù)，用于主被動遙感數(shù)據(jù)協(xié)同觀測；DEM數(shù)據(jù)用于控制地形影響；土地利用數(shù)據(jù)用于掩膜不透水面；降雨量和土壤數(shù)據(jù)用于信息復(fù)合。以雷達數(shù)據(jù)序列為基準數(shù)據(jù)集，數(shù)據(jù)源見表1。

表1 數(shù)據(jù)源一覽表Table 1 Data source list

1.3 研究方法

本研究的總體技術(shù)方案如圖1所示，主要技術(shù)流程為：1）預(yù)處理流程，包括數(shù)據(jù)導(dǎo)入、精軌數(shù)據(jù)更新、主影像選擇、DEM主輔影像配準、ESD配準，burst拼接，子條帶拼接等；2）差分干涉流程[21]，包括ROI截取、干涉相位計算、去除平地和地形相位、生成差分干涉圖等；3）PS點提取流程，包括累計形變點（即PS點）候選點提取、構(gòu)建Delaunay三角網(wǎng)、選擇參考點求取PS點的參數(shù)信息、求解二次差分相位模型、殘余相位解纏、大氣相位估計、累計形變量計算、地理編碼等[9]；4）多尺度濾波流程，包括空間域的不透水面掩膜、穩(wěn)定地表分離，時間域的畸變?yōu)V波、線性濾波等；5）侵蝕特征分析流程，包括土壤侵蝕的空間特征、時間特征等。

需要指出的是，PS-InSAR方法監(jiān)測到的地表累計形變，屬于多種人類擾動和自然因素的混合效應(yīng)，如何從混合效應(yīng)中分離出土壤侵蝕的貢獻，將監(jiān)測目標錨定于土壤侵蝕，是方法有效性的決定因素。本研究對PS-InSAR方法進行了針對性的改進，基于土壤侵蝕的物候特征，設(shè)計了空間域和時間域的多尺度濾波環(huán)節(jié)，定向提取土壤侵蝕造成的累計形變。

1.3.1 數(shù)據(jù)預(yù)處理

首先，生成 SLC數(shù)據(jù)，同時利用各子條帶及 burst的重疊信息生成mosaic SLC和mosaic MLI（Multi Looks Intensity），多視參數(shù)默認為 10和 2。其次，利用 ESA（European Space Agency）提供的精軌 POD數(shù)據(jù)，（https://qc.sentinel1.eo.esa.int/）修正軌道信息查詢列表，減小定時及軌道誤差對干涉相位的影響。再次，選擇同時位于時間基線中部，及時序影像多普勒質(zhì)心頻率中心的時段的數(shù)據(jù)作為主影像（表2）。然后，將 DEM與主影像配準，得到每個像元的三維坐標信息。之后，利用ESD方法生成主影像與副影像的幾何關(guān)系（查詢列表），確保配準精度達到1/1 000像素的要求[22]。最后，將分塊存儲的Sentinel-1A數(shù)據(jù)根據(jù)研究區(qū)的位置提取出來，以便縮小處理數(shù)據(jù)量加快處理時間。此外，對光學(xué)遙感數(shù)據(jù)、DEM、土地利用、土壤、地面觀測站數(shù)據(jù)等，執(zhí)行校正和配準等處理步驟，生成多源可比數(shù)據(jù)。

表2 合成孔徑雷達數(shù)據(jù)序列Table 2 Time series synthetic aperture radar data set

1.3.2 差分干涉處理

主影像分別與輔助影像生成17期差分干涉圖，并結(jié)合衛(wèi)星軌道參數(shù)、成像模型、參考DEM，去地形相位和平地相位。每一個干涉對都需要進行配準和重采樣，形成干涉相位。利用軌道和外部的DEM去除平地效應(yīng)和地形效應(yīng)，生成時間序列的干涉相位。具體包括感興趣區(qū)選擇、干涉處理、地形相位去除、差分模型構(gòu)建 4個步驟。主輔影像干涉相位可表示為公式（1）

式中φtopo為地形相位；φdef為形變相位；φatm為大氣延遲相位；φflat為平地相位；φnoise為噪聲相位。

1.3.3 提取永久散射體

PS候選點選?。哼x擇后向散射較強，并且在時序上較穩(wěn)定的地物目標，如建筑物與巖石等，將相干性差或者一些噪聲點去除，留下高相干且穩(wěn)定的點進行分析。研究區(qū)有150個行政村，并有蒙牛、蒙草等大型企業(yè)，建成區(qū)的人工建筑可以部分替代強散射點的作用。為了保證PS候選點密度，本文采用平均相干性閾值法選點，包含一些中等相干性點目標，以增加PS候選點的選點數(shù)量。

構(gòu)建三角剖分網(wǎng)：利用Delauany三角化方法構(gòu)建三角剖分網(wǎng)絡(luò)，以解決PS離散點的圖像重構(gòu)問題。同時，利用三角網(wǎng)絡(luò)的冗余測量邊，以區(qū)域網(wǎng)絡(luò)平差方法控制形變測量誤差，削弱大氣相位影響。Sentientl-1A時空基線較短，軌道誤差較小，通過構(gòu)建三角網(wǎng)及相鄰PS點的相位雙差模型等方法，可以削弱大氣相位的干擾。

計算二次差分相位：基于PS點的三角剖分網(wǎng)，計算網(wǎng)絡(luò)上每條邊的相鄰PS 點的二次差分相位。通過求解二次差分相位模型，獲得PS點對的形變速率差和DEM修正差。

相位解纏：相位解纏后可得每個PS點的絕對形變速率和DEM修正量，然后可得差分相位的形變相位和殘余地形相位。

大氣濾波：大氣延遲相位在時間域上表現(xiàn)為隨機高頻信號，在空間域上則表現(xiàn)為平滑低頻信號；非線性形變在時間域上表現(xiàn)為低頻信號，在頻率域表現(xiàn)為低頻信號；噪聲相位則為時域高頻，頻域高頻。根據(jù)各相位成分的特點，通過時空域高通、低通濾波分離各相位成分，從而得到形變相位，再由地理編碼得出累計形變量。

1.3.4 多尺度濾波

地表累計形變是多種因素共同驅(qū)動的疊加效應(yīng)，包括人類擾動和自然變化。自然變化又包含植被生長、土壤侵蝕、滑坡、泥石流、凍土、巖石風(fēng)化搬運、枯落物堆積、地質(zhì)運動等多種類型。本研究設(shè)計了空間域和時間域的多尺度濾波，剝離人類擾動和自然變化影響，將監(jiān)測對象錨定為土壤侵蝕引發(fā)的地表形變。

空間域濾波又分為植被、穩(wěn)定地表濾波、人類活動濾波。首先，通過計算主末影像的相位標準差，加入植被指數(shù)的限制因素，振幅小的信號通過，振幅大的剔除，濾除茂密植被，將感興趣區(qū)限制于植被稀疏區(qū)域[23]。雷達C波可以穿透稀疏植被，直接作用于地表。其次，設(shè)置穩(wěn)定地表值域，對永久散射體的累計形變信號進行高通濾波，濾除無變化或變化極小的區(qū)域。然后，根據(jù)土地利用數(shù)據(jù)掩膜城市村鎮(zhèn)的不透水面及活躍耕地，并設(shè)置閾值分割上限，對永久散射體的形變信號進行低通濾波，濾除人類活動干擾。

時間域濾波分為線性濾波和畸變?yōu)V波。不同自然因素驅(qū)動的地表形變在時間序列上有其規(guī)律性，如土壤侵蝕表現(xiàn)為雨季的非線性形變；滑坡、泥石流表現(xiàn)為突發(fā)性形變；地質(zhì)隆起或沉降表現(xiàn)為平滑線性形變；凍土表現(xiàn)為雙向形變等。根據(jù)各種自然因素在時間序列上的形變方向、形變區(qū)間、形變幅度差異，可以初步區(qū)分主要形變類型；再結(jié)合高分光學(xué)遙感數(shù)據(jù)的光譜信息、幾何信息及動態(tài)變化進行檢驗。選取侵蝕特征顯著的樣本點，統(tǒng)計其時序變化參數(shù)，迭代更新初始參數(shù)，修正分離區(qū)間和分離閾值，進而實現(xiàn)土壤侵蝕的定向識別。

2 結(jié)果與分析

2.1 形變速度場

基于時序差分干涉結(jié)果，提取后向散射系數(shù)較高、并且在時序上較穩(wěn)定的永久散射體。對PS候選點進行迭代更新、解纏等處理后，在研究區(qū)共計提取到 1 007 958個累計形變點，樣本量較大，信息豐富。為抵消衛(wèi)星在軌運動時產(chǎn)生的多普勒效應(yīng)，PS-InSAR運行過程均在斜距-零多普勒坐標系中完成，最后再通過地理編碼進行可視化處理，將坐標系統(tǒng)轉(zhuǎn)換為地理坐標系，生成可比數(shù)據(jù)。對1 007 958個累計形變點的平均形變速度進行IDW（Inverse Distance Weighting）插值，計算本區(qū)綜合形變速度場，如圖2所示。由于累計形變點在空間上屬于非均衡分布，導(dǎo)致相干條件差的區(qū)域提取到的PS點數(shù)量較少，

插值效果不佳，主要為和林縣北西側(cè)的集約化農(nóng)墾區(qū)域?？傮w上看，縣域大部分區(qū)域PS點密度較高，反演的年均形變速度場能如實反映區(qū)域地表形變的一般性規(guī)律。研究區(qū)綜合形變速度場分布特征為東南側(cè)隆升、西北側(cè)沉降，構(gòu)造隆升與侵蝕夷平并存。構(gòu)造上主要受NE方向的塊斷構(gòu)造控制（《1∶25萬呼和浩特市幅（K49C004003）區(qū)域地質(zhì)調(diào)查報告》），屬于陰山-燕山板內(nèi)造山帶的西段，鄂爾多斯盆地東緣。本區(qū)自燕山運動以來處于階段性隆升態(tài)勢[24]，全新世至今山體不斷隆升，形成了山前廣為發(fā)育的河流階地和沖洪積平原。其中，新近紀古夷平面的發(fā)現(xiàn)，證實了本區(qū)地表形變的主因是構(gòu)造隆升與侵蝕夷平的交互驅(qū)動。由圖2可知，沉降區(qū)主要分布在黃土丘陵及陸前平原，位于NE向高角度正斷層的北西側(cè)（該斷層分割了山區(qū)與平原），也是該縣的人口聚集區(qū)，沉降成因可能包含多種因素，諸如地質(zhì)凹陷、人類活動、滑坡、泥石流、土壤侵蝕等，需進一步細分形變機制，定向提取土壤侵蝕形變。

2.2 累計形變點篩選

雷達時序差分干涉在縣域內(nèi)共計提取累計形變點1 007 958個，包含不同區(qū)位條件、不同驅(qū)動因素的地表形變。從1 007 958個累計形變點中篩選出土壤侵蝕形變點，需要對全部累計形變點進行再次篩選。本文設(shè)計了多尺度濾波流程對累計形變點進行篩選，包括空間域濾波和時間域濾波。空間域濾波用于消除人類擾動影響。結(jié)合2017土地利用數(shù)據(jù)掩膜人類活動形成的不透水面及工作面；再使用振幅標準差濾波排除茂密森林植被及活躍耕地，使監(jiān)測范圍縮小至天然裸露（半裸露）地表；然后設(shè)置分割閾值，排除無形變或形變極小的穩(wěn)定地表，依據(jù)前人研究設(shè)置穩(wěn)定地表閾值為±1 cm/a[25]?？臻g域濾波結(jié)果如圖3所示。宏觀上，累計形變點集中分布于縣域東南部的山區(qū)丘陵地區(qū)，山區(qū)丘陵是地表形變的多發(fā)易發(fā)區(qū)域，同時相干條件較好，永久散射體密度較高。濾波區(qū)主要為縣鎮(zhèn)村的建成區(qū)、集約化農(nóng)業(yè)區(qū)、林區(qū)、公路鐵路及設(shè)施用地等。

時間域濾波用于分離地表自然形變。不同類型自然形變的物候特征不同，在時間序列上的差異性顯著，具體表現(xiàn)在形變方向、形變區(qū)間、形變幅度、形變趨勢差異等。形變方向分為抬升正向形變和沉降負向形變；形變區(qū)間分為全年形變、季節(jié)形變和突發(fā)形變；形變幅度分為平緩形變和劇烈形變；形變趨勢分為單向漸進形變和雙向波動形變。根據(jù)土壤侵蝕的物候特征，其形變特征表現(xiàn)為全年均有發(fā)生、在雨季驟增的漸近式負向形變。根據(jù)時間序列的形變差異，進行有監(jiān)督的閾值分割，可以初步圈定侵蝕發(fā)生地。再與同步觀測的Sentinel-2A高分光學(xué)遙感影像疊加，根據(jù)貝葉斯理論增加已知信息量[26]。檢驗形變點所處區(qū)域的光譜特征、區(qū)位特征、形態(tài)學(xué)特征和動態(tài)變化軌跡是否與典型侵蝕景觀一致。若否，則返回上一步，迭代更新初始閾值，直至確定最優(yōu)分離閾值。

經(jīng)由時間域濾波及多源數(shù)據(jù)信息符合，確定土壤侵蝕形變的最優(yōu)分離區(qū)間及約束條件，見公式（2）～（6）

式中Sc為年累計形變量（mm），Si為單月形變量（mm），i為月份，n為年份，R為區(qū)域降雨系數(shù)。時間域濾波的基本邏輯是，基于土壤侵蝕的物候特征及先驗知識，對累計形變點的形變方向、形變區(qū)間、形變幅度及形變趨勢加以限定，從而分離出符合侵蝕特征的形變點。經(jīng)由空間域和時間域濾波，從1 007 958個累計形變點中，提取出土壤侵蝕形變點 10 596個，占累計形變點總數(shù)的1.05%。

2.3 區(qū)域土壤侵蝕特征

2.3.1 空間特征

前人研究表明土壤侵蝕與植被[27]、地形[28]、降雨[29]、土壤[30]等因素具有顯著的相關(guān)性。將10 596個土壤侵蝕點與同步觀測計算的動態(tài)數(shù)據(jù)（NDVI分布圖、降雨量插值圖）及靜態(tài)數(shù)據(jù)（坡度圖、土壤可蝕性）進行信息復(fù)合，評估本區(qū)土壤侵蝕的宏觀分布特征（圖4）。其中NDVI由2017年7月的Sentinel-2A數(shù)據(jù)計算歸一化植被指數(shù)得出；降雨量由和林格爾縣氣象局提供的19個氣象觀測站的 2017—2018降雨數(shù)據(jù)，進行克里金插值得到；坡度分布由地理空間數(shù)據(jù)云下載的DEM數(shù)據(jù)計算得出；土壤可蝕性由中科院南京土壤所的土壤數(shù)據(jù)計算K因子。

圖4a可見侵蝕點云多分布于低植被或無植被覆蓋區(qū)域，植被群落以草本植物為主，且較為稀疏，光學(xué)遙感特征近似裸地。在NDVI介于0至0.5區(qū)間內(nèi)（NDVI大于0.5區(qū)段被濾除），土壤侵蝕點主要分布于NDVI零值附近，土壤侵蝕多發(fā)于無植被保護區(qū)域。圖4b可見縣域降雨量分布極不平衡，整體呈東高西低分布，東部山區(qū)的年降雨量是西北部平原的 3倍。本區(qū)土壤侵蝕主要分布于300 mm等降雨量線以下的干旱半干旱區(qū)域。圖4c可知平原與山區(qū)丘陵的地形分界線約為北東向 40°，線性構(gòu)造帶發(fā)育，從地質(zhì)地貌角度屬于土壤侵蝕多發(fā)、易發(fā)地帶。土壤侵蝕點主要分布于坡度 0～35°范圍內(nèi)，平原區(qū)的侵蝕類型主要為風(fēng)蝕，山區(qū)的侵蝕類型主要為水蝕。圖4d可見土壤可蝕性與侵蝕點云一致性較高。侵蝕點云聚集區(qū)多位于土壤可蝕性較高區(qū)域，在土壤可蝕性較低區(qū)域則分布較少且稀疏。總體而言，本區(qū)土壤侵蝕主要分布于山區(qū)平原過渡帶，侵蝕點云約呈北東向40°分布，貫穿縣域南北，在縣域西北、縣域東南也呈局部聚集態(tài)勢。

2.3.2 時間特征

對PS-InSAR觀測結(jié)果進行月度回溯，對全部土壤侵蝕點的累計侵蝕量按月取平均值，同時計算月度累計降雨量平均值。如圖5所示。

研究區(qū)土壤侵蝕總量逐月連續(xù)增加，即便在沒有降雨的12月、1月、2月，土壤侵蝕仍在發(fā)生，說明風(fēng)力是本區(qū)土壤侵蝕的重要營力。此外，在降雨集中的 6—9月，土壤侵蝕量顯著增加，變化梯度增高，表明降雨對研究區(qū)的土壤侵蝕具有重要驅(qū)動作用。研究區(qū)土壤侵蝕的主要營力為風(fēng)力水力混合，侵蝕類型為風(fēng)蝕水蝕混合。平均侵蝕速率為15 mm/a，風(fēng)蝕強度為中度?！锻寥狼治g分類分級標準SL190—2007》將本區(qū)劃分為內(nèi)蒙古高原草原中度風(fēng)蝕水蝕區(qū)，本文監(jiān)測結(jié)果與之一致。

2.4 監(jiān)測精度評估

PS-InSAR方法是一種直接觀測方法，比之光學(xué)遙感的米級分辨率，多時相SAR監(jiān)測地表形變的精度理論上可達厘米級甚至毫米級[8,12]，監(jiān)測量綱提升2個數(shù)量級。然而，受限于初始數(shù)據(jù)精度和研究對象的復(fù)雜程度，實踐中可能難以達到毫米級的測量精度。

InSAR方法的精度驗證多為與其他高精度GPS數(shù)據(jù)進行比對[31]、及實地調(diào)查形變徑跡[32]等方式，本文同時采用上述兩種方式驗證精度。依據(jù)《水土保持工程調(diào)查與勘測標準》（GB/T 51297—2018）規(guī)定，課題組赴研究區(qū)開展了野外調(diào)查驗證，共調(diào)查 103處土壤侵蝕點，采用侵蝕徑跡調(diào)查、高精度GPS測量、實地拍照的方式，驗證所提土壤侵蝕點與侵蝕徑跡的對應(yīng)準確度（如圖6）。

調(diào)查中發(fā)現(xiàn)，絕大多數(shù)調(diào)查點位于侵蝕溝（圖7a）、風(fēng)蝕柱（圖7b）、新鮮紅黏土（圖7c）等典型侵蝕徑跡區(qū)域。經(jīng)統(tǒng)計，土壤侵蝕點的提取準確度為86.4%。此外，收集了和林格爾縣境內(nèi) 142個四等水準點數(shù)據(jù)，依據(jù)最臨近原則對四等水準點及其最近的PS點的絕對高程進行比對，同時與SRTM DEM（30 m）在水準點位置的絕對高程進行比對（圖7d）。由于PS點與四等水準點在空間位置上不重復(fù)，故設(shè)置最大搜索距離為 10 m，相距超過10 m的點對不納入統(tǒng)計。統(tǒng)計后發(fā)現(xiàn)，四等水準點與其最近的PS點的絕對高程值非常接近，且接近程度與距離反相關(guān)。PS點和水準點相距1 m以內(nèi)的絕對高程相對偏差為7～8 cm；相距10 m以內(nèi)的相對偏差介于幾米至十幾米之間。尤其在水準點1 m以內(nèi)的監(jiān)測精度優(yōu)于SRTM DEM兩個數(shù)量級（后者為米級監(jiān)測精度）?？傮w來看，本文方法監(jiān)測土壤侵蝕形變具有可靠性。

3 討 論

國內(nèi)外使用InSAR方法監(jiān)測土壤侵蝕的研究較少，該領(lǐng)域整體研究程度不高，缺乏成熟可供借鑒的研究經(jīng)驗。本文旨在探索并發(fā)展PS-InSAR監(jiān)測土壤侵蝕的理論基礎(chǔ)和技術(shù)框架，并論證技術(shù)框架的可行性和有效性，為后續(xù)更深入的研究提供參考。方法運行結(jié)果表明，監(jiān)測到的土壤侵蝕形變在空間上與植被、土壤、降雨、地形因子相關(guān)，在時間上與降雨季節(jié)相關(guān)，在實地驗證中也發(fā)現(xiàn)了侵蝕徑跡，且與四等水準點的一致性較高，說明 PS-InSAR方法監(jiān)測土壤侵蝕具有可行性和可靠性。從精度驗證結(jié)果來看，所述方法具備厘米級監(jiān)測精度，可以識別中度以上風(fēng)蝕和極強烈水蝕，適用于中國北方干旱半干旱地區(qū)的土壤侵蝕動態(tài)監(jiān)測及分析評估。后續(xù)研究包括：1）采用分辨率更高的InSAR和DEM數(shù)據(jù)源，改進優(yōu)化技術(shù)框架，提高監(jiān)測精度至亞厘米級或毫米級，并開展“三同（同時、同地、同精度）”地面精密水準測量以驗證監(jiān)測精度；2）復(fù)雜地形地貌區(qū)域的土壤侵蝕監(jiān)測及侵蝕模數(shù)反演；3）基于InSAR監(jiān)測的土壤侵蝕模型構(gòu)建或現(xiàn)行水土流失模型修正等。

地表形變場是多重驅(qū)動力下的混合形變效應(yīng)，提取侵蝕形變的基本邏輯是在空間上圈定侵蝕易發(fā)地，在時間上符合侵蝕發(fā)生發(fā)展的物候規(guī)律。本文設(shè)計了多尺度濾波流程，對 1 007 958個累計形變點進行篩選，提取出侵蝕形變點10 596個，占累計形變點總數(shù)的1.05%。需要說明的是，本文的研究假設(shè)是各類形變在水準面上具有二維可分性（比照光學(xué)遙感分類），然而部分區(qū)域可能構(gòu)造運動與土壤侵蝕并存，土壤侵蝕量的一部分被構(gòu)造隆升量所抵消，從而低估實際土壤侵蝕量。本文的處理辦法是，在發(fā)現(xiàn)構(gòu)造隆升趨勢后，返回到提取永久散射體步驟，重新選擇了參考點，對每個參考點進行多次回歸迭代，更新高程誤差改正值，修正了最終形變結(jié)果。因此，在現(xiàn)今構(gòu)造運動活躍區(qū)域開展土壤侵蝕監(jiān)測，還需深入解析形變機制機理，疊加或分離現(xiàn)今構(gòu)造運動修正。

研究區(qū)位于黃土高原北坡，馬蘭期風(fēng)成黃土屬于自重濕陷性黃土，易發(fā)生侵蝕形變[33]。本區(qū)土壤侵蝕的一大特征是與人類活動相關(guān)，侵蝕嚴重區(qū)域的人類活動較為密集，約10%的侵蝕點位于耕地范圍。耕地范圍的侵蝕點主要位于撂荒地、坡耕地及大量井灌引發(fā)的地下水超采區(qū)域。土壤侵蝕對耕地耕層質(zhì)量退化作用表現(xiàn)為土壤性質(zhì)惡化、土壤質(zhì)量劣化、土地生產(chǎn)力衰退3個方面。應(yīng)注意避免陡坡開墾，加大退耕還林還草力度。在縣域西北極干旱地區(qū)，年均降雨量只有150 mm左右，卻是該縣的糧食主產(chǎn)區(qū)，且未見顯著的河流、湖泊等地表水補給源，集約化農(nóng)業(yè)大量依賴地下水，勢必引發(fā)地下水超采問題。地下水超采會降低土壤含水率，導(dǎo)致對土壤墑情敏感的植被群落衰亡，從而使局地生態(tài)環(huán)境逆向演替。建議加強對該區(qū)域地下水位的動態(tài)監(jiān)測，疏導(dǎo)產(chǎn)業(yè)結(jié)構(gòu)向低耗水量行業(yè)轉(zhuǎn)型。

4 結(jié) 論

本文利用土壤侵蝕主要特征是土層厚度剝蝕、且合成孔徑雷達對高程信息敏感的特點，研究PS-InSAR技術(shù)監(jiān)測土壤侵蝕的可行性?；谕寥狼治g的物候規(guī)律對技術(shù)框架進行了針對性的改進。通過多源數(shù)據(jù)融合增加已知信息量，構(gòu)建時空域多尺度濾波方法分離混合形變效應(yīng)，并給出了時間域濾波模型及參考值，使監(jiān)測對象錨定于土壤侵蝕形變。經(jīng)驗證，本文方法獲得了與實際侵蝕發(fā)生地一致的監(jiān)測結(jié)果，土壤侵蝕點提取準確度為86.4%，高程向總體監(jiān)測精度為厘米級。侵蝕點在空間上表現(xiàn)出與植被、降雨、坡度、土壤的強相關(guān)性，在時間上表現(xiàn)為在雨季驟增的非線性負向形變。侵蝕外營力為風(fēng)力水力混合，平均侵蝕速率為15 mm/a。本文方法或可為土壤侵蝕的宏觀面上研究提供參考。