多尺度土壤質地與光譜空間非平穩性關系探究

鄭崳珍，陳奕云，陳 敏，吳子豪，蔣江俊男

（1.武漢大學資源與環境科學學院/自然資源部數字制圖與國土信息應用重點實驗室/地理信息系統教育部重點實驗室，武漢 430079；2.廣州市城市規劃勘測設計研究院，廣州 510060）

土壤與人類的生存和發展有著密切聯系，在解決全球現實問題，尤其是糧食安全、水安全、氣候變化、生物多樣性等方面，土壤居于中心地位［1］。在全球可持續發展目標（Sustainable development goals，SDGs）中，有多個目標與土壤密切相關，人類社會已經認識到，科學合理地利用和管理土壤信息對人類自身的生存和可持續發展至關重要［2］。土壤質地是土壤重要的物理屬性，也是陸面過程模擬的關鍵參數，對土壤結構、孔隙狀況、保肥性、保水性、耕性等均有顯著影響［3］。有效地獲取、存儲、表達、傳輸與分析土壤質地信息，理解和表征土壤質地的空間變化信息，能為糧食安全、生態文明建設、鄉村振興、精準扶貧等國計民生提供重要決策支持［4，5］。

傳統土壤質地信息獲取具有周期長、成本高、過程復雜等顯著性缺點，難以進行大范圍、高覆蓋度的重復調查。基于高光譜技術，可以快速獲取土壤屬性性質及其空間分布特征。高光譜在土壤研究中的應用歷史可以追溯到19世紀20—30年代［6］。隨后，人們開始利用反射光譜研究土壤含水量、土壤組分、粒徑及有機質含量。自20世紀60年代起，研究者基于可見光-近紅外（Vis-NIR）光譜較為成功地反演了一系列土壤理化屬性。如今高光譜技術被廣泛地用于反演土壤中有機質含量、含水量、重金屬及土壤質地等方面的研究［7-13］。

目前，高光譜技術在土壤質地領域研究熱點多集中在光譜預測模型方面。學者們大多采用多元線性逐步回歸法、偏最小二乘回歸法、隨機森林算法、BP神經網絡、支持向量機等一系列線性和非線性方法進行土壤質地高光譜預測。張娜等［14］利用一元線性回歸、多元逐步回歸及BP神經網絡3種方法建立土壤光譜反射率與土壤沙粒、粉粒之間的預測模型并對比了預測精度。王德彩等［15］應用Vis-NIR光譜結合BP神經網絡預測土壤質地組分顆粒黏粒和沙粒，結果顯示BP人工神經網絡建立的模型精度更優。Hobley等［16］結合偏最小二乘回歸模型和隨機森林模型預測土壤質地，發現偏最小二乘回歸（PLSR）模型對黏粒和沙粒預測效果最好，而RF模型對粉粒組分預測效果較好。然而，這些方法忽略了土壤質地與高光譜關系的空間異質性，使得殘差存在空間結構。

地理加權回歸（GWR）是一種常用于估計目標變量的空間分布和探究解釋變量和響應變量之間非平穩關系的空間分析方法。由于GWR考慮了關系的空間異質性，能反映局部情況且具有更高準確性，已被廣泛用于土壤屬性空間預測［17］。在此基礎上，Fotheringham等［18］于2017年提出了多尺度地理加權回歸（MGWR）方法。MGWR不僅繼承了GWR反映空間對象局部效應的優點，還考慮了不同影響因子的差異化作用尺度，能更靈活地識別不同空間尺度下各自變量影響因變量的地理過程［19］。沈體雁等［20］應用MGWR識別出不同尺度下北京市2011—2017年二手住宅交易價格的空間分異特征。MGWR模型還能應用于探究土壤屬性和光譜關系空間非平穩性的問題。喬磊等［21］應用MGWR方法，很好地實現土壤有機質的空間預測，得到各環境協變量對有機質影響效應的空間變化，這種影響效應在不同環境協變量中有著不一樣的空間尺度。因此，使用MGWR可以更好地探究不同空間尺度下光譜與土壤質地之間的非平穩性關系。

云南省杞麓湖是典型的高原湖泊，湖濱平原大部分被開墾為農田，農田分布細碎化特征導致農民土地管理方式多樣化，進而使土壤質地分布可能具有異質性特征［22］。針對受人類活動影響較大的農耕區，土壤質地和土壤光譜關系是否具有空間非平穩特征更值得被探索。因此，選取云南省杞麓湖流域農田為研究區，以光譜的不同潛變量作為協變量，土壤質地的3種顆粒含量作為因變量，對比PLSR、GWR、MGWR 3種方法下土壤質地擬合模型效果，并應用MGWR探究多尺度效應影響下土壤質地和高光譜關系的空間非平穩特征，有益于加深對流域土壤屬性空間分布的認識和理解，以期為相似流域開展空間異質性研究提供普適性參考。

1 材料與方法

1.1 研究區概況

杞麓湖流域（102°33′48″—102°52′36″E，24°4′36″—24°14′2″N）位于云南省中部，隸屬云南省玉溪市通海縣，南北銜接曲江干流和江川星云湖，東西鄰近華寧龍洞河與玉溪大河（曲江上游段）。該流域屬亞熱帶季風氣候，干濕季分明，多年平均降水量約899 mm，年平均氣溫約15.6℃，年日照時間約2 300 h。流域地形為向南突出的新月形斷陷盆地，平均海拔約1 800 m。湖盆土壤類型主要有山地紅壤、紫色土和水稻土等。沿湖平原是重要的農耕區，土地利用方式主要是農業用地，水稻、小麥、烤煙、油料是主要農作物，蔬菜產業尤為發達。流域農田較為細碎，分布零星，70%以上現有耕地田塊面積小。由于受耕作、施肥等人類活動影響，土壤質地與光譜之間的關系可能具有非平穩性。210個土壤樣本均采集自杞麓湖流域周邊農田（圖1）。

圖1 土壤樣點空間分布

1.2 土壤機械組成測定和光譜預處理

采集土壤分成2份，分別用于土壤質地室內分析和室內光譜測量。樣本測定前先經自然風干、過篩、研磨等操作后，過2 mm篩孔，經H2O2-HCl-（NaPO3）6法預處理后，形成最終上機樣本。土壤質地測定方法采用激光粒度儀法，試驗儀器為馬爾文激光粒度儀（Mastersizer 3000），量程達0.01～3 500.00μm。土壤顆粒采用美國農業部土壤質地分類標準（USDA）進行分級：黏粒（＜0.002 mm）、粉粒（0.002～0.050 mm）和沙粒（0.050～2.000 mm）。

研磨土樣過20目篩孔后進行室內高光譜數據測量，試驗條件為避光室內環境。使用ASDFieldSpec 3地物光譜儀測定反射光譜，獲取的地物波長范圍為350～2 500 nm，輸出波段數可達2 151 nm。每次測量前及時進行標準白板校正，以保證測量效果穩定。光譜分析前去除邊緣噪聲，保留信噪比較高的400～2 350 nm波段。光譜預處理步驟包含Savitzky-Golay平滑（7點平滑數）、一階微分變換、均值中心化等，以上操作均在MATLAB中進行。

1.3 模型建立與評價

1.3.1 偏最小二乘回歸（PLSR） 偏最小二乘回歸由Wold等［23］首次提出，最早應用于化學工程領域，目前常用于定量光譜分析領域［16］。PLSR集成了多元線性回歸、主成分回歸和典型相關分析的優點，能有效解決光譜多波段間存在的多重共線性問題。PLSR的計算公式如下。

式中，y代表目標因變量；x代表土壤光譜參數，β（ii=1，2，……n）表示xi的回歸系數。使用PLSR建立多自變量（光譜）與因變量（土壤屬性）之間的回歸關系，PLSR因子即潛在變量（LVs）的最佳數量通常使用交叉驗證法來確定。模型建立采用Matlab 2014軟件和PLS_toolbox工具箱實現。

1.3.2 偏最小二乘-地理加權回歸（PLS-GWR） 地理加權回歸是探索潛在空間非平穩關系最為廣泛應用的方法之一［24-26］。GWR模型是對普通線性回歸模型的擴展，將數據的地理位置嵌入到回歸參數中。與普通線性回歸相比，GWR原理簡單、操作性強，還可提高模型擬合度和降低殘差空間自相關性。與傳統全局回歸模型相比，GWR不會對模型中的每個關系產生“平均”全局估計，而是允許自變量和因變量之間的關系因空間而異。PLS-GWR即偏最小二乘法與地理加權回歸法的結合，是將偏最小二乘方法中獲取的多個潛變量作為地理加權回歸的自變量，考慮了變量的空間位置，從而建立起自變量和因變量之間的回歸關系。

假設共有n個觀測樣點，每個觀測點i∈（1，2，……n）所屬的空間位置為(ui,vi)，GWR的線性回歸方程如下。

式中，(ui,vi)指第i個樣本的空間坐標；xij指i位置的第j個自變量；yi是i位置上的因變量；βj(ui,vi)是i位置(ui,vi)上的第j個自變量的回歸系數；εi是隨機誤差項。

空間權重矩陣和帶寬的確定是地理加權回歸建模計算和分析的關鍵，GWR分析中最常使用高斯函數（Gaussian）和雙平方函數（Bisquare）2種核函數確定空間權重矩陣。帶寬的確定方法仍以交叉驗證準則（CV）或修正赤池信息量準則（AICc）為主。模型各項參數選擇如下：建模采用雙平方自適應核函數，帶寬的選取采用默認黃金分割搜索法，以修正赤池信息量準則（AICc）作為信息評價準則。經典GWR模型中預設所有建模的過程皆在相同的空間尺度中操作，即所有自變量共用相同的有效帶寬和空間權重矩陣。

1.3.3 偏最小二乘-多尺度地理加權回歸（PLS-MG?WR） 通過放寬模型中所有空間變化過程在相同空間尺度上運作的假設，Fotheringham等［18］提出多重尺度地理加權回歸（MGWR），這是在經典GWR基礎上優化的一種改善有效帶寬選擇的地理空間回歸方法。不同于GWR模型中各自變量共用一個帶寬，MGWR模型在擬合過程中針對每個自變量選擇了差異性帶寬，這意味著每個空間位置上針對各變量都有特定的空間權重矩陣。PLS-MGWR是偏最小二乘方法與多尺度地理加權回歸方法的結合，是將偏最小二乘方法中獲取的多個潛變量作為多尺度地理加權回歸的自變量，可以看作對PLS-GWR的改進。MGWR模型的表達式如下。

式中，βbwj指第j個自變量的回歸系數；bwj是計算第j個自變量回歸系數時使用的帶寬；εi是隨機誤差項。

MGWR和經典GWR選擇的核函數和帶寬選擇準則相同，分別為自適應高斯函數和修正赤池信息量準則（AICc）。MGWR模型預設所有建模的過程皆在不同的空間尺度中操作，即針對每個自變量確定差異性帶寬和空間權重矩陣。后向擬合算法（Back-fitting）是校準MGWR模型的核心算法，其基本思想是模型中每一個加性項都使用單個光滑函數來估計，而每一個加性項都可以解釋因變量和自變量之間的因變關系。遵循廣義加性模型（GAM）的邏輯，MGWR模型還可表示如下。

式中，MGWR中的fj是指第j個加性項fj；ε是殘差項。

多尺度地理加權回歸主要基于線性可加模型和后向擬合法進行參數估計迭代，其主要做法是以平滑因子或者局部變量參數的地理加權回歸估計值作為初始值，結合后向擬合法和迭代過程，來實現對平滑因子或變量參數的估計［18］（圖2）。以經典GWR各參數和結果作為MGWR初始值，得到初始化殘差ε（y的實測值與初始估計值之差）。將初始化殘差ε與第一加性項f0相加得到新值，對第一個自變量x0和新值進行經典GWR回歸，得到有效帶寬bw0和新殘差，以此更新下一個加性項f（jβbwjx）j，隨后新殘差加第二個加性項f1與變量x1再進行GWR回歸，得到新帶寬bw1和新殘差，重復以上過程，直到與最后一個變量xm相關的局部參數都已估算完成即可。

圖2 MGWR中后向擬合算法原理

1.3.4 模型評價指標 根據土壤反射光譜和黏粒、粉粒、沙粒3種顆粒之間的關系，建立PLSR、PLS-GWR、PLS-MGWR模型，隨后進行模型質量評價。評價采用以下指標，包含決定系數（R2），平均絕對誤差（MAE），均方根誤差（RMSE），全局莫蘭指數（Moran′sI）。

式中，n是參與建模的土壤樣本數量是模型中黏粒（粉粒或沙粒）的擬合值；yi是實驗室黏粒（粉粒或沙粒）測量值是測量值均值。決定系數越大，平均絕對誤差和均方根誤差越小，表明模型相對越穩定，相對擬合程度越好。莫蘭指數介于［-1，1］，絕對值越大則表明空間自相關性越強；反之，空間自相關性越弱。使用Moran’sI來檢測模型殘差的空間結構［5］，若結果顯著，則說明殘差具有空間結構，模型的參數估計結果不可靠，反之，則說明殘差空間分布隨機，模型結果可靠。

2 結果與分析

2.1 杞麓湖流域農田土壤質地空間分布狀況

本研究土壤質地分類標準為美國制，其中土壤質地等邊三角形三頂點分別代表沙粒、黏粒、粉粒的0或100%的含量。杞麓湖土壤樣本的質地分類結果如圖3所示。

圖3 土壤質地三角形

由圖3可知，杞麓湖流域農田土壤樣本質地類型較為齊全，共有10類，其中以粉質壤土、粉沙質黏壤土、粉沙質黏土為主，少量為沙土、壤質沙土、粉沙土、壤土、沙質壤土、黏壤土、沙質黏壤土。土壤中黏粒、粉粒、沙粒含量波動范圍均較大，分別在0～57.51%、12.80～91.03%、0～87.20%波動（表1）。3種顆粒含量中粉粒均值最高，達68.09%，黏粒次之，沙粒最小。這一粒徑分布特征很可能與當地高強度的土地利用有關。

從空間分布（圖4）來看，黏粒含量大致呈現東高西低、南高北低的空間格局。其中，黏粒含量的高值集中在流域的東南部，這與杞麓湖流域南部土壤黏性較大的實際調研情況一致。粉粒含量在流域東西方向上分布較為平緩，在南北方向上呈現中部高、兩側低的趨勢。其中，粉粒的高值集中在杞麓湖南側。沙粒含量總體偏低，大致呈現西高東低、北高南低的空間格局。其中，沙粒的最高值出現在杞麓湖南側，根據實地調研發現由于當地農民為了改良土壤性質，在比較黏重的土壤中摻沙，以此保證土壤通透性。

表1 土壤質地的特征統計值

圖4 土壤質地的空間趨勢

2.2 模型比較

分別以土壤的黏粒、粉粒、沙粒為因變量，波段為400～2 350 nm的反射光譜為自變量，根據式（1）至式（7），運用PLSR、PLS-GWR、PLS-MGWR方法分別建立3種顆粒含量與反射光譜之間的關系，得到9種擬合模型。其中，PLSR模型中的最佳潛變量數（LVs）使用舍一交叉驗證法確定。經擬合，黏粒、粉粒、沙粒3種顆粒對應的LVs分別為5、2和3，利用評價指標對9種模型擬合效果進行評價，如表2所示。與PLSR擬合的土壤質地模型相比，MGWR和GWR各模型精度均有所提升，這是由于這2種方法均考慮了光譜在不同空間位置對土壤質地響應的差異。從整體上來看，MGWR對土壤質地（黏粒、粉粒、沙粒）的擬合效果（R2）和穩定性指標均優于PLSR和GWR，這是因為MGWR不僅考慮了土壤光譜與質地關系的空間異質性，還針對各光譜潛變量選擇了不同空間尺度，對每個自變量的回歸參數都進行局部估計，從而提升了模型精度。此外，PLSR方法建立的粉粒和沙粒模型殘差的莫蘭指數統計顯著，表明殘差具有空間依賴性，即存在空間結構，違背了經典統計方法中的殘差獨立性假設，使得模型精度較低；而GWR和MGWR擬合的各模型殘差的莫蘭指數在0.005水平上均不顯著，說明GWR和MGWR模型殘差為空間隨機，其擬合結果更為可靠。綜上所述，MGWR是一種適于擬合光譜與土壤質地關系的潛力模型。

表2 多尺度地理加權回歸、地理加權回歸和偏最小二乘回歸模型建立效果比較

2.3 MGWR標準化回歸系數空間格局分析

由于杞麓湖流域田塊較為分散，光譜的不同潛變量與土壤質地的關系在實際情況中可能隨著地理位置的變化而變化，呈現出空間非平穩性。多模型比較結果表明，MGWR方法針對實際情況的擬合效果較好，因此，選擇MGWR方法對黏粒、粉粒、沙粒分別進行建模，進一步探究每種模型中光譜各潛變量與3種土壤顆粒含量之間的關系。MGWR各系數統計描述見表3。多尺度地理加權回歸系數的空間格局如圖5所示。

表3 多尺度地理加權回歸系數統計描述

2.3.1 黏粒模型 潛變量LV1至LV5整體均具備正向響應。

1）LV1系數反映了光譜潛變量LV1對黏粒含量變化的響應程度。LV1回歸系數的取值在0.361～0.387，均值為0.373，標準差為0.004。高值主要集中在杞麓湖的北部區域，這說明在杞麓湖北部，隨著黏粒含量增加，潛變量LV1增大；低值出現在杞麓湖流域的最西部區域，在此處，黏粒含量變化較小，潛變量LV1對黏粒的響應較弱。從系數的絕對值來說，LV1在5個變量中對黏粒含量的響應程度最大。

2）LV2對黏粒含量產生的是正向的響應，在空間上呈現出一定的分層結構。高值集中在杞麓湖流域西北部，低值集中在東南部，從流域西北方向至東南方向，潛變量LV2對黏粒含量響應程度逐漸降低。回歸系數的取值在0.253～0.304，均值為0.265，標準差為0.013。從系數的絕對值上來看，其響應程度為中等程度。

3）對于潛變量LV3，同樣對黏粒含量具有正向響應，黏粒含量越大，潛變量LV3顯示出越大的趨勢。在空間分布上，高值和低值區域相銜接，分別集中在流域的北部和西北部區域。LV3回歸系數的取值在0.091～0.432，均值為0.213，標準差為0.070。從系數的絕對值上來看，潛變量LV3的響應強度最小。

4）潛變量LV4因黏粒含量具有正向響應，其系數取值在0.117～0.440，均值為0.277，標準差為0.108，整體呈現由南北兩端向雙側逐漸降低的趨勢，在流域西南部，LV4對黏粒的正相關響應程度達到最小。從系數的絕對值上來說，潛變量LV4的響應強度適中。

5）LV5回歸系數的取值在-0.097～0.611，均值為0.283，標準差為0.195，呈現出南高北低的空間格局。總體上，潛變量LV5對黏粒的響應為正向，正向效應最大值集中在杞麓湖南部區域，但北部有較小區域產生了負效應，這意味著南部區域潛變量LV5隨黏粒含量同向增減，北部負效應區域隨黏粒變化呈現相反趨勢。在5個潛變量中，潛變量LV5的響應強度略高于LV2。

2.3.2 粉粒模型

1）潛變量LV1系數數值變化較大，在-0.196～1.238波動，均值為0.426，標準差為0.335。整體上呈現出西低東高的格局，高值集中東南部，低值集中在西北部。東南部區域LV1的響應程度最強，此處粉粒的變化程度也最大。西部區域LV1表現出負效應，即LV1由于粉粒含量的提升而有所降低，反之亦然。粉粒模型中，LV1的絕對值遠大于LV2，對粉粒的響應程度極高。

2）潛變量LV2正向響應粉粒含量，即粉粒含量與LV2同步增減。LV2系數取值范圍在0.040～0.244，均值為0.180，標準差為0.073。在空間分布上，LV2系數在流域北、南、東3個方向數值偏高，西部地區普遍偏低。在最西部區域LV2對粉粒的響應程度最弱，南部區域響應程度最強。從系數的絕對值上來看，LV2對粉粒的響應程度偏弱。

2.3.3 沙粒模型 潛變量LV1、LV2、LV3對沙粒響應程度各不相同，總體上均為正效應。

1）LV1正向響應沙粒含量，沙粒含量越高，LV1越高。LV1回歸系數的取值在0.051～1.698，均值為0.422，標準差為0.323。空間上，系數高值低值交錯分布，南部區域沙粒的變化幅度達到最大，LV1對沙粒含量響應程度達到最大。從系數的絕對值上來講，LV1對沙粒響應程度為3個響應潛變量中最大。

圖5 多尺度地理加權回歸系數的空間格局

2）LV2回歸系數的取值在-0.245～0.713，均值為0.194，標準差為0.233。整體上高低值分布較為分散，局部為負，大半為正。負值主要集中在杞麓湖東南區域，沙粒含量隨著LV2的增加而降低。從系數的絕對值上來講，LV2的響應程度較低。

3）LV3系數空間分布上高低值區域參半，低值零散，高值集中在流域北部。LV3回歸系數的取值在-0.110～1.408，均值為0.199，標準差為0.243。負值集中在遠離湖周的西部，此處沙粒含量對LV3有負向影響，LV3隨沙粒的降低而增加。LV3對沙粒的響應程度與LV2相似，都為較弱程度。

3 小結與討論

本研究使用210個杞麓湖流域的實測土壤質地數據，揭示土壤質地（黏粒、粉粒、沙粒）的空間格局，并對比了偏最小二乘回歸、地理加權回歸和多尺度地理加權回歸3種不同方法擬合土壤質地和光譜關系的效果，最后探究差異化尺度作用下兩者關系的空間異質性。

土壤質地描述性統計結果顯示，研究區土壤質地類型以粉質壤土、粉沙質黏壤土、粉沙質黏土3種為主要類型，3種顆粒百分含量均值由大到小排序分別是粉粒（68.09%）、黏粒（22.14%）、沙粒（9.77%）。空間趨勢上，黏粒含量呈現東高西低、南高北低的格局；粉粒含量在南北方向上呈現中部高、兩側低的特點；沙粒含量呈現西高東低、北高南低的趨勢。模型評估結果顯示，MGWR的MAE、RMSE和R2指標均優于PLSR和GWR，且殘差為空間隨機，因此使用MG?WR擬合土壤質地與光譜間關系效果最好。MGWR模型回歸結果表明，光譜各潛變量隨著空間位置變化對土壤質地3種顆粒組分的響應效應顯示出不同的空間格局。從系數絕對值上來講，潛變量LV1對黏粒含量響應程度最大，系數為0.373±0.004；對粉粒含量響應程度最大的潛變量是LV1，系數為0.426±0.335；LV1對沙粒含量的響應程度最大，系數為0.422±0.323。

與GWR相比，MGWR建立了更加有效的空間模型，這主要是由于MGWR考慮了光譜空間上的不同潛變量的不同響應尺度，針對每個變量選擇特定帶寬作為其擬合過程中的空間尺度指標，從而避免更大誤差，使得模型結果更為穩健。本研究考慮了土壤質地與高光譜關系的空間異質性，使用MGWR更好地擬合了不同空間尺度下土壤質地和光譜各潛變量之間的關系，使得空間分析結果更加真實可靠，為開展土壤質地空間異質性研究提供方法借鑒。

由于條件限制，獲取的近地高光譜實測數據樣本數較少且為點狀，并非區域連續的面狀數據，因此需要通過大量樣本來驗證模型的普適性。未來結合機載高光譜、高分五號等遙感影像可以更好地實現大區域尺度的土壤質地空間分析，這將更好地服務于土壤質量風險性評價和生態環境保護。