基于全子集-分位數回歸的土壤含鹽量反演研究

張智韜 韓 佳 王新濤 陳皓銳 魏廣飛 姚志華

(1.西北農林科技大學旱區農業水土工程教育部重點實驗室， 陜西楊凌 712100； 2.西北農林科技大學水利與建筑工程學院， 陜西楊凌 712100； 3.中國水利水電科學研究院水利研究所， 北京 100038)

0 引言

土壤鹽漬化程度是影響灌區土地資源經濟效益和作物產量的關鍵因素，及時、準確地監測鹽漬化動態是進行科學防治的重要前提。與傳統的土壤鹽漬化監測方法相比，衛星遙感具有無創、動態、綜合、快速、宏觀等優點，已成為現階段大面積監測土壤鹽漬化的重要研究方向。

國內外學者在利用衛星遙感監測土壤鹽漬化程度方面取得了一定的進展。目前，常用方法主要有遙感圖像分類及解譯法[1]、光譜指數法[2-3]、特征空間分析法[4]。其中，光譜指數法是一種便捷高效的方法。陳紅艷等[2]在傳統的植被指數基礎上加入短波紅外波段，明顯提高了遙感反演土壤鹽漬化的精度。EL HARTI等[3]在鹽分指數(SI)的基礎上加入藍波段，構建了新的鹽分指數OLI-SR，提高了Morocco灌區鹽漬化反演的精度。但是，光譜指數易受到下墊面條件和大氣狀況等多種因素的影響，具有明顯的地域性和時效性。因此，篩選敏感光譜指數對于土壤含鹽量監測具有重要意義。目前，常用的變量篩選方法有灰度關聯法[5]、嶺回歸[6]、LASSO回歸[7]、變量投影重要性分析[8]等，這些篩選方法均為局部最優篩選。全子集篩選法列舉全部可能組合方式，建立全局最優模型，以包含最少自由變量的模型解釋因變量，進而消除共線性的影響。李長春等[9]在估測葉面積指數時發現，采用全子集篩選法建立的回歸模型優于PLSR、SVM和RF模型。馮克偉[10]經過全子集篩選，對野生二粒小麥和硬粒小麥的耐鹽指數分別擬合，預測決定系數達到0.95以上，進一步說明全子集篩選法對于篩選光譜指數具有一定的可行性。

在以往的土壤鹽漬化監測中，大多采用經典線性回歸模型[11-12](如最小二乘法和偏最小二乘法)和機器學習模型[13-14]。經典線性回歸對隨機誤差的分布特征要求嚴格，對于一些實際問題，難以得到無偏、有效的參數估計值[15]。人工神經網絡模型和支持向量機模型是常用的機器學習算法。人工神經網絡模型基于梯度學習方法，在層間傳播的非參數非線性模型，因其具有較強的自學習和自適應能力、非線性映射能力和容錯能力[16]，而被廣泛用于遙感反演中。支持向量機模型主要用于對高維度、小樣本的分類和回歸，將有限樣本信息在模型的復雜性和學習能力之間尋求最佳折中，在保證泛化能力的前提下達到最優學習效果[17]。與機器學習方法相比，分位數回歸模型所需參數少，訓練時間短，不易產生過擬合現象。分位數回歸可以在任意一個分位點準確地描述自變量對于因變量的變化范圍以及條件分布形狀的影響，適合處理波動性大、異質性強的數據[18]，廣泛應用于經濟、醫學等領域，但在監測鹽漬化方面的研究還很少[19]。

本文利用河套灌區解放閘灌域GF-1衛星遙感影像數據，計算多種光譜指數，將通過全變量及全子集篩選法得到的最優光譜指數組合方式作為自變量，分別構建不同深度下的人工神經網絡、支持向量機和分位數回歸3種模型，進而對比分析得到土壤含鹽量反演的最優模型，以期提高植被覆蓋條件下的土壤含鹽量反演精度。

1 材料與方法

1.1 研究區概況

研究區設在內蒙古自治區河套灌區解放閘灌域(40°30′～41°17′N，106°33′～107°31′E)，是河套灌區的第二大灌域，位于河套灌區西部，西接一干灌域，東鄰永濟灌域，灌溉面積1 420 km2，灌溉水主要來源于過境的黃河水。主要農作物包括玉米、小麥、葵花。該灌域地處干旱半干旱氣候區，年均降水量約為140 mm，年均蒸發量約為2 000 mm[20]。研究區內平均坡度為0.2%，地勢平坦使得排水系統的排水能力降低，地下水埋深較淺。加之降雨少、蒸發大、高灌低排，使得研究區內土壤鹽漬化問題較為嚴重。

1.2 遙感數據獲取及預處理

以國產高分一號衛星影像(GF-1 WFV相機)為數據源，通過中國資源衛星應用中心進行下載(www.cresda.com)。衛星影像的成像時間為2018年7月27日，與實測數據日期基本同步，高分一號衛星數據的重訪周期為4 d，空間分辨率為16 m，包括4個波段，分別為藍(0.45～0.52 μm)、綠(0.52～0.59 μm)、紅(0.63～0.69 μm)和近紅外波段(0.77～0.89 μm)。對下載的影像經過幾何精校正、輻射校正、大氣校正、鑲嵌、裁剪預處理過程。并對影像進行低通濾波處理[21]，保存圖像的低頻部分，使圖像平滑，在一定程度上消除了高頻隨機噪聲，提高數據的信噪比。

1.3 土壤樣本采集及化學分析

采樣點的選定考慮到灌域內鹽分分布、植被覆蓋種類(包括葵花35個、玉米22個、小麥8個、果蔬8個和荒地7個)，并顧及點位分布的均勻性，本次試驗在解放閘灌域共設置了80個不同土壤鹽漬化程度的采樣點，采樣點分布情況見圖1。采樣時間為2018年7月12—16日，根據該時期灌區內主要作物的根系活動層所在深度，每個采樣點按照0～20 cm、20～40 cm、40～60 cm分層采樣，并通過手持GPS記錄采樣點位置信息及周圍環境信息。采用五點法采樣，采樣單元為16 m×16 m，每個采樣點取5個土樣，編號后帶回實驗室。

將野外收集的土樣干燥研磨處理后，配置土水比為1∶5的土壤溶液，經攪拌、靜置、沉淀、過濾后，采用電導率儀(DDS-307A型，上海佑科儀器公司)測定土壤溶液電導率EC1∶5，對每個采樣點的5個土樣電導率取平均值作為該樣本樣點處的電導率，并通過經驗公式計算土壤含鹽量S(S=(0.288 2EC1∶5+0.018 3)×100%，剔除圖像中小麥收割后的8個樣本，剩余72個土壤樣本用于本次植被覆蓋條件下的灌區土壤含鹽量反演。將測得的土壤含鹽量從大到小排序，每隔2個樣本取出1個作為驗證集樣本。可保證建模樣本和驗證樣本范圍一致且分布均勻。土壤含鹽量數據的統計特征如表1所示。土壤鹽漬化等級劃分參照文獻[22]。表1中0～40 cm土壤含鹽量為0～20 cm和20～40 cm深度實測含鹽量的平均值；0～60 cm土壤含鹽量為0～20 cm、20～40 cm、40～60 cm深度實測含鹽量的平均值。

圖1 研究區地理示意圖及采樣點分布Fig.1 Schematic of study area and distribution of sampling points

表1 采樣點的土壤含鹽量統計Tab.1 Summary statistics of soil salinity of sampling points

1.4 光譜指數計算及全子集篩選

高分一號遙感影像自身的光譜分辨率不高，對于植被覆蓋條件下的土壤鹽漬化反演模型，僅用多個波段的反射率對比分析提取土壤鹽漬化信息有明顯的局限性，本研究利用高分一號衛星數據的4個波段反射率和12個光譜指數(包括5個土壤指數、3個鹽分指數和4個植被指數，見表2)建立遙感圖像與土壤含鹽量的定量關系。

表2 光譜指數Tab.2 Spectral index summary

注：RB、RG、RR、和RNIR分別為高分一號衛星數據對應4個波段的光譜反射率。

全子集篩選是基于不同自變量的所有可能的組合方式，對縮減后的變量組合通過最小二乘法進行擬合，在所有可能的模型中選擇一個最優模型。該方法簡單直觀，適合自變量較少的條件下使用。其主要計算步驟如下：記K為子變量數目(K=1,2，…，P)；擬合1～P個預測變量的模型；在1～P個模型中，根據“調整后R2最大”準則來選擇P個最優模型；根據驗證集調整后決定系數(Coefficient of determination，R2)、均方根誤差(Root mean squared error，RMSE)從P個模型中選擇一個最佳自變量組合。全子集篩選的模型建立和預測通過Matlab R2017b完成。

1.5 模型建立

1.5.1分位數回歸模型

分位數回歸研究自變量與因變量的條件分位數之間的關系，通過最小化離差絕對值的加權和，依據因變量的條件分位數對自變量進行回歸，進而得到所有分位數下的回歸模型[27]。相較于傳統的回歸方法，其優點主要表現在：①無需對模型中變量進行隨機擾動或進行正態變化，在非正態的干擾下，分位數模型精度更高。②模型中異常點對于模型整體精度影響較小，模型穩定性強。③可以給出任意分位點的參數估計，解釋自變量對不同水平因變量的影響[28]。分位數回歸模型的建立與預測在SPSS 23.0軟件中完成。

對于因變量Y的一個隨機樣本{y1,y2,…,yn}，通常τ分位數的樣本分位數線性回歸要求滿足

其中，yi為因變量，β為自變量的回歸系數，β(τ)為特定分位數τ的最終回歸系統，x′i為自變量。

對于任意的0<τ<1，求解參數估計值的公式為

(1)

求得參數x′iβ(τ)即為唯一的第τ回歸分位數。

1.5.2人工神經網絡模型

人工神經網絡模型從信息處理的角度對人腦的神經元進行抽象，建立某種簡單模型，按不同的連接方式組合不同的網絡。該模型拓撲結構由輸入層、隱含層、輸出層構成。經大量試算，本文輸入層為不同深度下全變量及經過全子集篩選后的敏感光譜指數組合，輸出層為土壤含鹽量，隱含層中的網絡層數為1。其中隱含層采用雙曲正切激活函數，輸出層采用恒等激活函數，可以在一定程度上防止梯度消失且結果可靠；訓練擬合目標誤差擬定為0.1%；培訓類型采用適用于小數據集且精度高的批處理，并選用對應的標度共軛梯度估計突出權重進行網絡訓練；為了消除量綱、數量級等對數據分析造成的影響，將光譜指數標準化處理。本文采用SPSS 23.0軟件中的多層感知器神經網絡進行模型的建立與預測。

1.5.3支持向量機模型

支持向量機模型是在結構風險最小化原理的基礎上，通過一個非線性變換將數據變換到一個高維特征空間，并在該特征空間建立線性模型來擬合回歸函數，很大程度上克服了“離散值多”和“過學習”等問題[29]。徑向基核(RBF)是應用最廣泛的核函數。懲罰參數C和RBF核參數δ是影響SVM的主要參數，其中參數C直接影響模型的穩定性，避免模型在學習和推廣過程中產生欠學習和過學習問題，決定了適應誤差的最小化和平滑程度；參數δ反映了支持向量之間的相關程度，直接影響支持向量之間聯系的松弛度，避免產生欠擬合和過擬合問題，決定模型預測的推廣能力和泛化性[30]。本文核函數采用RBF核函數進行計算，利用R3.5.1軟件進行支持向量機模型的建立與預測。

1.6 統計分析

本文通過決定系數R2、均方根誤差RMSE、赤池信息準則AIC[31](Akaike information criterion)和施瓦茨信息準則SIC(Schwarz information criterion)來綜合評價全子集篩選的效果。通過決定系數R2、均方根誤差RMSE綜合評價土壤含鹽量擬合效果。其中R2越接近于1，表示土壤含鹽量擬合效果越好，RMSE通過預測值和實測值的偏差度判斷模型的準確性，值越小，表示預測值和實測值越接近，AIC和SIC是建立在熵的概念基礎上，衡量統計模型擬合優良性的一種標準，有效避免過擬合現象的發生。AIC和SIC參數值越小表示該模型能夠以最少自由變量最好地解釋因變量[32]。其計算公式為

(2)

(3)

(4)

式中n——樣本數k——自變量數目

L——似然函數

RSS——殘差平方和

2 結果與分析

2.1 自變量與土壤含鹽量的相關性分析

從遙感圖像提取4個波段反射率和12個光譜指數，組成各個深度土壤含鹽量模型構建的樣本數據集，隨機選擇2/3的樣本數據組成估算數據集(48個樣本)，與土壤含鹽量進行相關性分析，其結果如圖2所示。參考相關系數檢驗臨界值進行變量的顯著性檢驗，當自由度為48，相關系數的絕對值大于0.365時，達到0.01顯著性水平。從圖2可知，B1、B2、B3、SI、SI1、SI3、S1、S2、RVI、NDVI、MSAVI、ARVI在各個深度下與土壤含鹽量的相關系數的絕對值均大于0.365，達到0.01顯著性水平。

圖2 自變量與土壤含鹽量的Pearson相關系數分析結果Fig.2 Results of Pearson correlation coefficient analysis between independent variables and soil salinity

2.2 基于全子集的最優自變量組合方式篩選確定

利用全子集篩選法列舉不同數目自變量的隨機組合方式；通過驗證集R2確定不同深度下不同數目自變量的最佳組合方式見表3；根據R2、RMSE、AIC、SIC4種不同的驗證指標確定不同深度下最優自變量組合方式。

從表3可以看出，在同一深度下隨著自變量數目的增加，R2呈先增大后減小的趨勢，AIC、SIC、RMSE呈先減小再增大的趨勢。在0～20 cm、20～40 cm下，RMSE隨著深度增加而減小，這是由于上層土壤容易受大氣、人為等外界環境影響造成土壤含鹽量變異性大，隨著深度增加外部環境的影響基本消除精度逐漸升高。總體來說，R2呈現先增大后減小的趨勢，在20～40 cm處達到最大值。結合表2、3可以看出，在相關性分析中B4和BI與不同深度土壤含鹽量沒有呈現極顯著的相關關系，經過全子集篩選后B4和BI的組合方式在不同深度下與土壤含鹽量均呈現良好的相關性，并成為各深度下自變量數目為2的最優自變量組合。在植被遙感中應用最為廣泛的歸一化植被指數NDVI與不同深度的土壤含鹽量均呈現極顯著的相關性，在20～40 cm、0～40 cm、40～60 cm、0～60 cm土壤深度下與B4、BI組合后成為該深度下自變量數目為3的最優自變量組合。在0～20 cm處，其最優組合方式為B4、BI、MSAVI，這是由于在0～20 cm處多是土壤和主要作物的側根，土壤背景和植被的相互作用通過改進型土壤調整植被指數MSAVI減小土壤亮度的影響[33]。當自變量數目增加到4，主要增加了由B2和B3計算得出的SI1、SI3兩種光譜指數，這是由于土壤一般都有很高的溶解性鹽分，在520～770 nm平均反射率最高[34]，且屬于確定不同鹽漬化過程中積鹽狀態特征的6個光譜區間之一[35]。隨著模型復雜程度的增加，AIC和SIC逐漸增大，但R2減小、RMSE增大，模型靈敏性降低。綜合分析全子集篩選的各個評價指標，確定在0～20 cm、0～40 cm處選擇B4、BI、SI1、SI3共4個自變量，在20～40 cm、40～60 cm、0～60 cm處選擇B4、BI、NDVI共3個自變量作為各深度下最優自變量組合方式。

表3 全子集篩選最佳組合方式結果統計Tab.3 Best combination result after total subset selection

注：** 表示0.01顯著性水平，*表示0.05顯著性水平。

2.3 人工神經網絡模型建立與分析

分別以不同深度下篩選前后的反射率及光譜指數為自變量、土壤含鹽量為因變量，運用人工神經網絡模型進行篩選前后不同深度下土壤含鹽量估算，其建模及驗證結果如表4所示。

表4 基于不同深度土壤含鹽量的人工神經網絡模型Tab.4 ANN model of soil salinity at different depths

2.4 支持向量機模型建立與分析

以不同深度下篩選前后的反射率及光譜指數為自變量，土壤含鹽量為因變量，運用支持向量機模型進行篩選前后不同深度下土壤含鹽量估算。其建模和驗證結果如表5所示。

表5 基于不同深度土壤含鹽量的支持向量機模型Tab.5 SVM model of soil salinity at different depths

2.5 分位數回歸模型建立與分析

對于不同深度選定的敏感植被指數建立分位數回歸模型，由于采樣點土壤含鹽量變異性較大，選取τ= 0.5分位點可以較好地解決最小二乘法中某些“離群值”影響回歸顯著性的問題。同時，由于0.5分位點處于因變量的中間位置，在對所有的數據進行擬合時較為適宜[36]。故本文以不同深度下篩選前后的反射率及光譜指數為自變量，以土壤含鹽量為因變量，運用分位數回歸模型中的0.5分位點進行篩選前后不同深度下土壤含鹽量估算，其建模及驗證結果如表6所示。

表6 基于不同深度土壤含鹽量的分位數回歸模型Tab.6 QR model of soil salinity at different depths

2.6 模型綜合評價

綜上所述，人工神經網絡模型的反演效果最差，分位數回歸模型和支持向量機模型在本次土壤含鹽量估算中就模型效果而言，二者基本相同。但是分位數回歸模型僅采用3個自變量，支持向量機模型采用16個自變量。因此，本文選擇擬合效果好、驗證精度高、模型穩定性強、簡潔高效的20～40 cm深度的全子集-分位數回歸模型作為本次最佳土壤鹽漬化估算模型。

圖3 土壤含鹽量空間分布圖Fig.3 Soil salinity distribution map

基于全子集-分位數回歸方法建立模型，進行解放閘灌域土壤含鹽量估測，結果如圖3所示。圖中白色部分為城鎮所在區域不參與本次含鹽量反演。在該時期，灌域內鹽漬化程度較輕，有利于作物的生長發育。灌域內輕度鹽漬化土所占比例最大為54%，非鹽土所占比重次之，為26%。這是由于該時期處于作物生長關鍵期，灌水次數多且灌水量大，土壤鹽分淋洗到下層土壤，土壤處于夏季脫鹽狀態。重鹽漬化土所占比重較少，為13%，對于灌水次數較少的地方，土壤仍產生積鹽現象。鹽土占總面積的比例最少，為7%，主要集中在灌域西部(多為鹽荒地)，其余鹽土零散分布于灌域內，主要包括部分鹽荒地和棄耕地。

3 討論

土壤鹽漬化問題在中國北方干旱灌區較為突出，已成為制約灌區農業可持續發展的重要因素。黃河流域上游的河套灌區是受鹽漬化影響的典型區，土壤鹽漬化面積約占總面積的69%[37]。利用衛星遙感對土壤含鹽量進行大面積監測，符合未來精準農業發展的要求。在本次植被覆蓋條件下，作物長勢可間接識別土壤根域鹽漬化水平，研究結果顯示20～40 cm深度下的土壤含鹽量與高分一號衛星數據相關性最高，土壤鹽分含量對作物主根系的影響明顯。史海濱等[38]對比各生育期內不同含鹽處理的各層主根重的變化情況發現，向日葵的苗期后期(7月27日)主根深度在30 cm之內，土壤鹽分含量對作物主根系層的生長影響明顯。與本文研究成果基本一致。由此可見，應用衛星遙感大范圍監測植被覆蓋條件下不同土層的鹽漬化水平具有一定的可行性。

采用相關性分析只能獲得單一自變量與含鹽量的相關關系，本研究采用全子集法列舉全部組合方式，得到不同深度下對于土壤含鹽量最優自變量數目的敏感變量組合，篩選后的結果與未篩選的結果在擬合優度和模型精度方面略有提升。表明全子集篩選過程簡單、高效、靈敏度強，適用于本次鹽分含量的估測。OKCU等[39]采用全子集法生成一個新的適用于河流和水槽數據的泥沙運移公式，精準預測泥沙運移的非線性關系。但是由于自變量數目的增加，全子集篩選的運行效率減弱，HOFMANN等[40]改進和擴展全子集回歸的計算方法，采用基于回歸樹的RadiusBBA和HBBA方法解決大規模模型選擇問題的窮舉和啟發式策略。全子集法和相關性分析法的篩選結果可作為自變量尋優的一種參考，而篩選結果與因變量是否有內在聯系，需要做更嚴謹的物理機制和數學推理工作。

本文對于不同深度選定的敏感指數建立的ANN、SVM、QR 3種模型進行土壤含鹽量反演，發現QR模型的反演精度最高。這是由于QR采用非對稱權重的方法估計參數，適用于離群值多、變量間關系較弱和需要了解因變量分布的情況，且無需對模型進行任何分布的假定，對于強變異性(變異系數約100%)的實測土壤鹽分具有抗耐性。AMAKOR[19]進行鈣質鹽堿土電導率模擬時，實測電導率變異系數在100%左右，采用QR模型決定系數R2達到了0.88。MUELLER等[41]、GERBER等[42]和CARSLAW等[43]采用QR模型，均很好地解決了實測數據重尾分布、變異性強的問題。τ=0.5分位數處于土壤含鹽量的中間位置，采用該分位數進行各個深度的土壤鹽漬化反演可以保證輕、重度鹽漬化土權重較大，非鹽土和鹽土權重較小，從而減小離群值對于模型精度的影響。前人也得到類似的結果，如王蕾等[36]發現在0.5分位數下進行冬小麥估測結果最為可靠。

研究植被覆蓋條件下不同深度的解放閘灌域土壤鹽漬化模型，對改善灌區土壤鹽漬化情況具有現實意義。然而最優土壤含鹽量反演模型會因作物的種類、生育期、天氣狀況、灌溉情況、甚至所使用的遙感平臺而異。本文所得的最佳含鹽量反演模型也僅限于本次測量結果，研究區內其他生育階段和其他地區還可嘗試更多圖像特征(如顏色特征[44]、地形特征[45])，尋找更優的含鹽量反演指標，從而建立精度更高且實用性更強的土壤含鹽量反演模型。

4 結論

(1)全子集篩選具有全局性，經過全子集篩選后敏感變量組合方式不盡相同。B4、BI、SI1、SI3是0～20 cm、0～40 cm處土壤含鹽量的敏感變量組合，B4、BI、NDVI為20～40 cm、40～60 cm、0～60 cm處土壤含鹽量的敏感變量組合。篩選后的敏感變量組合方式僅用3～4個自變量，與篩選之前的效果相差不多，因而全子集篩選具有靈敏度強、簡單、高效的優勢。

(2)在不同深度下的土壤含鹽量反演模型中，分位數回歸、人工神經網絡、支持向量機模型均取得良好的建模驗證精度。分位數回歸模型由于其抗異性較強，其精度與支持向量機模型相差不多，優于人工神經網絡模型。土壤深度對含鹽量反演模型精度具有一定的影響。在使用相同的建模方法時，20～40 cm深度下反演模型效果優于其他深度。

(3)在20～40 cm深度下建立的模型效果最為理想，且分位數回歸抗異性、穩健性強，在保證模型估測精度的前提下，比支持向量機模型和人工神經網絡模型操作簡單、結果可靠。因此，全子集-分位數回歸模型是本次估算植被覆蓋條件下土壤含鹽量的最優模型。