基于ENDVI-SI3特征空間的鹽漬化反演模型及風險評估

張思源， 岳楚， 袁國禮， 袁帥， 龐文強， 李俊

(1.中國地質調查局呼和浩特自然資源綜合調查中心，呼和浩特 010010； 2.中國地質大學(北京)地球科學與資源學院，北京 100083； 3.巴彥淖爾市現代農牧事業發展中心，巴彥淖爾 015000)

0 引言

土壤鹽漬化是氣候、地形、土壤和水文地質等自然條件共同作用下，導致水鹽運動變化，進而產生的生態問題。人類活動中的農業和水利等又對土壤鹽度演化產生重要影響，甚至導致土壤次生鹽漬化[1]。例如蓄水盆地、溝渠和土壤灌溉會引起當地水鹽狀況的變化，加重鹽漬化進而破壞土壤肥力，影響農作物生產，最終威脅糧食庫存。鹽漬化問題對干旱、半干旱地區人口的社會和經濟影響嚴重[2]。及時檢測土壤鹽分，監測和評估土壤鹽分脅迫的嚴重程度，在局部和區域尺度上就顯得非常重要[3]。

測量土壤鹽分的傳統方法是收集現場的土壤樣品并在實驗室分析，以確定其溶質濃度或電導率。然而，由于需要密集采樣才能充分表征該區域的空間變異性，既耗時又昂貴[4]。遙感(remote sensing，RS)技術通過快速、及時、相對廉價和可重復性的數據，對鹽漬化及時監測和反演有巨大優勢。地理信息系統(geographic information system，GIS)提供了在規模、時間、來源和結構方面集成不同性質數據的平臺優勢。通過RS和GIS配合適當的地面工作確定某一特定地區最適當的反演指標，能夠準確高效地提取鹽漬信息，目前已經成為國內外學者進行土壤鹽漬化定量檢測與分析的有效手段[5]。

目前已有許多國內外學者應用波段光譜組合和地表反照率等特征參量建立特征空間，對土壤鹽漬化開展定量研究，直接解譯鹽漬化土壤光譜與依據植被特征間接反演相結合，將鹽漬化信息量反映至最大[6]。研究主要集中于3種鹽漬區： 濱海鹽漬區、極端干旱漠境鹽漬區和半干旱綠洲農田鹽漬區。Thiam等[7]在塞內加爾的薩盧姆河三角洲濱海區通過電導率結合高程、離河距離、歸一化植被指數(normalized difference vegetation index，NDVI)、鹽分指數(salinity index，SI)和土壤調整植被指數，采用多元回歸分析建立了鹽漬化模型； 潘肖燕等[8]和樊彥國等[9]針對黃驊濱海鹽漬區和黃河口濱海三角洲鹽漬區構建NDVI-SI特征空間，對濱海鹽漬區進行空間分析及生態風險定量評價； 哈學萍等[10]基于地表反照率(Albedo)-SI特征空間對處于極端干旱漠境鹽漬區的于田縣進行鹽漬化的定量研究。針對半干旱綠洲農田鹽漬區的研究較多，多集中在新疆各地的綠洲農田。王飛等[11]和陳實等[12]建立基于NDVI-SI特征空間對塔里木南緣于田綠洲和天山北疆綠洲農區土壤鹽漬化進行時空變化分析，能夠顯著區分不同鹽漬化程度地類。在此基礎上，張添佑等[13]通過增加Albedo和土壤濕度(wetness index，WI)等反演指標，與SI和NDVI構建了生態指數對新疆瑪納斯河流域綠洲土壤鹽漬化進行監測與評價，在主成分貢獻中，植被和鹽分特征參量占85%以上，其他特征參量占比過小，造成信息冗余。在以往半干旱綠洲區域研究中，對于指數的選擇多選用NDVI與SI，引入擬合效果更優的指數特征參量將能達到更優擬合效果和地區適宜性。

黃河上游的內蒙古巴彥淖爾市河套灌區作為我國農業生產中重要的商品糧、油基地，是我國三大特大型灌區之一[14]。河套灌區農業生產受限于半干旱區較低的降水量，對灌溉的依賴性強，潛在的土壤次生鹽漬化生態環境問題顯著。對于河套灌區鹽漬化信息提取和定量分析，選取適合鹽漬反演指標的研究較少，僅盧晶等[15]利用SI-MSAVI特征空間對巴彥淖爾市磴口縣進行定量監測和相關性分析。巴彥淖爾市烏拉特前旗位于河套平原東部，是河套灌區的最下游，灌區內90%以上的高鹽高堿農田排水都注入旗內的烏梁素海，次生鹽漬化形勢嚴峻。同時由于烏拉特前旗原生鹽漬土地發育，鹽漬化防治難度大。掌握鹽漬土地空間分布和類型程度，綜合防治土壤鹽漬化已成為其農業發展和土地可持續利用的首要措施[16]。本文以烏拉特前旗為研究區，基于Landsat8 OLI數據和野外實地采樣數據，以引入短波紅外波段的增強型歸一化植被指數(enhanced normalized difference vegetation index，ENDVI)與半干旱區鹽分擬合效果最優的鹽分指數3(salinity index 3，SI3)為特征參量，組成ENDVI-SI3特征空間并構建改進型鹽漬化監測指數(improved salinization monitoring index，ISMI)，在此基礎上建立土壤特征變量與鹽漬化過程的定量關系，為區域尺度的土壤鹽漬化遙感監測和風險評估提供科學的參考依據。

1 研究區概況

烏拉特前旗位于內蒙古自治區巴彥淖爾市，河套平原東部(圖1)。地理坐標： N40°28′～41°16′，E108°11′～109°54′，總面積為7 476 km2。烏拉特前旗地貌類型多樣，旗區內平原、丘陵、山地、沙漠交錯，海拔范圍為938～2 326 m，自南向北烏拉山、白音查干山和查石太山3座山脈之間分別發育明安—大佘太川和小佘太川2個帶狀洼地，黃河在前旗境內自西向東貫穿烏拉山南部的山前洪積扇傾斜平原，即三湖河平原。中國八大淡水湖之一的烏梁素海位居旗區中央，烏梁素海西部為河套平原東緣。

圖1 烏拉特前旗地理位置示意圖Fig.1 Geographical position of Urad Front Banner

烏拉特前旗氣候劃分上屬溫帶大陸性氣候西部季風區，中亞半干旱荒漠區，冬冷夏熱，四季分明。年平均氣溫為5.3 ℃，年日照時數為3 208.1 h，年平均風速為3.4 m/s。年降水量介于200～250 mm之間，多集中于6—9月，年蒸發量為2 217.3 mm，蒸降比大[17]。

在地質構造上，烏拉特前旗所處的河套平原屬于湖相沉積的內陸封閉式斷陷盆地，盆地基底從東南到西北逐漸降低； 上更新統晚期，黃河形成并多次改道，導致盆內湖相沉積層上蓋了黃河沖積層，現代地形雛形基本形成。地面排水受限于南高北低的地形和位于南部的黃河位置，地下徑流受限于黃河流向導致含水層由西向東粒徑逐漸變細，地面和地下排水條件均不利于天然徑流排泄。研究區內主要為通氣透水的輕質土壤，埋深較淺的地下水在強烈的蒸發作用下，構成了有利于潛水蒸發排泄的基礎，這導致大量鹽分富集在淺層土壤中，鹽漬化問題顯著[18]。

2 數據來源與研究方法

2.1 土壤數據采集與化學分析

根據研究區鹽漬地類型差異性與空間均勻性等因素，2019年在野外采集66個土壤樣品，土壤樣品的采樣深度為0～20 cm，為減少單個樣品誤差，在采樣點周圍3～5處多點采集組合。野外采回的樣品在日光下曬干，干燥后的樣品用尼龍篩過篩剔除雜質便于化學分析，具體分析項目包括： 表層土壤含鹽量、容重、pH、有機質、全氮、有效磷、鈣離子、鎂離子、鉀離子、鈉離子、硫酸根、氯離子、碳酸氫根。

2.2 遙感數據處理

2.2.1 數據源與預處理

由于春季影像能反映植被自然生長趨勢，受人為灌溉和耕作影響小，且處于返鹽季節，土壤表層含鹽量高[19]，故遙感數據采用2019年4月30日的2幅Landsat8影像，行列號分別為 128/31 和128/32， 空間分辨率為 30 m，數據均來源于中國科學院計算機網絡信息中心地理空間數據云平臺(http: //www.gscloud.cn)。陸地成像儀(operational land imager，OLI)成像過程中會受到傳感器內部、太陽輻射和大氣影響，所以通過ENVI5.3軟件對原始影像數據進行輻射定標和FLAASH大氣校正等預處理來消減誤差； 將校正后的影像數據進行圖像融合、投影變換與裁剪處理，得到研究區影像，便于土壤鹽漬化的建模反演。

2.2.2 植被指數

土壤表面的植被發育狀況或鹽生植物能夠反映土壤鹽漬化問題。在正常條件下，不健康植被的光合活性較低，導致植被的可見光反射率增加，近紅外反射率降低[20]。在各種受鹽脅迫的植物中，這種模式被廣泛發現[21]，從而有可能從植被反射率間接反演鹽漬化程度。在此基礎上，植被指數常作為評價土壤鹽分的間接指標。常用于反演鹽分的植被指數包括NDVI、修正土壤調節植被指數(modified soil adjusted vegetation index，MSAVI)和增強植被指數(enhanced vegetation index，EVI)等[22]。

但是多數植被指數都僅含有近紅外和可見光中的2個波段信息，信息量較少。故在此基礎上引入短波紅外波段進行補充，構建含短波紅外信息的拓展植被指數，有利于減少波段間的多重共線性，提高土壤鹽漬化建模的精確度[23]。其中ENDVI為 NDVI引入短波紅外2波段得出，其數值大小隨表層含鹽量增加而減少，在植被類型復雜和稀疏程度不一的地區都有較好的反演效果，因此本文選擇ENDVI作為反演土壤鹽漬化的一個指標，其計算公式為：

(1)

(2)

式中NIR，R和SWIR2分別為遙感影像中的近紅外、紅光和短波紅外2波段的反射率值。

2.2.3 鹽分指數

鹽分指數作為直接反演土壤鹽漬化水平的特征指標，在土壤鹽漬化監測中能提供較高精度的反演結果。常見的鹽分指數包括歸一化鹽分指數(normalized difference salinity index，NDSI)、鹽分指數1(salinity index 1,SI1)、鹽分指數2(salinity index 2，SI2)、鹽分指數3(salinity index 3，SI3)、鹽分指數(salinity index，S2)、鹽分指數(salinity index，SI-T)等[24]。

Khan等[25]利用IRS-1B衛星的LISS-II傳感器提出了NDSI和SI，用于評估巴基斯坦的鹽漬化土地退化； Bannari等[26]從EO-1 ALI光譜波段提出了3種不同的鹽分指標SI1，SI2和SI3，用于區分摩洛哥的輕度和中度鹽漬土壤，其中SI3具有最高的相關性； Allbed等[24]發現在植被密度較低和裸地較多的地區， NDSI和SI-T與土壤鹽分的相關性最高。隨后王爽等[27]和張俊華等[28]在此基礎上對國內半干旱綠洲地區如寧夏銀北、新疆渭干河—庫車河三角洲等地的鹽分光譜特征進行研究分析，光譜反演模型中SI3為土壤鹽分反演的最佳指數，SI3與對各項鹽分指標尤其是Cl-和K+均呈極顯著相關性。基于相關性最強的原則，選取SI3作為構建特征空間的變量之一，其計算公式為：

(3)

式中G為遙感影像中綠光波段反射率值。

2.3 ENDVI-SI3特征空間構建

在ENVI5.3平臺下通過Band Math工具運算獲得ENDVI與SI3數據，為消除變量數據的量級差異，需對研究區內ENDVI與SI3進行統計分析得出最值，然后進行標準化處理，公式為：

(4)

(5)

式中：ENDVImax和ENDVImin分別為ENDVI′的最大值和最小值；SI3max和SI3min分別為SI3′的最大值和最小值。

標準化處理后，修正異常值，利用標準化后的ENDVI和SI3數據構建二維散點圖(圖2)，橫坐標為ENDVI，代表地表植被覆蓋度的變化，縱坐標為SI3，代表土壤表層含鹽量的變化。由散點圖和前人研究[29]可得，ENDVI與SI3大致呈負相關，即當土壤表層含鹽量逐漸增加時，地表植被受鹽分脅迫，其生長受到影響，植被覆蓋度降低。

圖2 ENDVI-SI3二維散點圖Fig.2 ENDVI-SI3 two-dimensional scatter plot

2.4 鹽漬化監測模型的構建

根據Verstraete等[30]的研究結果，垂直于特征軌跡曲線l的直線能夠將不同程度的鹽漬化土壤區分開，圖2中任一點到點A(1，0)的距離可以表示土壤表層含鹽量情況，點A為ENDVI為1，SI3為0的理想點，特征空間內的點距離點A越遠，土壤鹽漬化程度越高。基于此建立特征空間土壤鹽漬化監測模型ISMI為：

(6)

3 模型應用及驗證

3.1 ISMI精度驗證

利用ArcGIS軟件的Extract Values to Points功能，計算66個采樣點對應的ISMI，與實測表層土壤含鹽量作相關性分析。結果顯示，ISMI與實測含鹽量相關系數達到了0.82，并且通過0.01(雙側)水平的顯著性檢驗。與此同時將各項指數分別提取并與含鹽量做相關性分析(表1)，發現ISMI與含鹽量的相關系數不僅高于ENDVI和SI3與含鹽量的相關系數，也高于任意一項單獨指數與含鹽量的相關系數，即ISMI在反演鹽分的精確度上明顯優越于單項指數。在單項指數的比較中，ENDVI顯著優于其他植被指數和增強型植被指數，相較于未引入短波紅外的NDVI反演精度也有明顯提升； 在各項鹽分指數中，SI和SI3等指數與含鹽量的相關系數均超過了0.7，反演結果較優，SI3反演結果在鹽分指數中最優。從相關分析結果可以看出，ENDVI與SI3這2個特征參量的優越性高于其他指數，同時ISMI對于土壤鹽漬化的反演精度相較于單獨的NDVI，ENDVI和SI3大幅提升，說明ISMI能夠較好地反映研究區土壤鹽漬化現狀，能夠利用ISMI對研究區的鹽漬土壤進行快速定量監測。

表1 指數相關性分析結果Tab.1 Results of index correlation analysis

為了避免水體與山體對于鹽漬土壤提取的干擾，以中國科學院地理科學與資源研究所資源環境科學數據中心的中國1∶100萬比例尺地貌類型空間分布數據與土地利用數據為原始數據，將其平滑并修正提取出水體、山體和沙丘數據，通過掩模去除后，進一步對平原區不同程度的鹽漬化土壤進行統計分析，根據《中國鹽漬土》中半干旱、干旱區的鹽化分級[31]將鹽漬化程度分為5級(表2)。

表2 基于ISMI的鹽漬土壤分級Tab.2 Classification of saline soil based on ISMI

根據不同程度鹽漬土對應的ISMI進行重分類，得到烏拉特前旗土壤鹽漬化空間分布情況(圖3)，研究區內鹽漬土分布以烏梁素海為界，呈西高東低，不同程度的鹽漬土分布廣泛，各個鄉鎮的鹽漬化程度有較大差異。

圖3 基于ISMI的烏拉特前旗鹽漬土壤空間分布Fig.3 Spatial distribution of saline soil in Urad Front Banner based on ISMI

除含有大面積水體的烏梁素海漁場外，弱中度鹽漬土在各鄉鎮均有大面積分布，兩者占比在40%～65%左右。綜合分析表明，研究區內土壤鹽漬化程度高，ISMI均值為0.25。將水體和山體掩模后的平原地區ISMI均值為0.26，逼近中度鹽漬化的水平，需要重視鹽漬化現狀問題的鄉鎮為西小召鎮、新安農場、西山咀農場和蘇獨侖農場。

3.2 鹽漬化風險評估模型

3.2.1 鹽漬化風險評估因子處理

利用SPSS軟件對土壤分析項目做相關性分析，計算Pearson相關系數并進行顯著性檢驗。篩選出與ISMI存在顯著相關(P<0.05)的含鹽量、有機質、有效磷、鈣離子、鎂離子、鈉離子、硫酸根和氯離子8個變量指標，作為鹽漬化風險評估因子。

灰色關聯分析法為通過比較關聯度大小得出待識別對象對研究對象影響程度的數據分析方法，能夠衡量數據間的同步變化程度[32]。將ISMI作為灰色關聯分析法的母序列，將8個變量指標作為子序列，利用MATLAB程序分別計算出子序列的各項與母序列參數的灰色關聯系數最后得出灰色關聯度，將其標準化后得到評價因子對應的權重(表3)。利用ArcGIS軟件的Interpolation功能對8個評估因子選擇Kriging插值法進行空間插值，Kriging插值法廣泛應用于土壤科學中不規則分布樣點的空間插值[33]。同時將8個因子插值得到的評估中間層柵格，運用自然斷點法將其重分類為10級，通過ArcGIS軟件的Raster Calculator功能將8個評估中間層的重分類柵格計算加權得到研究區土壤鹽漬化風險評估值(salinization risk assessment value，SRAV)，公式為：

表3 鹽漬化風險評估因子相關性篩選及賦權Tab.3 Correlation screening and weighting of salinization risk assessment factors

(7)

式中：γ(xi0,xij)為點i在第j個評估因子上相較于與指標變量xi0的灰色關聯度；wj為第j個評估因子的權重系數；m為評估因子的個數，本文為8個。

3.2.2 鹽漬化風險評估模型精度評價

將土壤鹽漬化風險值SRAV與ISMI通過MATLAB軟件進行函數擬合，得出擬合程度最高的模型為Polynomial多項式逼近的quadratic polynomial二次多項式，Robust選擇最小絕對殘差(least absolute residua，LAR)模型，x為SRAV，y為ISMI，得其函數關系如圖4所示。

圖4 SRAV與ISMI的函數關系Fig.4 Functional relationship between SRAV and ISMI

函數關系式為：

(8)

結果顯示，確定系數R2達到了0.62，均方根誤差為0.06 ，誤差平方和為301.80，從擬合結果來看，基于二者存在較強相關性，擬合效果總體較好，表明基于ISMI構建后續鹽漬化風險模型有較高精度，但在SRAV高值區擬合效果較差，因為SRAV高值時，高值影響因子較多，高鹽漬化風險成因復雜，后續研究可以繼續通過引入更多元類型的鹽漬化風險評估因子來校正SRAV，達到較好的高值擬合效果。

3.2.3 鹽漬化風險空間分布與評估分級

鹽漬化風險評估值空間分布具有明顯的空間差異性(圖5)。對SRAV進行統計分析可得，SRAV范圍為0.11～0.86，平均值為0.35，總體鹽漬化風險較大，以烏梁素海為界呈西高東低的狀態。

圖5 烏拉特前旗土壤鹽漬化風險評估空間分布Fig.5 Spatial distribution of soil salinization risk assessment in Urad Front Banner

參考相關研究分等定級方法[32，34]及研究區實際情況，除去山體、水體和沙丘后，可將研究區土壤鹽漬化風險評估值分為5級(表4)。由風險分級表明，研究區內風險較大及以上的區域占比總和已達54.40%； 而一般風險以下的區域占比僅為45.60%。風險很大區域占比最大，比例為29.43%，總面積達到1 420.06 km2； 風險較大區域面積為803.66 km2，占比16.66%。綜合來看，基于ISMI的鹽漬化風險評估模型能夠有效評估鹽漬化生態問題的風險，能實現低成本高精度的定期鹽漬化風險評估，為農業生產和生態保護提供背景資料。

表4 烏拉特前旗土壤鹽漬化風險評估分級Tab.4 Grading for risk assessment of soil salinization in Urad Front Banner

4 結論與討論

1)由烏拉特前旗改進型鹽漬化監測模型可得，ISMI與實測含鹽量呈顯著相關，相關系數達0.82，在0.01(雙側)水平上通過顯著性檢驗，高于ENDVI與SI3的相關系數(-0.70和0.75)，同時也高于其他指數的相關系數。與之前研究中ENDVI與含鹽量的相關系數也高于NDVI的相關系數(-0.66)，半干旱區中SI3在鹽度指數里反演效果最好的結論吻合。相較于NDVI，ENDVI和SI3等單項指數，ISMI的反演精度顯著提升。

2)ISMI顯示研究區鹽漬化程度較高，平均值為0.25。平原地區的ISMI平均值達到0.26，接近中度鹽漬化水平。鹽土與強度鹽漬土ISMI在[0.41,0.56)和≥0.56的范圍的總面積達到了337.62 km2，占平原總面積的7.03%，嚴重影響農業生產，需積極采取排水、沖洗、改良等脫鹽農業措施，控制用水，實行保墑，井、渠雙排，井、溝雙排，防止土壤次生鹽漬化引發的一系列生態退化及農業減產等問題。

3)鹽漬化風險評估因子權重大小排序為： 鈣離子>有機質>含鹽量>鎂離子>有效磷>氯離子>鈉離子>硫酸根，鹽漬化風險評估因子插值賦權計算得到鹽漬化風險評估值，其值處于0.11～0.86之間，平均值為0.35，分級后研究區內風險較大及以上的區域占比總和已達54.40%，總體鹽漬化風險較大，以烏梁素海為界呈西高東低的狀態，原因是烏梁素海西部土地主要利用類型只有水澆地，且位于河套灌區，灌溉采用礦化度較高的黃河水淹灌且土壤粉粒含量高，水分入滲速度慢，鹽分常年滯留。