基于遺傳算法優化ＢＰ神經網絡的土壤鹽漬化反演

　　摘要：應用于土壤鹽分含量(SoilSalinityContent，SSC)反演的 BP神經網絡(BackPropagationNeuralNetwork，BPNN)較少關注對模型精度影響較大的結構參數和初始權重的優化。該文利用Landsat8OLI、Sentinel1SAR影像數據及SRTM 高程數據，基于谷歌地球引擎(GEE)平臺構建反演參數，并建立3種反演模型：先利用遺傳算法(GeneticAlgorithm，GA)同步優化輸入層反演參數子集和隱含層神經元數量，再優化初始權重的BPNN(GABP)模型;將變量投影重要性(VariableImportanceinProjection，VIP)算法分割閾值分別設為1和0.5，優化出兩組輸入層反演參數子集并將其分別代入GA優化隱含層神經元數量，再優化初始權重的BPNN(VIP1GABP、VIP2GABP)模型。在瑪納斯流域和三工河流域各選一靶區進行SSC反演，對比分析 GABP、VIP1GABP、VIP2GABP模型的反演精度，并統計各類鹽漬土的面積比例，結果表明：1)兩靶區3組模型反演精度由高到低排序均為 GABP、VIP1GABP、VIP2GABP;2)鹽分指數和植被指數在SSC反演中起到重要作用，同一模型篩選的反演參數存在空間分異性，但高程適用于不同的篩選模型，具有較強的魯棒性;3)兩靶區3組模型反演的SSC值域范圍與實際采樣點SSC值域范圍的差異均較小，各子區GABP反演的SSC空間分布地物輪廓最清晰，且地物內SSC的均質性最好;4)瑪納斯靶區和三工河靶區面積占比最大的鹽漬土類型分別為鹽漬土和中度鹽漬土。研究結果為構建具有一定推廣性的干旱區土壤鹽分含量反演模型奠定了基礎。

　　楊練兵; 鄭宏偉; 羅格平; 楊遼， 地理與地理信息科學 發表時間：2021-04-09 16:56 期刊《地理與地理信息科學》(雙月刊)創刊于1985年，由河北省地理科學研究所主辦。獲獎情況：全國《中文核心期刊要目總覽》核心期刊;河北省第六屆優秀科技期刊;中國科技論文統計源期刊。

　　關鍵詞：土壤鹽分含量;BP神經網絡;遺傳算法;同步優化;反演參數

　　土壤鹽漬化是土地退化和荒漠化的主要類型之一，其分布范圍廣，危害時間長，嚴重制約著干旱區農業生產，已成為全球性生態環境問題[1]。衛星遙感數據能宏觀、動態地監測土壤鹽漬化信息[2]，相關研究經歷了定性分類和定量反演兩個階段。在定量反演方面，眾多學者對獲得的遙感數據和非遙感數據進行數學變換以構建反演參數，并結合采樣點土壤鹽分含量(SoilSalinityContent，SSC)數據建立相應的反演模型，以獲得大范圍的土壤鹽分含量。

　　受自然條件和人類活動方式的影響，土壤鹽漬化成因復雜[3]，光譜波段、植被指數、鹽分指數、下墊面參數、特征空間、地形參數、物候特征等都被作為反演土壤鹽分含量的參數[4，5]。相關研究采用多元線性回歸(MLR)、偏最小二乘回歸(PLSR)、地理加權回歸(GWR)、分位數回歸(QR)等線性模型擬合反演參數與土壤鹽分含量的關系[6-8]。由于多數反演參數與土壤鹽分含量之間存在非線性關系[9]，BP神經網絡(BPNN)、支持向量機(SVM)、多元自適應回歸樣條(MARS)、隨機森林(RF)等非線性模型也被用于反演土壤鹽分含量[10-14]。由于部分機器學習算法自身沒有特征優選的能力，皮爾森相關分析[15，16]、最佳指數法[17]、灰色關聯分析[18]等過濾方法被用于篩選反演參數子集，這種特征篩選方式減少了信息冗余，但難以獲得最佳反演參數子集;同時，機器學習模型的參數影響算法的預測精度、運算速度和穩健性，已有研究多通過經驗方法選取模型參數，調參效率低且難以獲得最優參數;另外，應用于土壤鹽漬化反演的機器學習模型種類較多，以往多關注不同模型間的對比，對同一模型的優化改進較少研究。

　　BPNN具有優良的非線性逼近能力，是較早引入土壤鹽漬化反演的機器學習模型[13]，但其預測精度受結構參數和初始權重的影響[19]。為此，本研究建立3組優化模型：先利用遺傳算法(GA)同步優化輸入層反演參數子集和隱含層神經元數量，再優化初始權重的BPNN(GABP)模型;將變量投影重要性(VIP)算法分割閾值分別設為1和0.5，優化出兩組輸入層反演參數子集，并將其分別代入 GA 優化隱含 層 神 經 元 數 量，再 優 化 初 始 權 重 的 BPNN(VIP1GABP、VIP2GABP)模型。在新疆瑪納斯流域和三工河流域各選一靶區，基于谷歌地球引擎(GEE)平臺構建反演參數，對比分析兩靶區3組優化模型的反演精度，并統計各類鹽漬土的面積比例，以期為土壤鹽漬化信息的高效獲取提供支持。

　　1 研究區與數據

　　1.1 研究區概況

　　瑪納斯流域(85°01′～86°32′E，43°27′～45°21′N)位于天山北麓(圖1a)、準噶爾盆地南緣，由瑪納斯河、塔西河、寧家河、金溝河、巴音溝河和大南溝河組成，流域面積為3.35×104 km2，地勢東南高、西北低，高程為256～5242m，有典型的山地—綠洲—荒漠系統。該流域深居內陸干旱區，氣候干燥，光照充足，年均氣溫5～7℃，年降水量110～200mm，年蒸發量1500～2000mm[20]，水資源主要來源于冰雪融水和高山降水，農業生產起初的大水漫灌引發嚴重的土壤鹽漬化問題，改為滴灌后土壤鹽漬化程度有所減輕[21]。

　　三工河流域(87°47′～88°17′E，43°09′～45°29′N)位于天山北麓中段東部(圖1b)、準噶爾盆地南緣，主要由水磨溝河、三工河和四工河組成，流域面積為1.67×103km2，總人口約為11萬人[22]，高程為480～650 m，地 勢 由 東 南 向 西 北 傾 斜。該 流 域 降 水稀少(約220mm/a)，蒸發量大(約1399mm/a)，年均溫約為7.3℃[23]，以種植業為主，灌區水資源主要來源于冰雪融水，地下水礦化度高，地下水位抬升及強烈的蒸發作用容易導致地表積鹽[24]。

　　1.2 采樣點數據

　　2016年7-8月，在瑪納斯流域和三工河流域各選一靶區，用手持 RTKGPS進行采樣點定位，采樣點考慮不同的植被覆蓋、地貌類型和交通可達性，具有一定的代表性。在兩靶區環境因素相似的區域設置30m×30m樣方，在樣方內按照“五點采樣法”采集土壤樣品，采樣深度為 0～20cm;采集的土壤樣品經 過 自 然 風 干、磨 碎、過 篩 后，測 定 八 大 離 子(Ca2+ ，Mg2+ ，K+ ，Na+ ，CO32- ，HCO3- ，Cl- ，SO42- )的含量，對其求和獲得土壤樣品鹽分含量(SSC)，將樣方內土壤樣品的SSC均值作為實際觀測值。針對測量結果中的誤差，基于箱線圖分析用四分位差(IQR)檢測異常值點，最終在瑪納斯和三工河靶區獲得可用采樣點數量分別為97個和119個(圖1)。

　　在瑪納斯靶區和三工河靶區分別隨機選取總樣本集的26.80%和26.05%(約1/4)作為測試集，其余樣本數據作為建模集。由總樣本集、建模集、測試集的描述性統計特征(表1)可知：兩靶區總樣本集的均值和變異系數均介于建模集與測試集之間，表明兩靶區建模集和測試集樣本數據的范圍相對一致，一定程度上避免了模型構建和驗證中的偏差估計[25]。假設樣本數據能代替總體，按照變異系數(CV)評價準則[26]，兩靶區SSC均屬于中等變異;按照新疆鹽堿土標準[27]，從均值看，瑪納斯靶區和三工河靶區土壤分別屬于重度鹽漬化和中度鹽漬化。

　　1.3 影像數據與反演參數

　　谷歌地球引擎(GEE)是目前較為成熟的遙感大數據分析云平臺[28]，可提供高性能并行計算服務。本研究參考相關文獻，調用 GEE平臺中預處理好的Landsat8OLI影像數據(空間分辨率30m)、Sentinel1SAR 影像數據(空間分辨率5m×20m)及SRTM 高程數據(空間分辨率30m)計算出兩靶區的52個反演參數(表2)，并將廣義差分植被指數(GDVI)、土壤調節植被指數(SAVI)、增強植被指數(EVI)、綠色大氣阻抗指數(GARI)的反演參數設為經驗值;考慮到Landsat8OLI各波段的分辨率和特性，只選取可能對土壤鹽漬化具有表征能力的可見光、近紅外和短波紅外波段作為反演參數。由于本研究SSC采樣點數據來源于30m×30m樣方，故通過三次卷積內插將反演參數空間分辨率統一到30m。為接近SSC采樣時間且滿足靶區范圍內無云的條件，選取瑪納斯靶區和三工河靶區成像時間分別為2016年8月4日和8月29日的Landsat8OLI數據以及2016年8月2日和8月26日的Sentinel1SAR數據;SRTM 高程數據制作時間為2000年2月。

　　2 研究方法

　　本研究基于BP神經網絡(BPNN)，通過遺傳算法(GA)、變量投影重要性(VIP)算法共設計出3組優化的BPNN模型，具體研究流程(圖2)為：1)先利用GA同步優化輸入層反演參數子集和隱含層神經元數量，再優化初始權重的BPNN(GABP)模型;2)將 VIP算法分割閾值分別設為1和0.5，優化出兩組輸入層反演參數子集，將兩組輸入層反演參數子集分別代入GA，先優化隱含層神經元數量，再優化初始權重的BPNN(VIP1GABP、VIP2GABP)模型;3)對兩靶區3組模型反演SSC的精度、空間分布、篩選的反演參數進行對比，并統計各類鹽漬土的面積比例。

　　2.1 遺傳算法(GA)優化的BP神經網絡(GABP)

　　GA是 一 種 模 擬 生 物 進 化 過 程 的 啟 發 式 算法[40]，其將可能的解轉換為種群中個體的染色體(用二進制、實數、十進制、格雷、符號等編碼方式以符號串形式表示)，染色體的基因位取值區間需根據具體問題設定;GA 有選擇、交叉、變異3種遺傳操作算子，以適應度函數為評價指標，通過遺傳操作算子不斷更新種群，解碼最終種群中適應度值最高個體的染色體為問題的最優解。GABP結構參數染色體(圖3)中，犻犿(犿=1，…，狀，狀為待篩選反演參數的數量)代表輸入層神經元的基因位，取值為1(表示對應的反演參數參與建模)或0(表示反演參數不參與建模);犺表示隱含層神經元數量的基因位。本研究采用兩階段方式構建 GABP：先用 GA 對BPNN的結構參數進行優化，然后對BPNN 的初始權重進行優化。

　　(1)BPNN 結構參數的優化。BPNN 包括輸入層、隱含層、輸出層3層網絡結構[41]，基本組成單位為神經元，通過在隱含層和輸出層設置激活函數解決非線性擬合問題。當BPNN僅含一個隱含層且隱含層神經元數量較多時，其具備很強的函數逼近或映射能力[42]，故本研究用于SSC預測的BPNN只設一個隱含層(圖4)。本文中BPNN結構參數為輸入層神經元(輸入層反演參數子集)、隱含層神經元數量和輸出層神經元數量，由于輸出層神經元數量為1，故只需對輸入層反演參數子集和隱含層神經元數量進行優化。BPNN隱含層神經元數量與輸入層神經元數量有關[43]，而 GABP各結構參數染色體代表的輸入層神經元數量不同，故隱含層神經元數量的基因位取值區間需動態更改，研究中依據經驗公式(式(1)-式(3)[43]、式(4)[44])聯立取最大值和最小值確定。為實現隱含層神經元數量基因位取值區間的動態更新，本研究在 GABP中對結構參數染色體的交叉算子和變異算子進行改進：在結構參數染色體進行交叉和變異運算時，先對輸入層神經元基因進行運算;然后根據其代表的神經元數量，推導出隱含層神經元數量的基因位取值區間;最后在該取值區間中隨機生成一整數并將其作為隱含層神經元數量的基因位數值。通過遺傳算子不斷更新BPNN結構參數種群，解碼最終結構參數種群中適應度值最高個體的染色體，即可獲得最優的結構參數(圖5)。

　　式中：狆、犿、狀分別為隱含層、輸出層、輸入層的神經元數量;狇用于調整隱含層神經元數量，取值范圍為1～10之間的整數。

　　(2)BPNN 初始權重的優化。初始權重染色體如圖6所示，狑狆(狆=1，…，狇，狇為初始權重數量)代表初始權重的基因位。確定BPNN最優的輸入層反演參數子集和隱含層神經元數量后，利用 GA 的遺傳算子不斷更新初始權重種群，解碼最終初始權重種群中適應度值最高個體的染色體，即可獲得最優的初始權重(圖7)。

　　2.2 變量投影重要性(VIP)算法

　　VIP算法是基于偏最小二乘回歸的特征篩選方法，可評價自變量對因變量集合的解釋能力，公式為：犞犐犘犼= 槡犽/犚犱(犢;狋1，狋2，…，狋犿)∑ 犿犺=1犚犱(犢;狋犺)狑2犺犼 (5)式中：犽為自變量數量;狋犺 為從自變量集合犡=(狓1，狓2，…，狓犼，…，狓犽)中提取的主成分;犢 為因變量集合;狑犺犼用于衡量自變量狓犼 對狋犺 的邊際貢獻;犚犱(犢;狋犺)為狋犺 對犢的解釋能力，犚犱(犢;狋1，狋2，…，狋犿)為狋1，狋2，…，狋犿 對犢的累積解釋能力。犞犐犘犼 大于1，表示狓犼 對犢非常重要;犞犐犘犼 在0.5～1之間，表示狓犼 對犢 的重要性不明確，需要根據其他條件進行判斷或增加樣本;犞犐犘犼 小于0.5，表示狓犼 對犢不重要[45]。VIP1GABP和VIP2GABP的建模流程(略)與GABP相似，區別為：GABP中輸入層反演參數子集和隱含層神經元數量的優化同步進行，而 VIP1GABP和VIP2GABP對這兩者的優化是分開進行的。

　　2.3 精度檢驗指標

　　本研究采用均方根誤差(犚犕犛犈)、平均絕對百分誤差(犕犃犘犈)和相對分析誤差預測偏差(犚犘犇)對模型性能進行評價，各指標的計算公式如式(6)-式(8)所示。犚犕犛犈 越 小，表 示 模 型 的 預 測 精 度 越 高;犕犃犘犈越接近于0，表明模型預測的相對誤差越小;犚犘犇<1.4，說明模型不可靠，1.4<犚犘犇<2，說明模型精度一般，犚犘犇>2，說明模型精度較高。

　　式中：狀為樣本數量;狔犻 為樣本實測值;^狔犻 為樣本預測值;犛犇 為數據集的標準差。

　　3 反演參數篩選及模型參數設置

　　在兩靶區以建模集及其對應的反演參數為數據源，基于SIMCA軟件進行VIP分析;根據犞犐犘 值對解釋變量的重要性意義，設置分割閾值為1和0.5，將犞犐犘 值大于或等于1的反演參數集記為 A組，大于或等于0.5的記為B組(圖8)。在5%的顯著性水平上，瑪納斯靶區共有22個反演參數的犞犐犘 值大于1(包含三工河靶區犞犐犘 值大于1的反演參數)，13個在0.5～1之間;三工河靶區反演參數的犞犐犘 值均大于0.5，13個大于1。

　　本 研 究 采 用 MATLAB 遺 傳 算 法 工 具 箱(GAOT)設計 GA，GA 的編碼方式為實數編碼，建立優化的BPNN模型時，將測試集樣本數據平均絕對誤差的倒數作為 GA 優化 BPNN 結構參數(GA對 VIP1GABP、VIP2GABP優化的結構參數均僅為隱含層神經元數量)和BPNN 初始權重的適應度函數。兩靶區3組模型的 BPNN 訓練次數均設為1000，訓練目標設為0.02，學習速率設為0.03，初始權重基因位取值范圍設為-1～1。為避免各維度數據間數量級的差別，兩靶區訓練集和測試集及其對應的反演參數數值均歸一化至0～1之間。tansig能將數值映射到-1～1之間，便于BPNN對輸入和輸出數據的非線性擬合，logsig能將輸出數值映射到0～1之間，便于BPNN 后續的反歸一化處理，故將兩靶區3組模型的BPNN隱含層神經元激活函數均設為tansig，將BPNN 輸出層神經元激活函數均設為logsig。

　　4 結果分析

　　4.1 模型的預測精度和反演參數

　　隨著GA優化BPNN結構參數的進化代數逐步增加，兩靶區各組模型種群中個體適應度值也逐漸增大，最后趨于穩定(圖9);兩靶區最優適應度值從大到小排序分別為 GABP>VIP2GABP>VIP1GABP和 GABP>VIP1GABP>VIP2GABP，可見 VIP1GABP、VIP2GABP最優適應度值之間的差值均較小。

　　對兩靶區各組模型最終結構參數種群中適應度值最高個體的染色體進行解碼，確定 VIP1GABP、VIP2GABP的隱含層神經元數量及 GABP的輸入層反演參數子集(C組反演參數)和隱含層神經元數量，兩靶區 VIP分析確定的 A 組、B組反演參數分別作為 VIP1GABP和 VIP2GABP的輸入層反演參數子集;之后采用GA對初始權重進行優化，建立3組優化的 BPNN 模型。由兩靶區3組模型的SSC反演精度(表3)可以看出：瑪納斯靶區3組模型建模集間犚犕犛犈、犕犃犘犈 的差異均較小，而測試集間的差異均較大，其中 GABP模型在測試集上的犚犕犛犈 和 犕犃犘犈 均 最 小 (分 別 為 11.56g/kg、32.78%)、犚犘犇 最大(1.64)，表明 GABP模型的反演精度最高;在三工河靶區，GABP模型的反演精度也 最 高，其 犚犕犛犈、犕犃犘犈、犚犘犇 分 別 為 3.60g/kg、25.21%、1.47。綜上，3組模型反演精度由高至低 依 次 為 GABP、VIP1GABP、VIP2GABP。由 VIP算法原理可知，B組比 A組多出的反演參數對被解釋變量SSC的重要性不明確，信噪比較小，會造成 VIP2GABP較嚴重的過擬合，一定程度上解釋了兩靶區 VIP1GABP的反演精度均高于 VIP2GABP。GABP同步篩選了反演參數子集和隱含層神經元數量，考慮了反演參數之間及反演參數與反演模型間的相互關系，最終篩選的反演參數子集與BPNN耦合性高;而 VIP算法雖對反演參數進行重要性排序，但未考慮參數間的相互關系及反演SSC的特定建模方法[46]，故兩靶區 GABP的反演精度均高于 VIP1GABP和 VIP2GABP。

　　上文已對 A組、B組反演參數進行對比分析，且基于 A組反演參數的模型反演精度更高，故僅需對各靶區 A 組、C組反演參數進行對比分析(表4)。兩靶區 A 組、C組反演參數共有43種，鹽分指數最多(14種)，其次為植被指數(13種)，說明鹽分指數和植被指數在SSC反演中發揮著重要作用;另外，不同模型篩選的反演參數差異較大，同一模型篩選的反演參數存在區域異質性。瑪納斯靶區 A 組、C組反演參數數量均為22，共有的反演參數為 RVI、S1、Albedo_short、AVI、Elevation、Roughness;三工河靶區 A組、C組反演參數數量分別為13、18，共有反演參數為GARI、S2、S3、S6、Elevation、Swir2。兩靶區 A組、C組共有的反演參數僅為 Elevation，說明高程參數不僅適用不同的篩選模型，而且適用不同區域鹽漬化土壤研究，原因在于：地表徑流為地勢較低的土壤表層帶來鹽分，且地勢較低的地方潛水埋深較淺，潛水更易通過蒸發作用使地表積鹽[47];同時，高程不易受其他因素影響，且短時間內變化較小。

　　4.2 反演的SSC空間分布及各類鹽漬土面積比例

　　由于本研究僅對土壤進行鹽分反演，故對兩靶區3組模型反演結果中的水域、建筑用地進行掩膜處理(圖10、圖11)。可以看出，兩靶區3組模型反演的SSC空間分布存在一定差異，且圖斑的破碎度較大，而SSC值域范圍與實際采樣點SSC值域范圍的差異均較小。瑪納斯靶區中，GABP 和 VIP2GABP預測的SSC值域范圍(分別為0.20～70.00g/kg、0.18～69.99g/kg)更接近實際采樣點SSC值域范圍(0.18～70.00g/kg)，三工河靶區中，GABP預測的SSC值域范圍(2.80～26.39g/kg)更接近實際采樣點SSC值域范圍(2.62～26.56g/kg)。

　　土壤屬性的空間分布具有空間自相關和空間分異性[48]。為比較3組模型反演SSC的局部特征，在兩靶區隨機選取3個子區進行對比(圖12、圖13)，發現各子區3組模型反演的SSC空間分布均存在較大差異(可能與模型篩選的反演參數有關)，與各模型反演的SSC整體空間分布結果對應;各子區中GABP反演的各類地物輪廓最清晰，且地物內SSC的均質性也最好。

　　由上述分析可知，GABP反演精度最高，故參照新疆鹽堿土分類標準[27]，基于 GABP 反演的SSC將兩靶區土壤劃分為5類(圖14，彩圖見附錄3)，并統計各類鹽漬土的面積比例(圖15)。在瑪納斯靶區，鹽漬土面積占比最高(55.87%)，其次為非鹽漬土(16.60%)，二者多呈塊狀分布，重度鹽漬土、中度鹽漬土、輕度鹽漬土多呈點狀分布。在三工河靶區，中度鹽漬土面積占比最高(51.02%)，其次為輕度鹽漬土(16.60%)，鹽漬土面積占比最小(1.54%)，各類鹽漬土多呈塊狀分布，鹽漬土、重度鹽漬土主要分布在靶區北部。從鹽漬土分類結果及各類鹽漬土面積占比看，瑪納斯靶區土壤鹽漬化程度較三工河靶區嚴重，這與兩靶區采樣點SSC的統計結果吻合。

　　5 結論與討論

　　在新疆瑪納斯流域和三工河流域各選一靶區，基于 Landsat8OLI、Sentinel1SAR 影像數據和SRTM 高程數據構建反演參數，通 過 VIP、GA、BPNN建立3組優化模型，進行SSC反演并統計各類鹽漬土的面積比例。結論如下：1)模型反演精度由高到低排序為 GABP、VIP1GABP、VIP2GABP，表明同步優化反演參數和模型參數的特征篩選方式效果最好，這與 Xu等[11]的研究結果相似。2)鹽分指數和植被指數在土壤鹽漬化反演中起著重要作用，同一模型篩選的反演參數存在區域異質性，但高程適用不同的篩選模型，具有較強的魯棒性，與王飛等[5]的研究結果對應。3)兩靶區3組模型反演的SSC整體空間分布圖圖斑破碎度均較大，SSC值域范圍與實際采樣點SSC值域范圍的差異均較小;各子區中3組模型反演的SSC空間分布均存在較大差異，其中GABP反演的SSC空間分布地物輪廓最清晰，且地物內 SSC 的均質性最好，這與朱阿興等[48]的土壤數字制圖理論相符。4)瑪納斯靶區土壤鹽漬化程度較三工河靶區嚴重，兩靶區面積占比最大的鹽漬土類型分別為鹽漬土和中度鹽漬土。

　　本研究利用GABP進行SSC反演，取得了不錯效果，但GA優化的反演參數子集和模型參數可能是局部最優，今后可借鑒統計學中置信度檢驗的思想，設 定 衡 量 優 化 的 參 數 為 全 局 最 優 可 信 度 指標[49，50];研究中將一些反演參數設為經驗值，算法GA、BPNN也有部分參數是人為設定的，加之本研究沒有兩靶區實際的土壤鹽分含量分布圖，只能通過測試集數據判斷結果優劣，具有一定的不確定性;運用 GEE云計算功能提取反演參數集時，僅考慮VV極化方式的微波物理量，后期可以綜合考慮多種極化方式的微波物理量，采用 Python或JavaScript語言，并移植到 GEE云計算平臺中，進行大范圍智能化的土壤鹽漬化反演研究。