基于機器學習和多光譜遙感的銀川平原土壤鹽分預測

魏慧敏， 賈科利， 張 旭， 張俊華

（1.寧夏大學地理科學與規(guī)劃學院，寧夏 銀川 750021；2.寧夏大學生態(tài)環(huán)境學院西北土地退化與生態(tài)恢復國家重點實驗室培育基地，寧夏 銀川 750021）

土壤是人類進行社會生產(chǎn)生活的基本條件，土壤退化牽制著農(nóng)業(yè)的發(fā)展，監(jiān)測土壤性質(zhì)是一個重要過程，對如何可持續(xù)利用起著重要作用。全球六大洲的土壤中均有鹽漬化現(xiàn)象發(fā)生（除南極洲未有調(diào)查數(shù)據(jù)外），并以每年約1.5×104km2的速度增加［1］。鹽分預測指的是將某些指標與其相對應的地表信息間構(gòu)建關(guān)系模型，眾多學者對鹽漬化預測展開了相關(guān)研究［2-5］。徐紅濤等［6］、姜紅等［7］對新疆不同區(qū)域進行土壤鹽漬化預測模型的構(gòu)建，研究表明機器學習的方法使得模型的預測精度高于傳統(tǒng)的回歸法預測模型。馬國林等［8］基于Sentinel-2A影像，使用XGBoost 選擇有效變量后結(jié)合機器學習進行土壤鹽分反演模型訓練，模型預測效果明顯提高，使得鹽漬化制圖精度得到提升。楊厚翔等［9］選取多個影響鹽堿化的因子，借助BP網(wǎng)絡(luò)模型對黑龍江省地區(qū)進行鹽堿化危險度監(jiān)測，整體準確度高于95%。劉全明等［10］利用微波雷達數(shù)據(jù)對鹽分進行反演，所得BP 神經(jīng)網(wǎng)絡(luò)（BPNN）模型能夠在短時間內(nèi)進行大范圍區(qū)域的鹽分監(jiān)測且優(yōu)于多元線性回歸模型和地理加權(quán)回歸模型。楊練兵等［11］在BPNN 模型基礎(chǔ)上使用遺傳算法優(yōu)化輸入?yún)?shù)子集和隱含層神經(jīng)元數(shù)量并對模型的初始權(quán)重進行優(yōu)化，確定了適用于區(qū)域的反演模型，使得反演結(jié)果中土壤鹽分含量的均質(zhì)性最好。BPNN 模型的使用，提升了人們對于土壤鹽漬化變化及預測的精度，豐富了鹽漬化研究在時空方面的意義。

人工神經(jīng)網(wǎng)絡(luò)在預測方面已經(jīng)有了較多的應用成果，章龍管等［12］將貝葉斯網(wǎng)絡(luò)應用于施工風險管理中，結(jié)果證明該方法對于施工風險及發(fā)展趨勢能夠有效預測，可為現(xiàn)場施工防范管理提供理論指導。畢春光等［13］對玉米病害進行預警研究，證明貝葉斯神經(jīng)網(wǎng)絡(luò)（BNN）模型較BP 神經(jīng)模型的精準度提升了5.49%。將BNN模型應用到土壤鹽漬化預測中，探討不同模型對土壤鹽分的擬合效果，本文以銀川平原為研究區(qū)，基于遙感指數(shù)和鹽漬化影響因子，結(jié)合實測土壤鹽分數(shù)據(jù)構(gòu)建預測模型，選出最優(yōu)模型，預測研究區(qū)不同土壤深度的鹽分分布，為銀川平原土壤監(jiān)測和鹽漬化防治提供理論支持。

1 研究區(qū)概況

銀川平原（37.83°～39.38°N，104.28°～107.65°E）西、北、東三面被騰格里沙漠、烏蘭布和沙漠、毛烏素沙漠環(huán)繞，南部與黃土高原接壤。銀川平原屬于溫帶干旱區(qū)，多年平均氣溫5～9 ℃，年平均降水量約200 mm，蒸降比約10:1，年日照時數(shù)2500～3100 h，年太陽輻射5800×106～6100×106J·m-2。銀川平原主要是由沖洪積平原構(gòu)成，地勢西南高東北低。黃河貫穿銀川平原，優(yōu)越的地理位置使得平原能夠充分發(fā)展水利灌溉，灌溉渠系眾多，農(nóng)牧業(yè)發(fā)達［14］，主要作物有水稻、小麥、玉米。灌淤土和草甸土是主要的農(nóng)業(yè)土壤，熟化程度較高。地下水礦化度為0～6 g·L-1，平原地區(qū)地下水埋深較淺，土壤鹽分累積較重。

2 數(shù)據(jù)與方法

2.1 土壤樣本采集與鹽分測定

本文大致以5 km×5 km 的格網(wǎng)設(shè)置采樣點（圖1），于2021 年3 月7—28 日采樣。采樣時每個樣點采用梅花形采樣法，在采樣點30 m×30 m 范圍內(nèi)采集0～20 cm 和20～40 cm 土壤混勻，運用四分法留約500 g土壤裝進密封袋帶回，并用GPS記錄采樣點的坐標信息及周圍環(huán)境信息，共采集197 個樣本點。采集的土壤樣本剔除雜質(zhì)后自然風干，研磨后過2 mm 篩子，以水土比5:1 配置提取液［15］，每個樣品配置3 組作為對照組，每個樣點取3 組數(shù)值的平均值作為樣點最終pH 和電導率，運用電導率法［16］計算土壤樣本的含鹽量。

圖1 銀川平原采樣點分布Fig.1 Distribution of sampling points in Yinchuan Plain

式中：S為土壤含鹽量（g·kg-1）；EC 為測量所得土壤提取液電導率（S·m-1）。

對測量結(jié)果進行分析，刪除異常樣本點后獲得166個樣點供用。將樣本的鹽漬化程度分為5級［17］，樣本的統(tǒng)計性描述如表1 所示。0～20 cm 土壤樣本中共有91 個鹽漬化樣本，約占總樣本的54.82%。20～40 cm土壤樣本共有65個鹽漬化樣本，約占總樣本的39.16%。各層樣本土壤含鹽量總體變異系數(shù)均超過132%，變異性較強，表明樣本離散程度較高，具有普適性。

表1 銀川平原土壤樣本統(tǒng)計Tab.1 Statistics of soil sample in Yinchuan Plain

將土壤樣本按含鹽量由低到高排序，依據(jù)模型訓練集和結(jié)果驗證集2:1 的比例等間隔選擇110 個樣本進行模型訓練，56 個樣本作為驗證集用于模型驗證［18］。

2.2 遙感影像獲取與處理

研究選用的Landsat 8 OLI遙感影像來源于美國地質(zhì)調(diào)查局官網(wǎng)（https://earthexplorer.usgs.gov/），行列號為129/033 和129/034，空間分辨率為30 m，成像時間為2021年3月17日，與采樣時間同期。利用ENVI 5.3對所選的影像進行大氣校正等預處理流程后，進行各類鹽分指數(shù)的計算和提?。?9］，如表2所示。

表2 鹽分指數(shù)計算公式Tab.2 Calculation formulas of salt indices

2.3 土壤鹽漬化預測模型

2.3.1 鹽漬化影響因子組參數(shù)選擇依據(jù)已有學者對土壤鹽漬化成因條件的相關(guān)研究［26-27］，結(jié)合銀川平原實際情況，選擇8 個環(huán)境變量和3 個人為活動變量（表3）。其中，土壤pH 和含水率結(jié)果由采樣點數(shù)據(jù)計算后進行Kriging 空間插值而得，設(shè)置輸出的分辨率為30 m；地下水數(shù)據(jù)來源參考文獻［17］，將數(shù)據(jù)處理為30 m 分辨率的柵格數(shù)據(jù)；銀川平原數(shù)字高程模型從地理空間數(shù)據(jù)云平臺（http://www.gs-cloud.cn/）獲取，空間分辨率為30 m；地表溫度、增強植被指數(shù)和水體指數(shù)通過Landsat 8 OLI 影像在ENVI 5.3 軟件中計算獲得；土地利用數(shù)據(jù)從中國科學院資源環(huán)境科學與數(shù)據(jù)中心（http://www.resdc.cn/）獲取，強度指數(shù)采用莊大方［28］利用分級方法，土地利用強度依據(jù)不同的土地類型分別確定為不同的值，未利用地、草地、林地、耕地、居民用地、水體分別賦值為1、2、2、3、4、1。

表3 土壤鹽漬化影響因子Tab.3 Influence factors of soil salinization

2.3.2 數(shù)據(jù)處理為消除數(shù)據(jù)獲取過程中人為或者自然因素導致數(shù)據(jù)不標準的問題，先將數(shù)據(jù)進行歸一化處理以消除不同指標之間維度的影響。歸一化公式為：

式中：Amax為樣本數(shù)據(jù)的最大值；Amin為樣本數(shù)據(jù)的最小值；A為預警因子數(shù)據(jù)集；A*為歸一化后的數(shù)據(jù)。

2.3.3 預測模型構(gòu)建

（1）支持向量機（SVM）

SVM模型基于結(jié)構(gòu)風險最小原理，依靠有限的樣本來檢索全局最優(yōu)解，擁有對未知點有較好的泛化效果的優(yōu)點［29］，規(guī)避小樣本學習和局部極值出現(xiàn)的缺點。本文選用RBF（Radial basis function）為SVM 中核函數(shù)類型，懲罰參數(shù)（c）與核參量（g）由樣本進行測試計算得到。

（2）BP神經(jīng)網(wǎng)絡(luò)（BPNN）

BPNN 模型基于反向傳播誤差的方式對數(shù)據(jù)集進行訓練，達到誤差最小化的目的［30］。該模型擁有較強的自適應和自學習能力，可精確逼近任意的非線性關(guān)系［31］。本文基于MATLAB 2020軟件，將鹽分指數(shù)和影響因子作為模型的輸入數(shù)據(jù)，實測鹽分數(shù)據(jù)作為輸出數(shù)據(jù)，創(chuàng)建研究區(qū)土壤鹽分的BPNN 預測模型。經(jīng)多次訓練，在參數(shù)一致的情況下，將模型的輸出值和實際值進行擬合，以誤差作為選擇隱含層節(jié)點的標準，確定最佳訓練效果的隱含層節(jié)點。模型訓練目標的最小誤差設(shè)置為0.00001，訓練次數(shù)設(shè)置為1000次。

（3）貝葉斯神經(jīng)網(wǎng)絡(luò)（BNN）

BNN 模型通過為神經(jīng)網(wǎng)絡(luò)的權(quán)重分配引入隨機性進行正則化，也相當于任意權(quán)重的多種神經(jīng)網(wǎng)絡(luò)的預測［32-33］。與傳統(tǒng)BPNN 模型不同的是貝葉斯估計得到的是節(jié)點權(quán)重的后驗分布概率P( )w|D，并不是固定的值，因此可以為神經(jīng)網(wǎng)絡(luò)的預測添加隨機性。由給定數(shù)據(jù)輸出預測值的分布，公式為：

式中：給定數(shù)據(jù)、預測值分別為m和n；w為每個節(jié)點的權(quán)重；D為訓練集。下同。

所得到的預測值y是P( )w|D中的所有可能的預測值的期望。依據(jù)貝葉斯理論，樣本數(shù)據(jù)決定P(D)。

根據(jù)參考文獻［34］，獲得一個高斯先驗的P(w)：

式中：wj為網(wǎng)絡(luò)的第j個權(quán)重為先驗的第一個和第二個混合成分的方差，且σ1＞σ2；N()為wj以均值為0評價的高斯密度。

依據(jù)所給的先驗分布，運用貝葉斯定理得到后驗概率，若某參數(shù)能夠使后驗概率達到最大化，則該參數(shù)為模型的最佳參數(shù)。BNN 模型能夠克服過擬合及出現(xiàn)局部最小化的現(xiàn)象。

2.3.4 模型評價指標為量化土壤鹽分預測模型的預測效果，本文選擇決定系數(shù)（R2）和均方根誤差（RMSE）2 個常用指標來進行綜合評價。R2的值越高，RMSE的值越低，說明模型擬合效果越好。

3 結(jié)果與分析

3.1 銀川平原土壤鹽漬化空間分布特征分析

使用ArcGIS的空間分析工具對166個采樣點的土壤含鹽量進行插值，根據(jù)鹽漬化程度劃分為5 個等級，得到銀川平原不同深度土壤含鹽量的空間分布特征（圖2）。由圖2 可見，0～20 cm 重度鹽漬化土壤及鹽土主要分布在銀川平原北部地區(qū)，南部地區(qū)鹽漬化程度較輕。20～40 cm 主要以非鹽漬化和輕度鹽漬化土壤為主，銀川平原西北部土壤鹽漬化現(xiàn)象嚴重，部分區(qū)域土壤含鹽量大于6 g·kg-1。整體上研究區(qū)0～20 cm和20～40 cm土壤鹽漬化呈現(xiàn)出南低北高的趨勢，重度鹽漬化土壤和鹽土主要分布在石嘴山地區(qū)。

圖2 土壤鹽分空間分布Fig.2 Spatial distribution of soil salinity

3.2 不同預測模型對比分析

3.2.1 相關(guān)性分析對鹽漬化影響因子、鹽分指數(shù)與實測土壤鹽分進行相關(guān)性分析，所得結(jié)果如表4所示。由表4 可見，影響因子中，0～20 cm 土壤除土壤pH、土壤含水率、人口密度外，其余因子與含鹽量具有較強的顯著性；20～40 cm 土壤含鹽量與土地利用強度、增強植被指數(shù)、水體指數(shù)、地表溫度、地下水埋深、地下水礦化度在0.01 水平上具有極強的顯著性。鹽分指數(shù)中，0～20 cm 土壤含鹽量與S1、S2、S3、S6、SI2、NDSI、SI-T 的相關(guān)性均通過0.05 檢驗。20～40 cm 土壤含鹽量與S1、S2、S4、S5、SI3、NDSI、SI-T 的相關(guān)性均通過0.05 檢驗。故0～20 cm 土壤選擇通過0.01 顯著性水平檢驗的影響因子和鹽分指數(shù)分別作為輸入層參數(shù)進行建模，20～40 cm 土壤選擇通過0.05 顯著性水平檢驗的影響因子和鹽分指數(shù)分別作為輸入層參數(shù)進行建模。

表4 模型參數(shù)與土壤鹽分的相關(guān)性Tab.4 Correlation between model parameters and soil salinity

3.2.2 模型構(gòu)建與對比分析選擇SVM、BPNN 和BNN 模型，以土壤含鹽量為輸出層，相關(guān)性強的影響因子和鹽分指數(shù)為輸入層數(shù)據(jù)，對0～20 cm和20～40 cm 土壤鹽分含量進行模擬，結(jié)果如表5所示。以驗證集樣本對模型進行驗證，結(jié)果如圖3。由表5可見，0～20 cm土壤影響因子變量組的3個模型R2均大于0.76，其中BNN 模型建模效果最好，BPNN 模型次之，SVM 模型最差。由驗證結(jié)果表明，BNN 模型的驗證效果與SVM 模型和BPNN 模型相比結(jié)果最佳，其模擬值與實測值的R2值最大（圖3）。鹽分指數(shù)組的R2整體較低，但均通過0.01 顯著性檢驗。其中，該分組中BNN 模型的驗證集RSME 為最低，BPNN模型擬合誤差最大。對比6 個模型的驗證效果，發(fā)現(xiàn)基于影響因子組建立的BNN模型是0～20 cm土壤鹽分預測最佳模型，故選擇此模型對銀川平原0～20 cm土壤鹽分進行預測。

20～40 cm 土壤鹽分預測模型中（表5），影響因子組模型驗證集的RMSE 均大于2，BNN 模型的R2大于SVM 和BPNN 模型，BNN 模型的RMSE 小于SVM 和BPNN 模型，因此BNN 模型檢驗結(jié)果效果優(yōu)于SVM 和BPNN 模型（圖3）。鹽分指數(shù)組驗證集的R2相比于影響因子組有明顯提升，鹽分指數(shù)組BNN模型是3個模型中誤差最小、擬合程度最佳的模型，可以將該模型用于研究區(qū)20～40 cm土壤鹽分預測。

圖3 土壤鹽分與模型預測值的散點圖Fig.3 Scatter plot of soil salinity and model predictions value

表5 基于不同變量組的機器學習模型Tab.5 Machine learning models based on different variable groups

3.3 銀川平原土壤鹽分預測

利用經(jīng)過訓練的BNN 模型對整個銀川平原0～20 cm、20～40 cm 土壤進行鹽分預測，而后按照鹽漬化程度對結(jié)果劃分等級，得到銀川平原土壤鹽分預測等級圖（圖4），由圖4 可知，0～20 cm 土壤中鹽土主要分布在銀川平原北部的大武口區(qū)和平羅縣，重度鹽漬化和中度鹽漬化土壤主要分布在鹽土四周及平原東部地區(qū)，非鹽漬化和輕度鹽漬化土壤主要分布在平原南部地區(qū)。0～20 cm 土壤主要是以非鹽漬化土壤為主，鹽土分布于大武口區(qū)及平羅縣西部地區(qū)，重度鹽漬化土壤零星分布在銀川平原內(nèi)。根據(jù)像元計算得到，0～20 cm 土壤主要是以非鹽漬化土壤為主，占銀川平原總面積的33.29%，其次為輕度鹽漬化土壤，面積占比為31.86%，中度鹽漬化和重度鹽漬化土壤面積占比分別為23.26%和8.07%，鹽土面積占比（3.52%）不足5%。20～40 cm 土壤中非鹽漬化土壤和輕度鹽漬化土壤面積占比分別為50.54%和30.07%，中度鹽漬化土壤面積占比為12.35%，重度鹽漬化土壤和鹽土面積占比（分別為4.89%、2.15%）均小于5%。

圖4 土壤含鹽量預測等級分布Fig.4 Grade distribution of soil salinity inversion

為進一步驗證預測結(jié)果與實驗室結(jié)果的一致性，將預測結(jié)果與實測鹽分插值結(jié)果進行柵格對比，統(tǒng)計相同等級像元數(shù)據(jù)（表6）。由表6 可知，0～20 cm 和20～40 cm 土壤鹽漬化預測所得的鹽漬化程度與實測含鹽量插值相同等級像元個數(shù)的相同率均大于73%，表明預測模型的結(jié)果與實測數(shù)據(jù)間具有較好的一致性，所選模型可達到檢驗要求并用于研究區(qū)的土壤鹽漬化預測。

表6 實測插值與預測像元統(tǒng)計Tab.6 Measured interpolation and predictive pixel statistics

4 討論

使用3種機器學習的方法建立了對于不同深度的土壤鹽分預測模型，發(fā)現(xiàn)不同深度的土壤對同種變量在相同模型下的建模效果差異較大，對所有結(jié)果進行精度評價后，選擇適用不同深度土壤的最佳模型進行土壤鹽分預測。

本文對多個影響因子與實測值進行Pearson 相關(guān)性分析，發(fā)現(xiàn)敏感因子為地下水礦化度、地下水埋深、土地利用強度等多個因子，與楊思存等［35］、蘇春利等［36］所得結(jié)論相符，后續(xù)可用于挖掘因素與土壤鹽分的隱含信息。0～20 cm 土壤鹽分與各影響因子的相關(guān)性明顯高于20～40 cm 土壤，這與劉繼龍等［37］研究結(jié)果一致，可能是土壤鹽分受到外界影響程度隨深度不同而改變，人們所進行的農(nóng)業(yè)活動主要是在0～20 cm 土壤，對20～40 cm 土壤的干預性較小。

邊慧芹等［38］、何寶忠等［39］對于不同地區(qū)的鹽漬化進行研究，表明鹽分指數(shù)在監(jiān)測應用方面具有一定的優(yōu)勢，選擇多種鹽分指數(shù)與實測鹽分進行相關(guān)性分析，選取相關(guān)性較高的鹽分指數(shù)作為模型參數(shù)進行預測模型訓練。本文采樣時正值春季，氣候干燥，蒸發(fā)強烈，土壤鹽分隨水分的蒸發(fā)而在地表集聚，使0～20 cm土壤含鹽量高于20～40 cm土壤，造成鹽分指數(shù)與兩層土壤鹽分的相關(guān)性差異。土壤是一個復雜的綜合體，常規(guī)的線性回歸模型難以達到較好的精度，機器學習能夠較好地克服這個難題，許多研究者將機器學習用于預測鹽分，所得模型精度高于線性模型［5,40-41］，BNN 算法可運用較少的數(shù)據(jù)得到穩(wěn)定的模型和各層參數(shù)的分布，達到較好的泛化目的，可為鹽分預測提供新思路。本文基于2 個變量組，對3種不同機器學習方法進行建模比較，發(fā)現(xiàn)BNN 模型精度優(yōu)于BPNN 模型和SVM 模型。且0～20 cm 土壤鹽分預測整體效果優(yōu)于20～40 cm 深度的模型效果，與楊寧等［42］所得結(jié)論基本一致。

本文通過對比不同深度和不同輸入?yún)?shù)下多個模型的預測效果，選擇最佳模型進行預測。同一地區(qū)不同季節(jié)不同年份的環(huán)境差異較大，春季銀川平原土壤中的鹽分不斷向地表層匯聚，銀川平原主要為耕地，農(nóng)民在耕種前會對農(nóng)田進行浸灌以此來降低農(nóng)田土壤鹽漬化程度，后續(xù)研究應該考慮不同季節(jié)的環(huán)境變化，選擇多年多季節(jié)數(shù)據(jù)進行預測模型的訓練及驗證。兩組輸入?yún)?shù)對土壤鹽漬化程度的預測精度不同，在下一步工作中將對影響因子與鹽分指數(shù)進行組合參數(shù)，進行區(qū)域土壤鹽分的預測模型訓練選優(yōu)，提高模型的預測精度。另一方面，不同地區(qū)土壤鹽漬化驅(qū)動機制不同，本文所建立的預測模型能否運用到其他地區(qū)還需進行下一步的研究和驗證。

5 結(jié)論

以銀川平原為研究區(qū)，將鹽漬化影響因子和鹽分指數(shù)作為輸入?yún)?shù)構(gòu)建鹽分預測模型，對土壤0～20 cm及20～40 cm土壤進行預測，得到以下結(jié)論：

（1）通過土壤鹽分預測模型對比發(fā)現(xiàn)，BNN 模型進行土壤鹽分預測效果最優(yōu)，影響因子對0～20 cm土壤鹽分預測的影響較大；鹽分指數(shù)對20～40 cm土壤鹽分預測的影響較大。對比分析所選算法的不同變量建模和驗證效果，表明BNN 模型是效果最好的預測模型，神經(jīng)網(wǎng)絡(luò)的引入對模型的訓練有一定的優(yōu)勢。

（2）根據(jù)預測結(jié)果可知，銀川平原0～20 cm 土壤主要是以非鹽漬化和輕度鹽漬化為主，面積占比分別為33.29%、31.86%，重度鹽漬化及鹽土面積占比為11.59%。20～40 cm 土壤中非鹽漬化和輕度鹽漬化土壤面積占比分別為50.54%和30.07%，重度鹽漬化及鹽土面積占比為7.04%。20～40 cm 土壤鹽漬化程度較0～20 cm土壤鹽漬化輕。