基于成土環境地理鄰域分析的歷史土壤圖訓練樣本篩選*

高 鴻 朱 娟 王良杰 趙玉國 張甘霖

（1 土壤與農業可持續發展國家重點實驗室（中國科學院南京土壤研究所），南京 210008）

（2 中國科學院大學，北京 100049）

（3 安徽師范大學國土資源與旅游學院，安徽蕪湖 241002）

（4 安徽省測繪局，合肥 230031）

（5 南京林業大學林學院，南京 210037）

20世紀我國進行的兩次土壤普查產生的歷史土壤圖已應用至農業、生態環境等多個領域，成為土壤資源管理、生態水文模型等研究的主要數據來源［1-2］。然而，受當時制圖者經驗水平差異、傳統制圖技術以及數據支撐等的影響與限制，歷史土壤圖中存在邊界錯置的現象，空間分布精度受限，加之隨著時間的推延土壤圖的時效性降低，因而利用數據挖掘模型獲取土壤—環境關系知識，進行土壤圖的更新的研究越來越多［2-5］。此外，土壤圖可以為缺少大量野外樣點的地區提供土壤與環境關系知識［6］，從土壤圖中挖掘土壤—環境關系知識需要對數據進行預處理和篩選［6-8］。而無論是土壤圖更新中數據挖掘模型訓練樣本選取還是土壤圖知識挖掘中對土壤圖數據的預處理，關鍵點和難點均是如何在歷史土壤圖上選取具有代表性的訓練樣本。因而，如何從歷史土壤圖中有效選取代表性的訓練樣本具有較為重要的意義與價值。

土壤圖更新和知識挖掘所需要的訓練樣本選取是指從不同土壤類型中篩選出樣本數量與土壤圖中面積大小相適應的、可代表土壤類型成土環境和圖斑空間分布的訓練樣本，包括樣本數量和樣本空間位置的確定兩個主要內容。對于樣本數量的確定，劉雪琦等［9］確定了一種基于土壤類型面積分級的訓練樣本數量確定的方法，Odgers等［10-11］在每個圖斑里選擇相同數目的訓練樣本。對于樣本空間位置的確定，Odgers等［10-11］利用隨機方法篩選樣本，黃巍等［4］、朱阿興等［6］、Qi和Zhu［7］通過建立環境因子直方圖，在峰值區間進行土壤圖樣本選取，Yang等［3］利用模糊聚類法篩選土壤圖樣本并進行土壤圖更新。

以上研究在確定樣本數量時主要針對的是土壤類型，并未針對土壤類型所具有的多個圖斑設計樣點，這樣可能會造成對較大面積圖斑代表性不夠或不同面積圖斑之間樣點數量失衡，而每個圖斑均具有一定的地理空間意義，因而需要針對每個圖斑設計樣點數量。在確定樣本位置時，通過環境因子直方圖和聚類進行樣本位置的選擇只考慮了環境因子屬性數值大小，忽略了環境因子在地理空間分布上的特征，而進行土壤類型空間預測所使用的樣本應具有空間代表性，對環境因子進行空間分布變化的分析有助于篩選出成土環境穩定、空間代表性高的訓練樣本。基于此，本文提出一種基于環境因子鄰域分析的土壤圖訓練樣本篩選方法，包括根據土壤圖圖斑單元面積分段線性縮放的方法確定樣本數量和基于環境因子地理空間分布特征來確定樣本空間位置。

1 材料與方法

1.1 研究區概況

本文選取安徽省宣城市旌德縣為研究區域（圖1），11 8°1 5′～11 8°4 4′E，30°07′～30°29′N，位于安徽省南部。四面環山，地貌類型包括中山、低山、丘陵和山間盆地，海拔100～1 272 m；氣候屬于北亞熱帶濕潤季風氣候；植被以亞熱帶常綠闊葉林為主，包括常綠闊葉林、落葉闊葉林、馬尾松林、灌叢、草叢等。旌德縣第二次土壤普查歷史土壤圖包含紅壤、黃壤、石灰巖土、潮土和水稻土等5種土類，9種土壤亞類，分別為山地黃壤（編號為S1）、棕色石灰土（編號為S2）、淹育型水稻土（編號為S3）、潛育型水稻土（編號為S4）、潴育型水稻土（編號為S5）、灰潮土（標號為S6）、紅壤性土（編號為S7）、黃壤性土（編號為S8）、黃紅壤（編號為S9），共1 364個土壤類型圖斑。本文基于土壤亞類級別在土壤圖斑單元中篩選訓練樣本。

1.2 數據來源

本文所使用數據包括第二次土壤普查形成的土壤圖，將其數字化成矢量數據格式。此外，由于地形要素基本上可以代表土壤形成與發展過程中的主要影響要素［6］，所以本文所使用的另一主要數據為數字高程模型。高程數據來源于地理空間數據云平臺（http://www.gscloud.cn），柵格空間分辨率為30 m×30 m，在SAGA中計算坡度，采用土壤形成過程中高程和坡度這兩個主要地形環境因子進行分析。本文采用ArcGIS 10.0和python2.7 GDAL進行空間數據處理。

1.3 樣本數量計算

歷史土壤圖存在許多大小不一的圖斑單元，其中面積較小的圖斑單元受制作者主觀影響較大，圖斑邊界及整體的可信度較低，因而本文設置面積閾值，將低于該閾值的圖斑單元剔除。通過設置面積閾值選取歷史土壤圖中面積較大的土壤圖圖斑，再對所選擇的圖斑根據面積大小進行分段，得到不同面積子區間，不同的面積區間對應不同的樣本數量區間，面積較小的區間樣本數量小，反之樣本數量大，在不同面積區間通過線性縮放法建立面積大小與樣本數量的線性映射，計算得出所有圖斑對應的具體樣本數量。當存在過大面積圖斑時，對面積進行分段可消除根據圖斑面積比例直接確定樣點數量造成的樣本數量差異懸殊；在分段后的區間內通過線性縮放法建立面積大小與樣本數量的線性映射，基于面積大小計算得出樣本數量，則會保證樣本數量與整個研究區不同面積圖斑相對應。

圖1 研究區位置Fig. 1 Location map of the study area

對于面積較大的山地黃壤、棕色石灰土、潛育型水稻土、潴育型水稻土、紅壤性土、黃壤性土、黃紅壤等土壤類型，圖斑單元面積閾值為50 hm2，面積較小的淹育型水稻土和灰潮土，圖斑單元面積閾值為10 hm2，面積閾值是結合不同土壤類型圖斑的面積和數目確定的，主要目的是刪除精度較低的細小圖斑和減小計算量，具體閾值可隨不同的研究數據而變化。篩選后圖斑總面積占原始總面積的79.7%，圖斑單元數占原始單元數的18.6%，通過面積閾值篩選保留了較少數量的分析圖斑，但同時分析圖斑的面積占研究區總面積的比重較大，減小了分析過程中的圖斑計算量。

根據面積閾值篩選后的所有圖斑面積分布，設置不同的面積區間和對應的樣本數量區間（表1）。表中圖斑面積分段區間值和對應的樣本數量區間值是作者多次嘗試確定的，具體的參數值可隨不同的研究區域和不同的研究者而異。

表1 面積區間與樣本數量區間對應表Table 1 Reference between polygon area range and sample quantity range

對各面積區間，其中不同面積圖斑單元所對應樣本數量的計算公式為：

式中，Ni表示不同區間第i個圖斑樣本數量；Areai表示不同區間第i個圖斑的面積；Nmax表示樣本數量區間的上限；Nmin表示樣本數量區間的下限；Areamax表示該區間圖斑面積的最大值；Areamin表示該區間圖斑面積的最小值。

1.4 樣本空間位置確定

對于樣本位置，本文假設對于某些主要成土環境因子，土壤圖圖斑訓練樣本的位置位于該因子空間變化穩定的位置。因為一方面當成土環境要素在地理空間某一尺度局部區域均質、穩定的分布時，其對土壤的影響也是持續、穩定的存在，更可能對應著特定的土壤類型；另一方面，從環境因子變化穩定的局部區域選擇位于其地理中心的位置作為樣本，樣本對該局部區域成土環境具有較高的空間代表性。在土壤圖圖斑內，采用較小空間尺度，選擇環境因子空間變化穩定的位置作為該圖斑內土壤類型訓練樣本的位置，篩選出可能是土壤類型所需成土環境并且空間代表性較高的多個空間位置。通過設置指標對環境因子進行柵格窗口鄰域分析，得出當前柵格窗口的尺度下環境因子空間變化程度的分布情況，進而篩選出環境因子空間變化程度最小也即空間變化最穩定的位置。具體思路是，在圖斑單元內，按照環境因子鄰域分析結果值從小到大的順序依次篩選象元作為樣本，第一次選出鄰域分析結果值最小的樣本，如果下一次選出的樣本與已選出樣本的距離在一個鄰域內，則舍棄該樣本繼續選出下一個樣本，當樣本數量達到計算得出的該圖斑單元目標數量時為止。

其中，鄰域分析時采用的反映環境因子在空間分布上變化大小的指標定義如下：

式中，M表示鄰域內環境因子空間變化的程度，N表示鄰域內象元總個數，Di表示第i個鄰域象元因子值，D0表示自身象元因子值，M表示某一位置象元環境因子在當前空間鄰域內的變化程度。鄰域分析結果M越小表示該位置的環境因子數值與周圍的差距越小，空間變化最穩定，因而該位置的土壤類型可能是圖斑所對應的土壤類型。

基于本文的上述假設，在一個圖斑內，土壤類型訓練樣本出現在環境因子鄰域分析結果最小的位置。但對于不同的環境因子，這些位置可能存在較大差異，例如圖2顯示了黃壤性土類型海拔分布較大的某一圖斑內，以3×3柵格窗口為鄰域分別對高程因子和坡度因子進行鄰域分析，空間變化最小5%、10%和20%象元的分布情況，其中空間變化最小5%是指圖斑單元內鄰域分析結果值小于5%分位數的象元，即占總數目5%的空間變化較小的象元位置，10%和20%同理。

圖2 不同因子空間變化最穩定象元位置分布Fig. 2 Distribution of the most stable cells in spatial variation relative to covariate

由圖2可知，高程因子空間變化最小的象元位于山脊、山腳等高程變化小、平緩的位置，空間分布高度集聚，圖斑內絕大部分空間位置未被覆蓋。對坡度因子進行鄰域分析得出坡度的空間變化情況，空間變化最小的象元位于山脊與山腳中部的地形變化穩定的位置，并且覆蓋到圖斑大部分范圍。所以，需要進行不同位置圖斑、不同環境因子的討論用以說明本文所提出方法的不同效果。

由于地形要素基本上可以代表土壤形成與發展過程中的主要影響要素［6］，本文對高程和坡度兩個主要的地形要素環境因子進行討論。本文基于高程和坡度兩個重要的環境因子，選取3×3和7×7兩種柵格窗口進行鄰域分析篩選樣本，由于研究區土壤圖各土壤類型圖斑最小面積的中數為36 250 m2，而環境因子空間分辨率為30 m，所以不宜采用較7×7更大的柵格窗口作為鄰域。

1.5 結果評價

評價樣本篩選結果的指標之一就是樣本的空間分布特征，不同數量的訓練樣本空間分布應覆蓋大部分圖斑以增加樣本對圖斑全局的代表性，即盡可能完備地篩選出包含一系列接近圖斑土壤類型所對應的成土環境特征。此外，所采集的訓練樣本數值應該各不相同、具有較大的信息量，以保證樣本對全局具有較高的代表性。

本文將各土壤圖圖斑所選出的訓練樣本進行位置分布制圖以分析訓練樣本的空間分布，計算差異比例和標準差兩個指標來評價樣本的信息量。差異比例是指一組樣本中除去重復值后的樣本個數與該組樣本總個數之間的比例，差異比例越大，說明樣本重復程度低，差異程度高。在數理統計中，方差可作為評價信息量的一個簡易指標，如在主成分分析中就是尋找方差最大的方向以保證較小的信息損失，方差越大表示樣本的信息量越大，因而本文采用與方差相關的標準差這一統計指標。

本文還采用已有研究中環境因子直方圖峰值區間的方法基于高程和坡度確定樣本，通過對比來評價說明本文提出的鄰域分析方法確定樣本方法。環境因子直方圖峰值區間的方法能夠篩選出訓練樣本，被廣泛用于從歷史土壤圖篩選訓練樣本，并進行歷史土壤圖知識挖掘［6-7］和歷史土壤圖的更新［3-4,9］，具有較高的更新精度。

2 結果與討論

2.1 各土壤類型目標樣本數量

通過上述面積分段線性縮放的方法，計算土壤圖各個圖斑單元的樣本數量，最后對圖斑所屬土壤類型進行匯總，得出各土壤類型的樣本數量（表2）。

表2 各土壤類型目標樣本數量Table 2 Quantity of target samples relative to soil type

其中原始面積和原始圖斑數是在面積閾值篩選之前的土壤類型面積和土壤類型圖斑單元數，采樣面積和采樣圖斑數是面積閾值篩選之后的土壤類型面積和土壤類型圖斑單元數，目標樣本數量即為所屬土壤類型所有采樣圖斑單元樣本數量的總和。面積閾值篩選用較少的圖斑單元數保留了較大的面積，而目標樣本數量的大小大致與土壤類型采樣面積和采樣圖斑數相匹配。對不同土壤類型的具體圖斑單元而言，則具有與其面積對應的樣本數量（表3），表中圖斑為后文樣本空間位置展示時所選的研究區典型土壤類型圖斑單元。

表3 典型土壤類型圖斑樣本數量Table 3 Quantity of samples of a typical soil map polygon

2.2 典型土壤圖斑單元樣本空間位置

對于樣本空間位置分布結果，本文選取研究區三種典型土壤類型的典型圖斑進行展示分析。由于研究區山體、丘陵分布較廣，所以選取研究區面積較大或者圖斑較多、出現在山體不同位置的土壤類型作為典型土壤類型，即選取隨海拔升高發育的土壤類型為序列，分別為潴育型水稻土、紅壤性土和黃壤性土。潴育型水稻土分布在河流兩岸、低平丘陵的底部區域；紅壤性土主要分布在600～700 m以下的低山地區；黃壤性土分布在600～700 m以上的山坡。三種典型圖斑內，基于高程和坡度的不同環境因子，采用不同方法確定樣本的空間位置（圖3）。

圖3 典型土壤圖斑單元訓練樣本空間分布Fig. 3 Spatial distribution of training samples in typical soil map polygon

相比于直方圖峰值區間方法，鄰域分析方法是根據環境因子空間變化，篩選出一系列成土環境穩定、對圖斑具有全局代表性的空間位置。直方圖峰值區間的方法只在直方圖峰值一個數值區間內篩選樣本，樣本相似度高，篩選樣本不全面，并且可能會出樣本的局部集聚，造成樣本冗余。采用鄰域分析方法確定樣本位置，當圖斑位于地勢平緩的區域時，如潴育型水稻土，基于高程因子和基于坡度因子確定的樣本位置空間分布差異較小，而當圖斑位置位于山區時，如黃壤性土，基于高程因子采用鄰域分析方法確定的樣本位置會出現集聚，而基于坡度因子采用鄰域分析方法確定的樣本處于地形變化穩定的位置，樣本分布在圖斑范圍內覆蓋更大，樣本對圖斑分布的代表性更高。因而，當圖斑位于山區時，基于坡度因子比基于高程因子鄰域分析篩選的樣本，對圖斑整體分布范圍具有更高的代表性。

2.3 典型土壤圖斑單元樣本信息量

同樣以3種典型土壤為例，每種土壤選取5個圖斑單元，分別計算不同方法基于高程和基于坡度確定各圖斑單元樣本高程信息和坡度信息的差異比例以及標準差，其中差異比例圖所示（圖4）。

各圖斑單元樣本標準差如表4所示：

圖4 不同土壤圖圖斑樣本差異比例Fig. 4 Odds ratio of difference in samples relative to soil map polygons

表4 不同土壤圖圖斑樣本標準差Table 4 Standard deviation of samples relative to soil map polygon

續表

無論是基于高程因子還是坡度因子，鄰域分析方法確定的樣本其差異比例和樣本標準差均大于環境因子直方圖確定的樣本，這是因為直方圖峰值區間范圍內的樣本的環境因子值比較接近，樣本差異比例較小，標準差也較小。鄰域分析方法是從周圍地理環境的角度確定樣本位置，未從某一數值區間內篩選樣本，因而確定的樣本相比于直方圖峰值區間樣本具有較大的差異比例和標準差，具有較高的信息量，全局代表性更高，其中鄰域大小對差異比例和標準差的影響隨圖斑而異。

但受歷史土壤圖空間分布精度的影響，從歷史土壤圖上篩選出的樣本精度各異，因而歷史土壤圖精度是影響土壤圖訓練樣本的主要因素。

3 結 論

對于歷史土壤圖圖斑單元采用面積分段線性縮放的方法確定樣本數量，可保證面積和圖斑數與訓練樣本數量之間的對應，消除了過大面積土壤類型確定樣本數量時造成的樣本數量差異懸殊，圖斑面積與樣本數量直接匹配，利于從圖斑直接篩選訓練樣本。采用鄰域分析方法，當圖斑位于地勢平緩的區域時，基于高程因子和坡度因子確定的訓練樣本空間分布差異較小；當圖斑位于山區時，基于坡度因子所確定的訓練樣本處于地形變化穩定的位置，對圖斑的空間分布也具有更高的代表性，較高程因子更適合訓練樣本位置的篩選。此外，鄰域分析方法確定的訓練樣本較環境因子直方圖方法確定的樣本具有更高的差異比例和標準差，樣本信息更豐富，樣本對全局的代表性更高。在土壤圖不同土壤類型的圖斑內，基于面積分段線性縮放法確定合適數量、基于環境因子鄰域分析法確定特定位置的訓練樣本后，可以進行后續土壤圖知識挖掘和土壤圖更新的研究，其中不同土壤類型篩選樣本時的環境因子的選取以及不同位置、大小圖斑所需要的空間鄰域大小的確定，均可進一步根據相關土壤地理學知識進行優化。

［ 1 ］ 席承藩，章士炎. 全國土壤普查科研項目成果簡介. 土壤學報，1994,31（3）：330—335 Xi C F，Zhang S Y. Brief introduction on achievements in national soil survey project since 1979 （In Chinese）. Acta Pedologica Sinica，1994，31（3）：330—335

［ 2 ］ 楊琳，Fahmy Sherif，Jiao You，等. 基于土壤—環境關系的更新傳統土壤圖研究. 土壤學報，2010，47（6）：1039—1049 Yang L，Fahmy S，Jiao Y，et al. Updating conventional soil maps using knowledge on soil environment relationships extracted from the maps （In Chinese）. Acta Pedologica Sinica，2010，47（6）：1039—1049

［ 3 ］ Yang L，Jiao Y，Fahmy S,et al. Updating conventional soil maps through digital soil mapping.Soil Science Society of America Journal，2011，75（3）：1044—1053

［ 4 ］ 黃巍，羅云，汪善勤，等. 基于傳統土壤圖的土壤—環境關系獲取及推理制圖研究，土壤學報，2016,53（1）：72—80 Huang W，Luo Y，Wang S Q，et al. Knowledge of soil-landscape model obtain from a soil map and mapping（In Chinese）. Acta Pedologica Sinica，2016，53（1）：72—80

［ 5 ］ Rad M R P，Toomanian N，Khormali F，et al.Updating soil survey maps using random forest and conditioned Latin hypercube sampling in the loess derived soils of northern Iran. Geoderma，2014，232/234：97—106

［ 6 ］ 朱阿興，李寶林，裴韜，等. 精細數字土壤普查模型與方法. 北京：科學出版社，2008 Zhu A X，Li B L，Pei T，et al. Fine digital soil census model and method（In Chinese）. Beijing：Science Press，2008

［ 7 ］ Qi F，Zhu A X. Knowledge discovery from soil maps using inductive learning. International Journal of Geographical Information Science，2003，17（8）：771—795

［ 8 ］ 楊琳，朱阿興，李寶林，等. 應用模糊c-均值聚類獲取土壤制圖所需土壤—環境關系知識的方法研究. 土壤學報，2007，44（5）：784—791 Yang L，Zhu A X，Li B L，et al. Extraction of knowledge about soil-environment relationship for soil mapping using fuzzy c-means （FCM） clustering （In Chinese）. Acta Pedologica Sinica，2007，44（5）：784—791

［ 9 ］ 劉雪琦，朱阿興，楊琳，等. 土壤圖更新中基于土壤類型面積分級的訓練樣點選擇方法. 土壤學報，2017，54（1）：36—47 Liu X Q，Zhu A X，Yang L，et al. Training sample selection method based on grading of soil types by area for updating conventional soil maps （In Chinese）.Acta Pedologica Sinica，2017，54（1）：36—47

［10］ Odgers N P，Sun W，Mcbratney A B，et al.Disaggregating and harmonising soil map units through resampled classification trees. Geoderma，2014，214/215：91—100

［11］ Odgers N P，McBratney A B，Minasny B. Digital soil property mapping and uncertainty estimation using soil class probability rasters. Geoderma，2015，237/238：190—198

［12］ Yang L，Zhu A X，Qi F，et al. An integrative hierarchical stepwise sampling strategy and its application in digital soil mapping. International Journal of Geographical Information Science，2013，27（1）：1—23

［13］ Silva S H G，Menezes M D D，Owens P R，et al.Retrieving pedologist’s mental model from existing soil map and comparing data mining tools for refining a larger area map under similar environmental conditions in Southeastern Brazil. Geoderma，2016，267：65—77

［14］ Malone B P，Minasny B，Odgers N P，et al. Using model averaging to combine soil property rasters from legacy soil maps and from point data. Geoderma，2014，232/234：34—44

［15］ Pahlavan-Rad M R，Khormali F，Toomanian N，et al. Legacy soil maps as a covariate in digital soil mapping：A case study from Northern Iran.Geoderma，2016，279：141—148