空間關聯隨機森林模型結合Sentinel-2 影像估算濰北地區裸土期土壤鹽分

彭遠新,王澤強,周忠科,宋曉寧,徐夕博

1.棗莊學院旅游與資源環境學院,山東 棗莊 277160

2.山東師范大學地理與環境學院,山東 濟南 250358

土壤是地表生物的各項生命活動正常運行的基礎,在生態系統的物質轉化和能量傳遞中起到關鍵的媒介作用[1-2]。 土壤鹽分含量(SSC,本文均指質量分數)是評價土壤質量的一項重要指標,其對綠色植物生長異常敏感。 土壤在發生鹽化后,鈉離子會發生富集進而限制植物細胞中的物質交換,導致植物生長受限甚至死亡[3-4]。 在全球變暖和人類活動的作用下,全球范圍內約有7%的土地受到不同程度的鹽化影響,1/3 的世界耕地存在鹽漬化的風險,土壤鹽化對干旱區域乃至世界的生態環境安全帶來極大的威脅[5-6]。 所以,為實現土壤鹽化的風險預警和治理恢復,快速精準地進行區域SSC 識別與分異狀態監測尤為重要。

近年來,遙感作為一種新型對地觀測技術,因其具有覆蓋廣、數據更新快和獲取成本低等優勢,被廣泛用于土壤制圖與環境監測領域。 彭杰等[7-10]在實驗室控制條件下利用土壤鹽分光譜定量分析技術,分析出SSC 在450 ～680 nm 和短波紅外的波譜范圍內存在確定的響應波段,并利用其作為模型輸入變量完成SSC 的光譜定量估算,證實SSC 光譜定量估算方式的可行性。 在此基礎上,LIU 等[11-14]分析了衛星影像波段反射率(Landsat-5、Landsat-8、Sentinel-2 和HJ-1 等)與土壤鹽分間的相關關系,探尋光譜波段吸收特征與石膏、碳酸鈣、硫酸鈉等鹽分礦物間微妙且穩定的聯系性,最后建立遙感估算模型,分析得出區域尺度SSC 分布規律。

鹽分在表層土壤中的積累是地貌、降水和人類活動等因素共同作用下的一種綜合表現,因而在利用統計模型擬合SSC 與衛星接收的光譜發射率間的關系時,僅通過簡單線性模型是很難取得較高精度的,而基于非線性回歸技術的機器學習算法在構建SSC 遙感估算模型中具有更大的優勢。 例如,王飛等[15-17]在衛星遙感數據基礎上采用機器學習技術(例如人工神經網絡、支持向量機和隨機森林等)完成了土壤鹽分在區域尺度上的反演制圖。 通常來說,機器學習過程中通過增多隱藏層數量、增加結構參數設置和決策樹節點數目,可以實現對土壤屬性到光譜信號間的信息轉換。 但是,輸入樣本數據在土壤環境、人類活動和采樣過程等因素影響下表現出的特異性和空間非關聯性增加了目標地物與光譜特征間映射關系的復雜度,使機器學習模型的泛化能力減弱和出現局部極值問題,限制了估算精度的提升和模型的可應用性。 針對此問題,方利民等[18-22]探索在遙感估算模型構建之前通過預先對樣本數據進行隨機選擇、含量分布、光譜特征和土壤類型劃分等方式減弱或消除特異性樣本的影響,選取最佳訓練樣本集,在確保驗證集精度的同時有效提升訓練效率和準確度。 上述研究對特異樣本數據的識別均是基于樣本的含量特征信息,然而,采集到的地理樣本除帶有描述統計的含量特征數據,還具有顯著的地理位置特征信息[23-24]。 根據地理學第一定律[25],地理事物之間距離越近,表現出的空間相關性越強,而特異性地理數據則呈現出弱的地理相關性并在空間上表現出非關聯性。 所以,在立足于決策樹集成學習的隨機森林中加入空間關聯函數將空間位置信息引入,用以探尋樣本數據在連續地理空間上的空間關聯度,并完成對空間特異數據的識別和最優訓練集的構建,減少樣本量和提升計算效率,在區域尺度下的SSC狀態識別與分布規律監測中具有一定的潛力。

濰北平原位于萊州灣南岸地區,土壤的母質及發育過程兼具海洋和陸地的特征,因此在此特定景觀環境下構建準確高效的鹽分估算制圖模型,一直是資源環境領域的重點和難點工作。 本文選取濰北平原為研究區,野外采集233 個土壤樣本并獲取同時相的Sentinel-2 多光譜影像,分別構建兩波段差值、比值、歸一值指數、相鄰三波段正切指數和單波段共135 個光譜參量;進一步采用隨機森林變量重要度評估技術,選擇SSC 敏感光譜參量為輸入自變量,測試得到的鹽分含量值為因變量,分別建立基于隨機森林和空間關聯隨機森林算法的遙感估算;最后完成區域尺度下的SSC 含量估算和數字制圖,以期為實現土地資源優化管理和土壤環境政策制定提供理論支持和方法依據。

1 研究方法

1.1 研究區概況

研究區位于萊州灣南岸濰北平原地區(圖1),地理坐標為北緯37°8'24″～37°16'36″和東經118°37'40″～118°46'24″,占地總面積約為59.14 km2。 作為典型的溫帶大陸性季風氣候區,該地受海洋環境影響明顯且干濕季分明。 據歷史氣象資料記錄,多年平均氣溫和降水量分別為12.7 ℃和608.5 mm。 研究區內整體地形平坦,以平原為主,適合農業種植,主要種植作物類型為玉米和小麥。 地質構造上屬遼冀臺向斜第四系沉積層,當前地質狀況較為穩定。 因距海較近,地下水受到海咸水入侵,在冬季強烈的蒸發作用下土壤表現為輕微的鹽堿化(1 g/kg ＜SSC ＜ 2 g/kg)[26]。

圖1 研究區及采樣點示意圖Fig.1 The study area and samplings sites

1.2 樣點數據與測試分析

在室內首先以ArcGIS 10.2 軟件的數字底圖為基礎,綜合考慮土地利用、水文地質和道路交通可達性等因素,完成233 處土壤樣點的預設。 采樣過程中,將預設樣點周邊10 m 內設定為采樣范圍,采用多點混合法將土壤綜合至1 kg 左右;其后裝入聚乙烯密封袋中,送往實驗室待測。 同時,采用手持式GPS 確定并記錄采樣點的真實地理坐標,采樣過程在2019 年1 月完成。 在實驗室內,首先去除土壤中小石塊、木棒和草根等明顯異質物;接下來,土壤樣品在室溫條件(25 ℃)下風干并過1 mm 篩,完成待測前處理;最后,運用質量法測定土壤樣品的全鹽含量,即吸取一定量基質水浸出液,并蒸干除去有機質,再烘干稱量測得全鹽量[27]。

1.3 影像獲取及預處理

研究所使用的Sentinel-2(哨兵二號)多光譜影像數據(2019 年1 月17 日)由美國地質調查局網站提供,數據獲取條件為晴朗無云天氣及地面無明顯積雪。 下載得到的影像數據為L1C 級的大氣表觀反射率產品,在鹽分估算中,要將L1C級數據轉換為L2A 的地表真實反射率產品。 需對數據進行輻射校正和大氣校正,具體實現過程由歐空局推薦的SNAP 軟件(已安裝Sen2Cor 大氣校正插件)完成[28]。 哨兵二號影像數據中的波段1、 波段9 和波段10 設計用于大氣和水分特征的檢測,不能用于土壤鹽分含量的反演,因此在接下來的計算過程中將其排除。 為確保各個光譜波段與地面信息間的匹配性,所有波段均進行重采樣操作(空間分辨率為10 m)。 并且影像與地面采樣數據進行地理校正,確保像元偏差不大于0.5 個像元。 此外,在土地覆蓋以裸地為主,部分存在植被的區域,參照裸土的判別標準(NDVI ＜0.2)[29-30],剔除存在光譜干擾的植被區域。 地面樣品和影像數據獲取時間正值冬季,降水較為稀少,因此土壤水分對光譜反射率的影響最小。

1.4 光譜指數構建

在土壤有機質、鐵和黏土的重疊吸收作用下,土壤鹽分在可見光-近紅外范圍內的光譜吸收特征通常具有非特異性,表現出覆蓋范圍廣、光譜信號較弱的特點[31]。 因此,除對多光譜影像的10個波段進行光譜分析以外,本研究還選擇光譜差值指數(D)、光譜比值指數(RI)、光譜歸一化指數(ND)和相鄰三波段正切指數(tan)共4 種方式用于挖掘土壤鹽分的特征光譜信號,計算過程:

式中:Bi和Bj分別表示哨兵數據的第i、j 波段的光譜反射率;tanpmq為相鄰的3 個波段(波段p、波段m 和波段q)之間夾角的正切值,用于表征光譜曲線形狀對土壤鹽分的響應能力;kpm是波段p 和波段m 之間連線的斜率;kmq是波段m 和波段q之間連線的斜率。 波段運算在ENVI 5.3.1 軟件中的IDL 二次開發平臺中完成。

1.5 空間關聯隨機森林模型

空間關聯隨機森林模型的構建是在隨機森林算法[32]的基礎之上(流程見圖2),在樣本數據輸入側引入空間權重函數[式(5)～式(7)],用以評估樣本數據間的關聯度和空間特異性,在完成輸入樣本數據的最優組合聚類后,進行回歸模型的擬合計算。 模型計算過程如圖2 所示。

圖2 空間關聯隨機森林模型流程圖Fig.2 Work flowcharts for the spatial random forest model

式中:F 表示因變量yi、自變量xi和位置信息(μi,νi)間的函數關系; Fi(φ) 為對應位置i 處的空間權重高斯函數; δ 和L 分別代表隨機誤差和土壤鹽分真實值;t 表示樣本數據的空間關聯度值,數值越小,樣本數據空間關聯性越強,依據此數據實現對輸入數據集的組合聚類。

決策樹數量(ntree)、總特征集上的特征子集(m try)、變量重要度值(VIS)和地理要素信息協同變量(Xi)是空間關聯隨機森林模型構建過程中的4 項重要參數。 隨機森林算法的核心思想是構建若干棵決策樹(ntree),從決策樹上的總特征集中,隨機選擇m 個特征子集(m try)用作預測變量,最后投票匯總結果;VIS 的大小用于衡量輸入自變量對因變量的影響強度和解釋能力的高低,在模型構建中將根據VIS 的大小選擇最佳的光譜參數用于SSC 估算;此外,為增加空間權重評估函數對樣本數據的位置信息解釋能力和穩定性,選擇典型且易于計算的5 個環境要素特征作為輔助信息,分別是地面高程、植被歸一化指數、土壤濕度指數、地表溫度和離海距離。 地面高程數據免費下載自地理數據共享平臺(http:/ /www.gscloud.cn/),NDVI、土壤濕度指數和地表溫度數據的求算參照文獻[33-34],獲取時間與地面采樣時間相同,樣點離海岸線邊界的距離的計算在ArcGIS 10.2 軟件中的歐式距離計算模塊中完成。在本研究中,ntree 和m try 分別設置為1 000 和3,模型的實現與計算在Python 3.8 中進行。

1.6 精度評估

采集到的土壤樣本總集隨機劃分為建模集(n=209)和驗證集(n=24),建模樣本集用于估算模型的訓練,估算模型精度的驗證在驗證集進行。 決定系數(R2)和相對分析誤差(RPD)是評價模型精度的重要統計指標。 R2表示自變量對因變量的信息解釋能力的強弱,RPD 表征反演模型比只使用均值模型性能提升的能力。 RPD 大于2.0 時,說明模型取得滿意的精度;RPD 處在1.4 ～2.0 之間時表明反演結果尚可接受,有待進一步提升;而RPD 小于1.40 時,意味著模型尚不具備反演能力[17]。 通常,R2和RPD 越大,表明估算模型的精度和穩定性越強。

2 結果與討論

2.1 土壤鹽分含量描述性統計特征

濰北平原表層土壤(0 ～20 cm)SSC 的描述性統計特征如表1 所示。 從表1 可以看出,土壤樣品總集的SCC 為0.49 ～21.22 g/kg 內,均值為2.29 g/kg,覆蓋范圍較大且分布不均勻。 SSC 較高的采樣點主要分布在中部和北部鹽田附近,其余大部分采樣點位于耕地,SSC 較低。 在本研究中,變異系數主要用來評估數據的離散程度[35],3個樣品數據集的變異系數均大于0.36,已處于高度變異,表明區域內SSC 分異較大。 SSC 數據存在空間變異會增加地表反演模型的訓練難度,但一定程度的數據變異會有利于模型收斂并對模型的計算效率產生積極影響,所以隨機抽取到的建模集和驗證集用于土壤鹽分的遙感建模反演是可靠的[36]。 峰度和偏度用于描述數據分布的集中和雙尾特性[37],3 個數據集的偏度值均大于2,表現出正偏,表明外部活動已經干擾了成土之初SSC 正態分布狀態。 土壤樣本在總集合上的峰度偏高,說明鹽分含量分布呈現出非集中的均勻分布態勢。 SSC 的統計特征變異性,一定程度上也反映了土壤鹽堿化的狀態趨勢。 濰北平原地區工農業生產發展造成地下水過度開采,引發海咸水倒灌,加之旱季蒸發強烈,可溶性鹽運移至表層土壤中產生結晶富集;土壤表現出輕微的鹽堿化,并帶來土地肥力下降等一系列的環境問題,這也是SSC 的數據分布表現顯著變異的潛在驅動因素[38]。

表1 土壤SSC 描述性統計特征Table 1 Statistical characteristics of soil salt contents in the surface soils

2.2 選擇光譜特征作為模型輸入

通過式(1)至式(4)共計算得到125 個光譜指數,包括39 個D、39 個RI、39 個ND 和8 個tan。 對SSC 最敏感的前10 個單波段和光譜指數的變量重要度(VIS)如圖3 所示。 B11 是最重要的光譜波段,VIS 達到35.1%,其次是B3 和B8,其余波段的VIS 均小于10%。

圖3 光譜參數變量重要度Fig.3 Variable importance scores of the spectral parameters

通過兩兩波段變換得到的光譜指數,對SSC信息的響應能力提升明顯,波段比值指數的光譜信息提升能力最強;例如,B3 和B4 的變量重要度(VIS)分別為11.6%和4.4%,但經過波段比值變換后,RI34的VIS 達到27.5%,有效增強了土壤鹽分的光譜特征。 B11 是重要度最高的單波段,在RF 變量重要度評估算法輸出中的10 個最優光譜指數中,B11 參與了6 個光譜指數的構建,在鹽分遙感反演模型和光譜指數的構建中應給予重點關注。 基于此,選擇3 個波段(B11、B3 和B8)和4個光譜指數(RI34、RI711、ND611和D45)作為SSC 的特征光譜參數,并在SSC 估算模型構建中充當輸入自變量。

土壤鹽化通常是鹽分離子在蒸發作用下發生運移,在土壤表層生成白色晶體聚合物和鹽殼的過程。 鹽分分子及其基團在晶格結構上振動的合頻與倍頻,會產生特異性的吸收特征峰[39]。WANG 等[3]和SIDIKE 等[40]在對SSC 的光譜吸收特征的室內控制實驗分析時發現,波長在560 nm 和1 640 nm 附近波段與SSC 具有顯著的相關性。 FOURATI 等[41]在利用衛星影像構建估算模型對土壤鹽化水平進行識別時,得出了同樣的結果,即處在短波紅外范圍內的2 個波段(波長在1 600 nm 和2 200 nm 附近)對SSC 表現出極顯著的相關性。 而本研究篩選得到的特征波段B3 和B11 的中心波長分別處在560 nm 和1 610 nm 附近,證實了本研究所選波段的典型性和估算模型構建的可靠性。 此外,張俊華等[8]在分析寧夏銀北地區鹽化土壤的原位光譜分析時指出,600 nm附近是SSC 的一個明顯吸收峰,并且640 ～870 nm光譜區間內波段的光譜反射率與碳酸根鹽分離子顯著相關,這與哨兵影像的波段8 的中心波長范圍高度一致。 總之,SSC 估算模型的輸入特征波段有著明確的物理意義,能夠反映地表土壤鹽分的光譜反射率特征。

2.3 估算模型構建與精度對比

不同的特征光譜參數(特征波段和最優光譜指數) 構建的估算模型的預測結果如圖4所示。

圖4 不同估算模型的土壤鹽分含量估算值與實測值對比散點圖Fig.4 The scatter plots of estimated SSC values versus measured SSC values based on the different model

圖4(a)和圖4(b)分別是以3 個特征波段(B11 、B3 和B8)和4 個光譜指數(RI34、RI711、ND611和D45)作為輸入自變量,測得的土壤SSC作為因變量,訓練得到的基于RF 算法的反演模型的結果。 可以看出,2 個估算模型的驗證精度RPD 均小于1.40,實測值和估算值偏離1 ∶1 線較大,表明構建模型尚不能夠進行土壤鹽分的準確估算。 圖4(d)和圖4(e)表示的估算模型是以空間關聯隨機森林算法為基礎構建,對比RF 算法構型的模型[圖4(a)和圖4(b)],以空間關聯隨機森林算法為基礎構建的估算模型[圖4(d)和圖4(e)] 的驗證精度RPD 分別提升至1.17 和1.60,以4 個光譜指數(RI34、RI711、ND611和D45)為輸入變量的空間關聯隨機森林模型,基本能夠較準確地實現土壤鹽分的反演[RPD ＞1.40,圖4(e)]。

圖4(f)是組合3 個特征波段和4 個光譜指數作為輸入自變量,得到的空間關聯隨機森林模型估算效果,土壤鹽分實測值與估算值的對比的散點基本在1 ∶1 線兩側,估算值與實測值偏差較小,反演結果最佳(R2=0.89 和RPD=2.04);而采用相同輸入變量的RF 估算模型[圖4(c)],則尚不能準確進行SSC 的遙感估算。 總之,僅用特征波段或最佳光譜指數來構建SSC 估算模型均不能取得滿意精度,空間關聯隨機森林模型的精度和穩定度要高于隨機森林模型,再將特征波段和最佳光譜指數組合,輸入到空間關聯隨機森林算法中去,建立得到的反演模型能夠較準確地完成土壤鹽分的遙感估算。

2.4 土壤鹽分空間分布

以特征波段和光譜指數作為輸入自變量的空間關聯隨機森林模型在不同空間關聯度(t)下的精度變化如圖5 所示。 當t 為4 時,共有197 個數據樣本(n 表示輸入樣本數)參與估算模型的建立,R2為0.77;t 繼續減少,選取得到的數據樣本間的空間關聯度增強,t 為2 時,共有185 個土壤樣本參與模型構建,取得了最高的估算精度(R2=0.89);隨著樣本空間關聯度要求提升,樣本量減少,模型逐漸不夠完全訓練,精度隨之下降。

圖5 不同t 時空間關聯隨機森林模型精度Fig.5 The perform ance changes of spatial random forest model with the different t value

為進一步驗證構建的土壤鹽分反演模型的穩定性和可靠度,將構建得到的最優反演模型推廣至整個研究區,計算得到的鹽分空間分布如圖6所示。

圖6 土壤鹽分含量分布圖Fig.6 The spatial distribution of the soil salt contents

從圖6 可以看出,土壤鹽分含量在中部地區略高于北部和南部地區,鹽分含量高值區與鹽田的分布位置密切相關。 此外,研究區土地利用主要為農用耕地,地質地貌類型也較為單一,土壤鹽分含量的變異并沒有在大的空間尺度上表現得很明顯。 但是,SSC 在田塊尺度上表現出顯著變化,這與當地的微地貌和海水倒灌密切相關。 研究區是國內知名的蔬菜生產基地,生產蔬菜的耕地需要進行大量抽水灌溉,地下水儲備不足時引發海水倒灌,可溶性鹽分運移至表層土壤產生一定程度富集,所以不同田塊間的SSC 存在明顯差異。

在遙感估算中,地理環境過程時刻在變化,土壤鹽分的光譜特征僅強調瞬時信號是遠遠不夠的。 下一步要將遙感反演的經驗統計模型與土壤鹽分動態和環境的交互作用結合起來,從而最大化地減弱水汽、粗糙度和植被殘積物對反演的影響,最終建立起穩健的遙感反演模型,實現動態大尺度層面上的土壤屬性監測。 其次,計算得到的光譜參數的重要度存在一定差異(即光譜參數對SSC 的響應能力不同),但在模型輸入中并未考慮不同光譜參數的權重,若能將其對土壤鹽分的不同光譜參數的響應強度考慮進估算模型,將有助于選取最佳分裂節點和決策樹數目,提升模型的運算精度和效率。 最后,混合像元問題是遙感估算中難以忽略的問題[12,42]。 混合像元的存在,會減弱土壤鹽分的目標光譜特征,增大信號噪聲和模型誤差。 因此,在研究中選用空間分辨率(10 m)較高的Sentinel-2 影像作為數據輸入源,土壤樣點的采集也均在純像元中進行,并且影像處理中掩膜掉道路、植被和建設用地等非裸土區域,以此最小化混合像元對反演結果的影像;并利用光譜指數構建法凸顯SSC 的光譜信息,同時減弱環境弱信號。 但是,實現像元解混合采用超高空間分辨率影像可能是解決混合像元問題的最佳方案[43-44],在下一步研究中需要參考改進。

3 結論

1)哨兵影像的B3、B8 和B11 是SSC 的特征光譜波段,B11 對SSC 的重要度值最高。 兩兩波段間比值變換能夠有效綜合多個波段的光譜信息,增強衛星光譜信號對SSC 的響應,最優的光譜指數分別為RI34、RI711、ND611和D45。

2)空間關聯隨機森林模型的精度和計算效率優于隨機森林遙感估算模型,以3 個特征波段和4 個最優光譜指數作為輸入自變量時,空間關聯隨機森林算法支持下的遙感估算模型的精度指標R2和RPD 分別達到0.89 和2.04,能夠較準確地完成SSC 的反演制圖。

3)SSC 的空間分布整體上中部略高于南部和北部,高值區主要受中部地區的鹽田分布影響;農業活動和微地貌使SSC 的分布在田塊尺度上分異趨勢顯著,建議采用農業措施、水利措施、生物措施相結合的方式進行土壤鹽化治理。