基于PS-InSAR技術(shù)監(jiān)測(cè)土壤侵蝕可行性研究

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

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

0 引 言

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

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

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

目前，國(guó)內(nèi)外少見(jiàn)PS-InSAR監(jiān)測(cè)土壤侵蝕的研究，有學(xué)者指出土壤侵蝕造成了上海市的地面沉降[2]，也有學(xué)者注意到浙江四明山自然地表在雨季發(fā)生沉降，推測(cè)為土壤侵蝕所致[18]；國(guó)外在伊朗西部[19]、意大利中部[20]有應(yīng)用InSAR技術(shù)監(jiān)測(cè)土壤侵蝕的研究案例。本研究擬利用哨兵一號(hào)衛(wèi)星平臺(tái)數(shù)據(jù)（Sentinel-1A），基于改進(jìn)后的PS-InSAR方法，對(duì)內(nèi)蒙古和林格爾縣進(jìn)行土壤侵蝕的面上觀測(cè)，結(jié)合多源數(shù)據(jù)對(duì)監(jiān)測(cè)結(jié)果進(jìn)行評(píng)估和驗(yàn)證，檢驗(yàn)PS-InSAR技術(shù)監(jiān)測(cè)土壤侵蝕的可行性，為該方法在土壤侵蝕監(jiān)測(cè)領(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)蒙古高原和黃土高原的過(guò)渡地帶，東南部山區(qū)、黃土丘陵區(qū)的面積占縣域面積的77.7%；西北部屬土默川平原，面積占縣域面積的22.3%。研究區(qū)超過(guò)一半?yún)^(qū)域?qū)儆邳S土高原。黃土高原是世界上土壤侵蝕最嚴(yán)重區(qū)域，每年入黃泥沙量達(dá) 16億 t，且每年由于人類活動(dòng)新增水土流失面積900～1 100 km2。自大規(guī)模農(nóng)耕活動(dòng)開(kāi)始，人類活動(dòng)就已經(jīng)與水力、風(fēng)力并列，成為土壤侵蝕的主要營(yíng)力[1]。根據(jù)《土壤侵蝕分類分級(jí)標(biāo)準(zhǔn)SL190—2007》分類體系，研究區(qū)既屬于西北黃土高原水力侵蝕區(qū)，又屬于內(nèi)蒙古高原草原中度風(fēng)蝕水蝕區(qū)，水力侵蝕與風(fēng)力侵蝕并存。

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

監(jiān)測(cè)數(shù)據(jù)為2017至2018年的17期Sentinel-1A雷達(dá)數(shù)據(jù)，每期1景，用于監(jiān)測(cè)土壤侵蝕。Sentinel-1A的軌道類型為近圓形的太陽(yáng)同步軌道，其軌道高度為693 km，傾角為98.18°。衛(wèi)星重訪周期為12 d，內(nèi)有175個(gè)軌道。為定向監(jiān)測(cè)土層厚度變化，選用干涉寬模式（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波段屬于長(zhǎng)波，可以穿透稀疏植被直接作用于地表。SLC數(shù)據(jù)處理后的距離向分辨率為2.3 m，方位向分辨率為13.9 m。

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

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

1.3 研究方法

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

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

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

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

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

1.3.2 差分干涉處理

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

式中φtopo為地形相位；φdef為形變相位；φatm為大氣延遲相位；φf(shuō)lat為平地相位；φnoise為噪聲相位。

1.3.3 提取永久散射體

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

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

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

相位解纏：相位解纏后可得每個(gè)PS點(diǎn)的絕對(duì)形變速率和DEM修正量，然后可得差分相位的形變相位和殘余地形相位。

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

1.3.4 多尺度濾波

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

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

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

2 結(jié)果與分析

2.1 形變速度場(chǎng)

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

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

2.2 累計(jì)形變點(diǎn)篩選

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

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

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

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

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

2.3.1 空間特征

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

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

2.3.2 時(shí)間特征

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

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

2.4 監(jiān)測(cè)精度評(píng)估

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

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

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

3 討 論

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

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

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

4 結(jié) 論

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