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

楊 練 兵，鄭 宏 偉,4*，羅 格 平,4，楊 遼,4

(1.中國科學院新疆生態與地理研究所荒漠與綠洲國家重點實驗室,新疆 烏魯木齊 830011;2.新疆維吾爾自治區遙感與地理信息系統應用重點實驗室,新疆 烏魯木齊 830011;3.中國科學院大學,北京 100049;4.中國科學院中亞生態與環境研究中心，新疆 烏魯木齊 830011)

0 引言

土壤鹽漬化是土地退化和荒漠化的主要類型之一，其分布范圍廣，危害時間長，嚴重制約著干旱區農業生產，已成為全球性生態環境問題[1]。衛星遙感數據能宏觀、動態地監測土壤鹽漬化信息[2]，相關研究經歷了定性分類和定量反演兩個階段。在定量反演方面，眾多學者對獲得的遙感數據和非遙感數據進行數學變換以構建反演參數，并結合采樣點土壤鹽分含量(Soil Salinity Content,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×104km2，地勢東南高、西北低，高程為256～5 242 m，有典型的山地—綠洲—荒漠系統。該流域深居內陸干旱區，氣候干燥，光照充足，年均氣溫 5～7 ℃，年降水量110～200 mm，年蒸發量1 500～2 000 mm[20]，水資源主要來源于冰雪融水和高山降水，農業生產起初的大水漫灌引發嚴重的土壤鹽漬化問題，改為滴灌后土壤鹽漬化程度有所減輕[21]。

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

1.2 采樣點數據

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

圖1 兩靶區樣點分布

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

表1 兩靶區采樣點SSC統計描述

1.3 影像數據與反演參數

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

表2 SSC反演參數

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 模型設計流程及研究內容

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

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

圖3 結構參數染色體示意

(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)。

圖4 BPNN結構圖示

(1)

p=2n+1

(2)

p=log2n

(3)

(4)

式中：p、m、n分別為隱含層、輸出層、輸入層的神經元數量；q用于調整隱含層神經元數量，取值范圍為1～10之間的整數。

圖5 GA優化BPNN結構參數示意

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

圖6 初始權重染色體示意

圖7 GA優化BPNN初始權重示意

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

VIP算法是基于偏最小二乘回歸的特征篩選方法，可評價自變量對因變量集合的解釋能力，公式為：

(5)

式中：k為自變量數量；th為從自變量集合X=(x1,x2,…,xj,…,xk)中提取的主成分；Y為因變量集合；whj用于衡量自變量xj對th的邊際貢獻；Rd(Y;th)為th對Y的解釋能力，Rd(Y;t1,t2,…,tm)為t1,t2,…,tm對Y的累積解釋能力。VIPj大于1，表示xj對Y非常重要；VIPj在0.5～1之間，表示xj對Y的重要性不明確，需要根據其他條件進行判斷或增加樣本；VIPj小于0.5，表示xj對Y不重要[45]。

VIP1-GA-BP和VIP2-GA-BP的建模流程(略)與GA-BP相似，區別為： GA-BP中輸入層反演參數子集和隱含層神經元數量的優化同步進行，而VIP1-GA-BP和VIP2-GA-BP對這兩者的優化是分開進行的。

2.3 精度檢驗指標

本研究采用均方根誤差(RMSE)、平均絕對百分誤差(MAPE)和相對分析誤差預測偏差(RPD)對模型性能進行評價，各指標的計算公式如式(6)-式(8)所示。RMSE越小，表示模型的預測精度越高；MAPE越接近于0，表明模型預測的相對誤差越小；RPD<1.4，說明模型不可靠，1.42 ，說明模型精度較高。

(6)

(7)

RPD=SD/RMSE

(8)

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

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

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

圖8 SSC反演參數變量投影重要性排序

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最優適應度值之間的差值均較小。

圖9 GA優化BPNN結構參數的最佳適應度曲線

對兩靶區各組模型最終結構參數種群中適應度值最高個體的染色體進行解碼，確定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組模型建模集間RMSE、MAPE的差異均較小，而測試集間的差異均較大，其中GA-BP模型在測試集上的RMSE和MAPE均最小(分別為11.56 g/kg、32.78%)、RPD最大(1.64)，表明GA-BP模型的反演精度最高；在三工河靶區，GA-BP模型的反演精度也最高，其RMSE、MAPE、RPD分別為3.60 g/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。

表3 兩靶區SSC反演精度統計

上文已對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 不同模型反演參數子集

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

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

圖10 瑪納斯靶區3組模型反演的SSC分布

圖11 三工河靶區3組模型反演的SSC分布

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

圖12 瑪納斯靶區子區的位置及SSC分布

圖13 三工河靶區子區的位置及SSC分布

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

5 結論與討論

在新疆瑪納斯流域和三工河流域各選一靶區，基于Landsat-8 OLI、Sentinel-1 SAR影像數據和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)瑪納斯靶區土壤鹽漬化程度較三工河靶區嚴重，兩靶區面積占比最大的鹽漬土類型分別為鹽漬土和中度鹽漬土。

圖14 兩靶區各類鹽漬土的空間分布

圖15 兩靶區各類鹽漬土的面積比例

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