EN中文

Journal of Geo-information Science >

2015 , Vol. 17 >Issue 7: 798 - 803

DOI: https://doi.org/10.3724/SP.J.1047.2015.00798

Orginal Article

Study on Acquiring Appropriate Scales of Ground Features Based on Monte Carlo Simulation

  • ZHU Junxiang , 1, 2 ,
  • WANG Juanle , 1, 3, *
Expand
  • 1. State Key Laboratory of Resources and Environmental Information System, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China
  • 2. University of Chinese Academy of Sciences, Beijing 100049, China
  • 3. Jiangsu Center for Collaborative Innovation in Geographical Information Resource Development and Application, Nanjing 210023, China
*Corresponding author: WANG Juanle, E-mail:

Received date: 2014-12-31

  Request revised date: 2015-02-22

  Online published: 2015-07-08

Copyright

《地球信息科学学报》编辑部 所有
Fold

Abstract

The semi-variogarm has been widely applied in many fields such as mining, soil science, and environmental science to acquire the impact range of spatial process. In remote sensing, it could be used to analyze the spatial structure of remote sensing images or obtain appropriate scales for ground features in the images. Nevertheless, images with huge sizes make the application of semi-variogram on remote sensing different from other disciplines. The computers, on which the semi-variogram curves are calculated and fitted, need more memory and stronger CPUs, which is seemingly impossible to always meet the requirement. A common solution is to decrease the data volume by random sampling, under this circumstance. Specifically, a small amount of samples are selected randomly from the population and analyzed for supposed result instead of using the whole population. This method does decrease the requirement of analysis for computing capacity and memory of computers. However, the accuracy of analysis drops simultaneously, because the result derived from samples in one single sampling is likely to contain errors. To solve this problem that related to spatial structure analysis, this study proposed a method based on Monte Carlo simulation in which a small amount of samples are taken from the huge-volume population for enormous times and then analyzed respectively. In this way, the amount of computation in a single simulation can be reduced to the level that an average computer could tolerate, meanwhile the accuracy can be guaranteed. Also, the parallel computing technology was introduced in this study in order to minimize the time needed for simulation. The parallel computing of semi-variogram was executed in MATLAB whose parallel computing service is simple and easy to manipulate. The experimental area is a rectangular part of An'sai County of Shaanxi Province, China, with an area of 41.32 square kilometers. In this area, there are many types of ground features such as forest, grass, water, farmland and built area. Among all these features, grass is the dominant one. The simulation result shows that, this method could acquire appropriate scales for the common ground features while keeping the estimation errors low.

Key words： Monte Carlo simulation; semi-variogram; appropriate scale; grid

Cite this article

ZHU Junxiang , WANG Juanle . Study on Acquiring Appropriate Scales of Ground Features Based on Monte Carlo Simulation[J]. Journal of Geo-information Science, 2015 , 17(7) : 798 -803 . DOI: 10.3724/SP.J.1047.2015.00798

1 引言

半方差函数（也称半变异函数）作为克里金插值的前处理,是一种非常重要的分析方法。目前,半方差分析不仅在地质采矿领域应用广泛,而且在遥感调查、生态评价、经济分析、土壤监测等领域已得到运用和发展[1-4]。
在诸多应用领域（如地质、生态、土壤等）中受多方条件限制,样本量一般保持在103以内[5],然而在遥感应用领域,影像中一个像元就是一个数据,影像数据量可达到106量级甚至更大。庞大的数据量对计算机的计算、存储性能和数据处理算法,提出了较其他领域更高的要求。尤其近年来,更多高空间分辨率数据的出现、发展,使得遥感影像数据量呈几何增长态势,致使这一矛盾更加突出。
遥感信息提取适宜尺度的判别是半方差函数在遥感领域的重要应用之一。由于遥感数据量大,在用半方差函数分析遥感影像空间结构时,往往需要对样本量进行限制。传统做法是从影像中直接随机选择若干像元进行分析,如Woodcook在1988年将半方差分析初次应用于影像空间结构分析时,从影像中选择了600个点[6-7]。这种处理方式能降低对计算机性能的要求,但其以一次随机抽样的模拟结果来代表总体,难免会出现以偏概全的情况。
针对该问题,借鉴Monte Carlo模拟思想,本文通过对大影像以小样本大量随机重复采样,分别计算、模拟各次抽样的半方差函数,绘制所有模拟结果的概率分布曲线,并用统计分析的方法来获取地物的适宜尺度。

2 试验区和数据

（1）研究区位于陕西省延安市安塞县沿河湾镇（109°20'16''~109°25'13''E,36°43'37''~36°46'38''N）,面积约41.32 km2。研究区内的地物类型主要有森林、草地、水体、耕地、建筑等。延安市第二大河延河由西北向东南流经该研究区,河流两侧分布着研究区内大部分的耕地和建筑。森林整体分布比较稀疏,局部较为集中,草地为优势地类,占试验区面积很大比例（图1）。
Fig. 1 The test site location

图1 试验区位置图

（2）采用的数据为Rapideye L1B遥感影像及由其解译的土地覆被图。影像空间分辨率为6.4 m,大小为1134×887像元,使用前已经过正射校正、大气校正。Rapideye影像包含蓝、绿、红、红边和红外5个波段,因红色波段对植被较为敏感,故本研究采用其作为分析波段[7-8]。土地覆被图分类总体精度为89.6%,Kappa系数为0.88,用于从原始影像中提取出单一地物的影像,以减少其他地类对半方差分析的干扰。

3 地物适宜尺度分析

3.1 提取原理和方法

本研究所涉及的原理和方法主要包括Monte Carlo模拟、半方差函数和并行计算。
Monte Carlo法最为突出的特点是大量重复随机试验,其应用模式可归为3类：最优化、数值积分和绘制概率分布。虽然其形式多样,但模拟过程基本服从以下范式：（1）界定输入范围;（2）从输入的概率分布中产生随机数;（3）根据随机数进行抽样并计算初步结果;（4）对初步结果进行统计汇总[9-11]。
本研究以遥感影像为输入,从中随机选择若干像元,根据像元值及其空间位置计算、拟合半方差函数,并对半方差函数拟合所得变程进行统计分析。模拟结果的精度与模拟次数有关,次数越多,所得结果与真值的近似程度越高。根据Monte Carlo检验中置信度与模拟次数的关系,本研究进行999次模拟,以期估计变程99%置信区间。
计算遥感影像的半方差函数如式（1）所示[12]。
γ ( h ) = 1 2 N ( h ) | ∑ { Z ( s i ) - Z ( s i + h ) } （1）
式中, N ( h ) 代表相互间距离为h的点对的数量; Z s i 和 Z s i + h 分别为点 s i 和点 s i + h 处的值。
图2是一个典型的半方差曲线[13],主要参数包括基台值(s)和变程(h)。基台值代表变量的最大变异程度,是对真实方差的一个估计[14]。变程代表变量的影响范围,在分析遥感影像的空间结构时,代表了地物的尺寸大小。
Fig. 2 An ideal semi-variogram

图2 理论半方差图

半方差计算和拟合,需要设定系列参数,主要包括：选用的模型、最大采样距离、拟合点数及样本量。本研究经过多次试验,选用指数模型,最大采样距离设为500个像元,拟合点数设为50,样本量选择10 000。除此之外,假设各地类均为各向同性分布,即其空间分布仅与距离有关,与方向无关。
Monte Carlo模拟需进行大量重复运算,使用多核计算机或集群技术并行计算可提高数据处理效率[15]。MATLAB平台为用户提供了高性能的交互式计算环境,广泛应用于工程计算,其从7.0版本开始,提供并行计算工具箱、支持并行计算[16-17]。
本研究借助MATLAB 2012b,在台式多核计算机上,实现Monte Carlo模拟的并行计算。MATLAB在单一计算机上最高支持12个工作核,可同时进行12项模拟[18]。本研究所用台式机主要配置：主频为3.1 GHz的i5-2400四核处理器,14 G内存。

3.2 模拟结果和分析

3.2.1 实验流程
以解译的分类结果数据作为输入,分别对森林、草地、建筑、农田、水体等主要类型,进行半方差分析的Monte Carlo模拟（图3）。模拟结果分析主要包括：异常值检验、概率分布特征、正态性检验,以及适宜尺度的区间估计。
Fig. 3 Experiment workflow

图3 实验流程

（1）异常值检验。将位于 [ Q 1 - 1.5 × ( Q 3 - Q 1 ) , Q 3 + 1.5 × ( Q 3 - Q 1 ) ] 范围之外的值列为异常值,通过重复检测直至消除所有异常值。其中,Q1是下四分位点,Q3是上四分位点。
（2）概率分布特征。概率分布特征主要由概率分布图和箱线图来表征。箱线图由最小值、下四分位、中值、上四分位点、最大值,以及2个极值和四分位点之间的虚线构成,可反映数据的对称性和集中性等整体分布情况[19]。
（3）正态性检验。正态性检验用QQ图（Quality-Quality Plot）完成。QQ图是一种通过画出分位数来比较,2个概率分布的方法,也可用来检验一组数据的经验分布和理论分布是否一致。如果被比较的2个分布比较相似,则其QQ图近似地位于 y = x 直线上[20]。
3.2.2 试验结果和分析
模拟耗时约36 h,获得建筑、耕地、森林、草地和水体5种类型共4995个拟合变程值。
对变程进行异常值检测,从建筑、耕地、森林、草地和水体的拟合结果中分别检查出14、12、36、19和13个异常值,将这些异常值设为空值（NAN）。自定义一种结合了概率分布和箱线图的三维图（图4）。在该三维图中,X轴为地物类型,Y轴为变程,Z轴为变程的频率。X-Y平面内绘制有箱线图,X-Y-Z空间内绘有5条表征各类型变程分布特征的概率分布曲线。从概率分布图可较为直观地看出模拟结果分布的正态性特征,从箱线图的中值（短红线）可看出各地物类型适宜尺度之间的差异。
Fig. 4 The probability distribution of simulation ranges

图4 模拟变程概率分布

为更直观地观察各类型变程的概率分布特性,分别绘制每种类型的QQ图（图5）。图5中,纵坐标为排好序的模拟结果（经验分位点）,横坐标为这些数据正态分布的理论分位点,蓝色十字为实际对比点,红色虚线为理论参考线。由图5可看出,各条曲线除两端稍有偏离外,基本都落在参考线上,由此可判定各类型变程正态分布。
Fig. 5 The normality test of all simulation results

图5 各地物类型模拟结果的正态性检验

表1为拟合变程的统计分析结果,包含中值、均值和方差等。由表1看出,中值与均值相近,因而可认为各组数据基本对称分布,2个统计量都可用来表征适宜尺度。
Tab. 1 The statistical analysis of simulation results

表1 模拟结果统计分析

建筑(m) 农田(m) 森林(m) 草地(m) 水体(m)
中值 159.3 334.9 307.1 294.1 252.0
均值 160.6 335.7 308.4 294.4 253.9
方差 277.9 366.6 695.0 437.8 1664.3
以均值为适宜尺度,则建筑、农田、森林、草地和水体的适宜尺度分别为160.6、335.7、308.4、294.4和253.9 m。此处所指的建筑、农田、森林等并不是单栋建筑、单片农田或单棵树木,而是建筑、农田和树木的集合体。
若以方差作为不确定性指标可看出,建筑、农田、森林和草地的不确定性基本在同一水平,而水体的方差远大于其他几组数据,为最小方差的6倍左右,表明其不确定性较高,其可能的原因是水体（主要是河流）特征与本研究假设不符。本研究假设地物各向同性,而河流具有明确的方向性（从高到低）,即各向异性。因此,本研究的方法不太适合道路、河流等细长地物表达适宜尺度的提取分析。
置信区间长度与采样次数的开方成反比,模拟次数越多,区间长度越小,估计越精确。从表2可看出,建筑、农田、森林、草地和水体的适宜尺度的99%置信区间分别为[159.2,161.9]、[334.1,337.2]、[306.2,310.5]、[292.7,296.1]和[250.6,257.3],无论是95%还是99%置信区间,其区间长度都较小。
Tab. 2 The 95% and 99% confidence intervals of all ground features

表2 各地物类型95%、99%置信区间

建筑（m） 农田（m） 森林（m） 草地（m） 水体（m）
95%置信上限 161.6 336.8 310.0 295.7 256.5
95%置信下限 159.5 334.5 306.8 293.1 251.4
99%置信上限 161.9 337.2 310.5 296.1 257.3
99%置信下限 159.2 334.1 306.2 292.7 250.6
本文得出的适宜尺度是Monte Carlo多次模拟结果的均值,代表了同种类型大多数对象的尺度特征,可用于遥感影像选择。从该结果可见,农田的适宜尺度相对最大,这与安塞以农业生产为主的特点对应;森林和草地的适宜尺度较大,这与本区域作为国家退耕还林还草的重点区域,林地和草地的斑块较为规则、连片有关;水体和建筑的适宜尺度略小,这与该地区地处黄土高原腹地,其水系发育受地形影响大,且城镇分布与水体分布区域相近有关,因而二者的适宜尺度都较小。
Woodcock认为,小于影像像元的物体不能被检测到,小于2-3倍影像分辨率的物体不能被有效地检测[6]。以试验区内建筑为例,分析结果意味着分辨率大于160.6 m的遥感影像无法表现建筑特征,分辨率大于53.5 m（适宜尺度三分之一）而小于160.6 m的影像不能有效地体现建筑特征,而分辨率小于53.5 m的影像可较有效地体现建筑特征。但要注意,半方差函数分析所得适宜尺度与场景相关,具有地域性,某一区域的结果不能直接应用于其他区域,但尺度的模拟和研究方法具有普遍性,能应用于不同区域。
半方差分析结果与采用的影像波段相关。遥感影像各波段对地物的敏感度不同,因此,不同波段下同一场景影像的空间结构存在差异,用相同方法分析其他波段时,得出的结果不同。最后,减小Monte Carlo模拟误差在于选择合理的模拟参数——最大采样距离、拟合点数、样本量及模型。这些参数需要控制变量多次试验来确定。4个参数中,最大采样距离对抽样影响较大,很大程度上决定了模拟精度。过大的最大采样距离会造成较多的跨区域点对（2个像元同属一种类型,但位于相离区域内）,尤其是同类地物区域较为密集时。因同属一类,这些像元的DN（Digital Number）值较为接近,会造成距离增大、相关性增强的假象,在半方差图上表现为曲线中后段呈下降趋势。前期试验中确定的最大采样距离应大于变程,同时半方差曲线后半段应较为平稳。总之,分析过程中,像元的选取是随机的,像元分布于各个相离的同类地物区域,一个合适的最大采样距离,可减少跨区域点对,使得大部分点对在同一区域内,这对减少模拟误差较为关键。

4 结论

针对遥感信息提取地物适宜尺度研究中的抽样样本量大、运行效率低的难题,本文提出一种以Monte Carlo模拟来估算影像空间结构的方法,结合并行计算,能够高效地通过半方差函数运算和分析获取特定遥感影像中的地物适宜尺度。在我国中部的黄土高原安塞县沿河湾镇的实验表明,其能高效获得建筑、农田、森林、草地和水体的适宜尺度区间,具有较高的精度一致性。通过本研究,预期为遥感地学分析中的大数据处理量与有限计算能力之间的矛盾、地物适宜尺度传统大样本数据分析时,以偏概全等问题的解决有所启示。
致谢：感谢国家地球系统科学数据共享平台提供遥感数据。

The authors have declared that no competing interests exist.

References

Publishing order | Descend order by publishing year | Descend order by cited within
[1]
郑袁明,陈煌,陈同斌,等.北京市土壤中Cr, Ni含量的空间结构与分布特征[J].第四纪研究,2003,23(4):436-445.

[2]
岳文泽,徐建华,谈文琦,等.城市景观多样性的空间尺度分析——以上海市外环线以内区域为例[J].生态学报,2005,25(1):122-128.

[3]
沈国状,廖静娟.基于半变异函数的多极化SAR图像地表淹没程度分析[J].遥感技术与应用,2006,20(6):569-573.

[4]
梅志雄. 基于半变异函数的住宅价格空间异质性分析——以东莞市为例[J].华南师范大学学报:自然科学版,2008(4):123-128.

[5]
WEBSTER R, OLIVER M A.How large a sample is needed to estimate the regional variogram adequately?[M]. Geostatistics Tróia'92, Springer, 1993:155-166.

[6]
WOODCOCK C E, STRAHLER A H, JUPP D L.The use of variograms in remote sensing: I. Scene models and simulated images[J]. Remote Sensing of Environment, 1988,25(3):323-348.

[7]
WOODCOCK C E, STRAHLER A H, JUPP D L B. The use of variograms in remote sensing: II. Real digital images[J]. Remote Sensing of Environment, 1988,25(3):349-379.

[8]
ATKINSON P M, CURRAN P J.Choosing an appropriate spatial resolution for remote sensing investigations[J]. Photogrammetric engineering and remote sensing, 1997,63(12):1345-1351.

[9]
ROBERT C P, CASELLA G.Monte Carlo statistical methods[M]. Springer, 1999.

[10]
曲双石,王会娟.Monte Carlo方法及其应用[J].统计教育,2009(1):45-55.

[11]
METROPOLIS N.The beginning of the Monte Carlo Method[J]. Los Alamos Science, 1987(15):125.

[12]
JUPP D L, STRAHLER A H, WOODCOCK C E.Autocorrelation and regularization in digital images. I. Basic theory[J]. Geoscience and Remote Sensing, IEEE Transactions on, 1988,26(4):463-473.

[13]
CLARK I.Practical Geostatistics[M]. Applied Science Publishers London, 1979.

[14]
LEVESQUE J, KING D J.Airborne digital camera image semivariance for evaluation of forest structural damage at an acid mine site[J]. Remote Sensing of Environment, 1999,68(2):112-124.

[15]
陈国良,孙广中,徐云,等.并行计算的一体化研究现状与发展趋势[J].科学通报,2009(8):1043-1049.

[16]
SHARMA G, MARTIN J.MATLAB®: a language for parallel computing[J]. International Journal of Parallel Programming, 2009,37(1):3-36.

[17]
CHOY R, EDELMAN A.Parallel MATLAB: Doing it right[C]. Proceedings of the IEEE, 2005,93(2):331-341.

[18]
刘维. 实战MATLAB之并行程序设计[M].北京:北京航空航天大学出版社,2012.

[19]
庄作钦. BOXPLOT——描述统计的一个简便工具[J]. 统计教育,2003(1):16.

[20]
宗序平,姚玉兰.利用QQ图与PP图快速检验数据的统计分布[J].统计与决策,2010(20):151-152.

Options
Outlines

/

〈 〉