多尺度格網的耕地信息無損提取與表達
——以山西省芮城縣為例

王鵬，畢如田，徐振，高陽

山西農業大學資源環境學院，山西 太谷 030801

資源環境信息通過格網化方法可表達其空間分布與分異規律，可匹配、融合多源數據，實現資源環境數據空間模型的構建和表達（王卷樂等，2015）。格網化方法已在土地利用（Koga et al.，2016）、城市管理（李德仁等，2007）、地理國情監測（李德仁等，2014；王妍程等，2016）、林業（吳東亮等，2007）和生態環境管理（左偉等，2003）等領域得到廣泛應用。尺度是生態學研究的基本概念之一，也被廣泛應用于地理學研究中（岳文澤等，2005；劉學軍等，2007；王曉明等，2005）。在地理學中，尺度指地理范圍、詳細程度及時間頻率（張遷遷等，2015）。多尺度格網是指用不同尺度格網將指定范圍劃分成不同大小層次的網格，實現對全區域的無縫覆蓋，每個層次的網格，在范圍上具有上下層涵蓋關系（李德仁等，2005）。由于氣候、地貌、土壤、土地利用等不同類型的資源環境信息在不同空間尺度上存在著多尺度效應（張遷遷等，2015；畢如田等，2012；高艷等，2010），傳統的格網化方法會導致屬性信息的缺失（白燕等，2010）。Galton（2003）系統地研究了在矢量數據格網化過程中空間屬性隨尺度變化的敏感程度，某些屬性的丟失，實質上是一種面積誤差。連世忠等（2017）提出了一種顧及空間自相關的地理國情信息統計格網尺度選擇方法。因此探索不同尺度下空間信息的無損提取和表達對于準確掌握資源環境信息具有重要作用。

耕地資源是資源環境最重要的組成部分之一，耕地信息的無損提取和格網化表達對耕地資源管理具有重要意義。對耕地的格網化表達研究發現，傳統的格網化方法多以面積占優法（Rule of Max Area，RMA）為基礎，重點研究格網尺度對格網化表達精度的關系，但一定程度上忽略了RMA方法對于格網化表達精度的影響。采用RMA格網化方法的精確性依賴于最佳格網尺度的選擇，格網尺度大于或小于最佳格網尺度時，耕地會呈現擴張或壓縮的趨勢，無法準確表達耕地的數量和質量信息；也無法表達格網單元覆蓋范圍內耕地面積小于設定閾值的耕地面積及其分布情況。對于耕地資源單要素而言，使用RMA方法格網化近似于二值化處理，其結果只能反映某格網單元內是否為耕地，不能表達耕地資源空間分布的漸變規律。鑒于此，在耕地單要素屬性信息無損格網化提取的基礎上，提出了以格網中耕地面積占格網面積的比值即格網耕地純度指數來表達耕地多尺度空間分布的方法，實現了自底向上的多尺度屬性信息快速提取，可根據不同需求從底層數據中提取生成多尺度的耕地格網數據，為耕地資源的精細化管理提供支持。

1 研究方法

1.1 格網劃分方法

地理格網的劃分依據主要有兩種：經緯度坐標格網和高斯-克呂格投影直角坐標格網，其中，前者受緯度影響較大，低緯度和高緯度地區的格網面積和形狀有明顯的差異。后者在同一投影下，所有格網的形狀和大小都相同，利于耕地格網化表達過程中面積的計算（陳彥清等，2014）。由于在耕地格網化表達和數據提取時，要求面積不變或在允許誤差范圍以內，因此本研究選用高斯-克呂格投影直角坐標格網進行分析。

1.2 耕地信息無損格網化表達方法

1.2.1 基礎格網尺度的選擇

在典型地貌土地利用現狀矢量數據柵格化尺度效應的研究中，通過RMA法對該土地利用數據進行不同尺度的柵格化處理，分析得出：當格網單元的大小為30 m×30 m時，在不同典型地貌中耕地面積誤差百分比均小于1%（王健欣，2013），可滿足精度的需求。借鑒全球等經緯度四叉樹剖分體系（Geographic Subdividing grid with One dimension coding on 2n-Tree，GeoSOT）中格網尺度的設置方法，該方法采用2n為各級格網尺度，在向下剖分時可避免小數帶來的計算困難，以實現精確四叉剖分。因此本研究選擇最接近30 m的25m作為基礎格網尺度，既能滿足格網化表達過程中的精度要求，又能為以后更高精度要求小尺度格網劃分提供向下連續剖分的數據基礎。

1.2.2 矢量數據預處理

對矢量數據進行屬性值規范化處理，參照GB/T21010—2007《土地利用現狀分類》，將土地利用類型分為：耕地、園地、林地、草地、城鎮村及工礦用地、交通運輸用地、水域及水利設施用地、其他土地。耕地的土壤類型根據《中國土壤分類與代碼》（GB/T17296—2009）將土類、亞類、土屬、土種進行規范化命名。耕地土壤養分屬性值規范化采用均值替換法（Mean Imputation）處理異常值和空值。

1.2.3 創建基礎格網

本研究在ArcEngine二次開發平臺支持下，基于高斯-克呂格投影平面直角坐標系統，以研究區行政界線外包矩形范圍創建大小為512 m×512 m的格網。采用四叉剖分方法將該格網向下剖分至 32 m×32 m的基礎格網，生成512 m×512 m、256 m×256 m、128 m×128 m、64 m×64 m、32 m×32 m的多尺度格網數據集，在創建格網的同時創建唯一標識碼，用于屬性提取運算。

1.2.4 耕地單要素屬性信息無損格網化表達及格網耕地純度計算

經數據預處理后提取出矢量耕地數據，將其與創建的32 m基礎格網進行疊加分析，計算每個格網單元中的耕地面積。保留格網單元的耕地面積屬性值，保證格網數據中耕地面積與矢量耕地數據的面積完全相同，實現耕地單要素面積屬性的無損格網化表達。根據格網大小和格網中耕地的面積計算格網耕地純度指數I，具體計算公式為：

式中，Ac為格網中耕地的面積；Ag為該格網的面積。通過耕地純度指數定量化表達格網覆蓋范圍中耕地的面積和空間分布。為了形成全區域連續無縫的格網并方便后續提取運算，對與耕地圖斑沒有相交關系的格網采取不刪除的方式，通過RMA方法為格網賦屬性值，為表達耕地空間分布提供參照。具體方法是：在保留耕地面積屬性字段的情況下，新增加 1個地類面積屬性，將該格網與土地利用現狀數據進行再次疊加，通過RMA方法確定格網的地類類型、土壤類型、土壤養分等屬性信息。

1.3 多尺度格網屬性信息提取方法

多尺度格網數據在空間范圍上主要表現為兩種形式：（1）不同尺度格網的尺度縮放因子成整數倍數關系，大尺度格網完整包含小尺度格網，不同層的數據完整覆蓋，不交叉；（2）不同尺度格網的尺度縮放因子不成整數倍數關系，大尺度的格網不能完整包含小尺度格網，不同層的數據存在交叉（張小虎等，2014）。

本研究中多尺度格網數據由頂層格網四叉剖分得到，大尺度的格網完全包含小尺度格網，格網大小成倍數關系，屬于第1種形式。因此在基于基礎格網數據進行屬性提取時，首先確定尺度縮放因子 n（n=2，4，····），n 目標格網的大小為 n×n 個基礎格網。如n=2時，目標格網屬性值由4個基礎格網提取得到，如圖1所示。

圖1 格網信息提取示意圖Fig.1 Sketch map of grid information extraction

相同空間位置上目標格網耕地面積等于 n×n個基礎格網的耕地面積之和，以保證耕地面積屬性的無損提??；目標格網的耕地土壤養分屬性采用地類面積加權平均的方法確定，目標格網的其余各屬性值通過面積占優法提取。面積加權平均計算公式為：

式中，n為尺度縮放因子；Ai為第 i個基礎格網的地類面積；Ni為第i個基礎格網的土壤養分值；N為目標格網的土壤養分屬性值。

根據基礎格網數據的順序編碼和行列數自左下向右上依次進行迭代運算，實現整個研究區數據提取。通過改變尺度縮放因子實現多尺度耕地格網數據提取。多尺度格網生成、多尺度耕地屬性信息提取和耕地純度指數計算均在ArcEngine二次開發平臺下運用C#語言編程實現。

1.4 分析方法

1.4.1 多尺度耕地格網面積誤差分析

分別對多尺度耕地格網數據和縣域矢量耕地數據的耕地面積進行匯總，并計算耕地面積誤差。具體計算公式為：

式中，Ag表示耕地格網化后得出的面積；Av表示耕地矢量數據的面積；E表示耕地面積誤差，正值表示耕地面積比實際大，負值表示耕地面積比實際小。

耕地面積誤差百分比的計算公式為：

式中，Y表示耕地面積誤差百分比，當Y＜0時，表示格網化表達后耕地面積收縮；Y＞0時，表示格網化表達后耕地面積擴張。

本研究通過對使用無損格網化法和使用面積占優法得到的耕地面積與研究區矢量數據的耕地面積進行面積誤差分析。

1.4.2 多尺度耕地格網數據空間分布的一致性檢驗

本研究采用標準差橢圓模型檢驗多尺度耕地格網數據與原始耕地圖斑空間分布情況。標準差橢圓方法（Standard Deviational Ellipse，SDE），是用于精確揭示地理要素的空間分布特征的方法（Warntz et al.，1960）。SDE中的平均中心（Mean Center）是以要素的某一屬性為權重對所有要素的x，y值進行加權平均。SDE的方位角反映分布的主趨勢方向，SDE的長軸與短軸長度為標準距離，分別表示要素在主趨勢和次要方向上的離散程度。首先計算平均中心：

式中，Wi為要素權重；xi，yi為要素的質心坐標；n為要素個數；為要素平均加權中心。其次，計算SDE，公式為：

式中，SDEx為標準差橢圓短半軸長度；SDEy為長半軸長度；n為要素個數。最后計算旋轉角θ，公式為：

本文以格網耕地純度指數和耕地圖斑的面積作為權重，分別計算多尺度耕地格網數據和矢量耕地數據的平均中心和空間分布方向角，比較二者的平均中心和空間分布方向角的差異，差異越小，說明該多尺度格網數據越能正確反映耕地空間分布。

1.4.3 基于耕地純度指數的耕地空間分布尺度特征分析

對各尺度耕地格網數據的耕地純度指數進行分級，統計各級純度指數范圍的格網數及其占整個耕地格網的比重，分析研究區耕地資源在不同尺度下的空間分布特征。對比分析RMA格網化方法與格網耕地純度指數在不同尺度下表達耕地資源空間分布的差異。

2 研究區與數據處理

2.1 研究區及數據來源

本研究以山西省芮城縣為例進行耕地信息無損提取與表達研究。芮城縣位于山西省南端，地處中條山南部，黃河北面，東西長，南北窄，北高南低，為典型黃土丘陵及階地臺塬區，地理坐標為110°16′～110°58′E，34°35′～34°51′N。基礎數據為山西省芮城縣1∶1萬土地利用現狀數據庫。耕地矢量數據從該縣1∶1萬土地利用現狀成果數據中提取得到。數據坐標系為西安80坐標系，高斯-克呂格3度分帶。

2.2 數據處理

（1）對土地利用現狀矢量數據進行屬性值規范處理，以該縣行政區界線的外包矩形為范圍生成512 m×512 m的省級格網，通過四叉剖分方法將該格網向下逐級剖分至32 m×32 m的基礎格網，生成多尺度格網數據集，其規模為 912×2096的共1911552個的格網單元。

（2）將32 m基礎格網與縣域耕地數據進行疊加分析，計算各格網單元中的耕地面積，其中與縣域耕地數據有相交關系的格網有476851個。計算該格網中各格網單元的耕地純度指數，完成耕地單要素無損格網化表達。為實現基于基礎格網的統一的屬性提取運算，保留與縣域土地利用現狀數據沒有相交關系的格網。

（3）將帶有耕地面積屬性的基礎格網和土地利用現狀數據進行再次疊加，得到與縣域土地利用現狀數據有相交關系的格網有1152074個。根據上文研究步驟進行基礎格網屬性賦值。為實現基于基礎格網的統一的屬性提取運算，保留與縣域土地利用現狀數據沒有相交關系的格網。

（4）在ArcEngine二次開發平臺支持下，基于32 m耕地格網數據，運用C#語言依次以縮放因子n=2，4，8，16逐級進行屬性信息提取和耕地純度指數計算，生成64、128、256、512 m多尺度耕地格網數據集。

3 結果與分析

3.1 多尺度耕地格網面積誤差分析

采用上述誤差分析方法進行尺度耕地格網面積誤差分析，結果如表1和圖2所示。通過無損格網化方法得到的多尺度耕地面積誤差均為0.1 m2，而使用RMA方法得到的不同尺度的耕地面積誤差隨格網尺度的變化而變化，格網尺度為32 m和64 m時耕地面積誤差百分比小于1%，格網尺度為128、256和512 m時耕地面積誤差百分比不斷增大，其中格網尺度為 512 m時面積誤差百分比達到了11.65%。

表1 耕地面積誤差統計表Table1 The statistical table of cultivated land area error

圖2 耕地面積誤差百分比Fig.2 The percent of cultivated land area error

3.2 多尺度耕地格網數據空間分布的一致性檢驗

根據式（5）～（11）計算該縣矢量耕地數據和各尺度耕地格網信息數據的平均中心坐標和標準差橢圓旋轉角度，得到該縣矢量耕地分布的平均中心坐標為（37464030.3，3839694.2），耕地空間分布標準差橢圓的旋轉角為80.30°，如圖3所示。多尺度耕地格網數據的平均中心和空間分布標準差橢圓的旋轉角如表2所示，與原矢量耕地數據的平均中心相差 0.1～0.8 m，旋轉角相差 0.01°～0.05°。

圖3 耕地空間分布圖Fig.3 The map about spatial distribution of cultivated land

表2 不同尺度耕地格網標準差橢圓模型計算結果Table2 The calculation results of SDE for cultivated land on different scales

表3 不同尺度格網單元耕地純度分布表Tale 3 The purity distribution of different scale grid cell

3.3 基于格網耕地純度指數的耕地空間分布尺度特征分析

對多尺度耕地格網數據集的格網耕地純度進行分級（0～0.2，0.2～0.4，0.4～0.6，0.6～0.8，0.8～1.0），統計各級格網耕地純度的分布情況，結果如表 3所示。隨著格網尺度從32～512 m逐級增大，純度指數在0.8～1.0的格網數量比重迅速減少，其他純度指數的格網數量比重呈增長趨勢。以研究區內耕地破碎與集中連片交接區域為例，分別制作基于面積占優方法的多尺度耕地空間分布專題圖（圖4）和基于耕地純度指數的多尺度耕地空間分布圖專題圖（圖5）。如圖4所示，在耕地破碎地區耕地的格網隨尺度逐級增大可能會被判定為非耕地，耕地集中連片區域周邊的非耕地格網隨尺度逐級增大可能被判定為耕地，從而導致誤差逐漸變大。如圖 5所示，隨格網尺度逐級增大，耕地純度指數呈現出一定的變化規律：在耕地集中連片區，格網耕地純度指數在原值上下小范圍浮動；在耕地破碎地區，格網耕地純度指數急劇下降，正向值趨于0，表示該區域有耕地分布但面積較小。

圖4 基于面積占優法的耕地空間分布示意圖Fig.4 Thesketch map of spatial distribution of cultivated land based on the Rule of Max Area

圖5 基于格網純度指數的耕地空間分布示意圖Fig.5 The sketch map of spatial distribution of cultivated land based on grid purity index

4 討論

格網化的誤差與格網單元大小、多邊形的結構和算法有很大關系（Shortridge，2004），矢量數據多尺度格網化過程存在精度失損問題（楊存建等，2001）。本研究由耕地面積無損提取方法得到的不同尺度的耕地面積誤差均為0.1 m2，該誤差來源于計算過程中的舍入誤差，相對于全縣488111225.1 m2耕地，可忽略不計。在格網尺度較小時，采用無損格網化方法和RMA方法處理的面積誤差相差不大，隨著格網尺度逐級增大，RMA方法會呈現誤差逐級增大的分布特征（劉明亮等，2001），說明RMA格網化方法的準確性依賴于格網尺度的大小。因此，基于多尺度耕地信息無損提取方法，可以在任何尺度下準確表達耕地的面積信息。

對該多尺度格網數據耕地面積進行提取匯總，能滿足統計上的一致性，在空間分布一致性檢驗中，陳彥清等（2014）通過標準差橢圓模型對不同尺度耕地資源空間分布相似性進行了驗證。標準差橢圓模型可以查看要素分布是否被拉伸，并以中心、長軸、短軸、方位角為標準差橢圓的基本參數（Zhao et al.，2012）。本研究中該縣域面積為1178.8 km2，從表 2可以看出，各尺度耕地格網數據與原矢量耕地分布的平均中心相差均小于1 m，偏差很?。桓鞒叨雀馗窬W數據與原矢量耕地空間分布標準差橢圓的旋轉角相差均小于0.1°，在0°～180°范圍之間，偏差可忽略不計。因此，該多尺度格網數據的耕地空間分布與原矢量耕地數據保持高度一致。

在耕地資源空間分布地圖表達過程中，傳統的面積占優法只能將耕地面積占優的格網表示為耕地，無法表達耕地面積不占優的格網中的耕地分布。耕地純度指數法可將耕地面積不占優的格網中的耕地（如圖5中耕地純度指數在0～0.4范圍內的格網）通過色調進行空間表達，形成了一個連續漸變的耕地純度指數面（圖5），不僅能通過不同色調和尺度特征反映該區域耕地空間分布漸變特征和集中連片情況，而且能通過耕地純度指數和格網大小準確計算格網中的耕地面積。

5 結論

本文以位于黃土丘陵區的山西省芮城縣為實例，對基于多尺度格網的耕地信息無損提取與表達方法進行了驗證，得出如下結論：

（1）通過耕地面積無損提取方法得到的不同尺度的耕地面積誤差均為0.1 m2（可忽略不計）。實現了多尺度格網的耕地面積無損提取，表明無損格網化方法能精確地提取耕地數量信息。

（2）經標準差橢圓模型（SDE）驗證，該多尺度耕地格網數據與原矢量耕地的平均中心相差小于1 m，標準橢圓角旋轉度相差小于0.1°，表明該多尺度耕地格網數據可以準確反映耕地資源的分布情況。

（3）通過基于格網耕地純度指數的耕地空間分布尺度特征分析，得出基于耕地純度指數的方法能夠實現不同尺度下耕地空間分布的定量化表達。對比RMA格網化方法處理結果，基于格網耕地純度指數的多尺度表達方法更能準確表達和展示耕地資源專題信息。

參考文獻：

GALTON A. 2003. Granularity-sensitive spatial attributes [J]. Spatial Cognition and Computation, 3(2-3): 97-118.

KOGA T, KAWAMURA A. 2016. Assessing impervious area ratios of grid-based land use classifications on the example of an urban watershed [J]. Procedia Engineering, 154(9): 609-616.

SHORTRIDGE A M. 2004. Geometricvariability of raster cell class assignment [J]. International Journal of Geographical Information Science, 18(6): 539-558.

WARNTZ W, NEFT D. 1960. Contributions to a statistical methodology for areal distributions [J]. Journal of Regional Science, 2(1): 47-66.

ZHAO Y, YANG Z Y, HAO L S, et al. 2012. The evolution of spatial displacement pattern of China's crude oil flow source-sink system [J].Acta Geographica Sinica, 67(4): 455-466.

白燕, 廖順寶. 2010. 矢量數據屬性信息無損柵格化的實現方法——以全國 1∶25萬土地覆被數據為例[J]. 地球信息科學學報, 12(3):385-391.

畢如田, 高艷. 2012. 典型地貌景觀指數的多尺度效應分析——以山西省運城市為例[J]. 地球信息科學, 14(3): 338-343.

陳彥清, 楊建宇, 張超. 2014. 基于網格的耕地質量分等成果省級匯總方法[J]. 農業工程報, 30(24): 280-285.

高艷, 畢如田, 曹毅. 2010. 空間粒度變化及土地利用分類對景觀指數的影響——以山西省聞喜縣為例[J]. 中國生態農業學報, 18(5):1076-1080.

李德仁, 李宗華, 彭明軍, 等. 2007. 武漢市城市網格化管理與服務系統建設與應用[J]. 測繪通報, (8): 1-5.

李德仁, 邵振峰, 丁霖. 2014. 地理國情信息的多級網格化表達[J]. 地理空間信息, 12(1): 1-9.

李德仁, 邵振峰. 2005. 空間信息多級網格及其功能[J]. 地理空間信息,3(4): 1-3.

連世忠, 丁霖, 陳江平. 2017. 顧及空間自相關的地理國情信息統計格網尺度選擇——以植被覆蓋信息為例[J]. 測繪通報, (3): 46-51.

劉明亮, 唐先明, 劉紀遠, 等. 2001. 基于1 km格網的空間數據尺度效應研究[J]. 遙感學報, 5(3): 183-190.

劉學軍, 盧華興, 仁政, 等. 2007. 論 DEM 地形分析中的尺度問題[J].地理研究, 26(3): 433-442.

王健欣. 2013. 典型地貌土地利用現狀矢量數據柵格化的尺度效應及誤差分析[D]. 山西: 山西農業大學: 58.

王卷樂, 孫九林. 2015. 格網化資源環境綜合科學調查的若干問題思考[J]. 地球信息科學學報,17(7): 758-764.

王曉明, 劉瑜, 張晶. 2005. 地理空間認知綜述[J]. 地理與地理信息科學, 21(6): 1-10.

王妍程, 蔡列飛, 侯繼虎, 等. 2016. 基于GeoSOT模型的地理國情監測多級網格信息統計[J]. 地理空間信息, 14(1): 8-12.

吳東亮, 張旸. 2007. 林業網格化管理系統模型及實現的研究[J]. 北京林業大學學報, 29(S2): 18-21.

楊存建, 張增祥. 2001. 矢量數據在多尺度柵格化中精度損失模型探討[J]. 地理研究, 20(4): 416-421.

岳文澤, 徐建華, 談文琦, 等. 2005. 城市景觀多樣性的空間尺度分析——以上海市外環線以內區域為例[J]. 生態學報, 25(1): 122-128.

張遷遷, 畢如田, 張吳平, 等. 2015. 多尺度 EGLSN的構建及耕地地力模擬研究[J]. 中國生態農業學報, 23(3): 345-353.

張小虎, 鐘耳順, 王少華. 2014. 多尺度空間格網數據的索引編碼研究[J]. 測繪通報, (7): 35-38.

左偉, 張桂蘭, 萬必文, 等. 2003. 中尺度生態評價研究中格網空間尺度的選擇與確定[J]. 測繪學報, 32(3): 267-271.