CN1651860A - 用不同尺度遥感数据估计面积变化的对称系统抽样技术 - Google Patents

Abstract

一种以不同尺度遥感数据为数据源，采用对称系统抽样技术，估计农作物种植面积、森林覆盖面积、城市用地面积等地物面积要素变化的方法，主要应用于遥感工程、农林业工程、城市规划、国土资源调查和统计领域。利用中低分辨率遥感影像，如MODIS、TM，构造抽样基本单元；用面积变化大小的初略值为辅助变量进行排序；采用对称系统抽样技术选取代表性强的样本点；以此样本信息推算研究对象在整个区域的面积变化情况。可以广泛应用于各种地物面积要素变化的监测。

Description

用不同尺度遥感数据估计面积变化的对称系统抽样技术

技术领域

本发明涉及一种以不同尺度遥感数据为数据源，采用对称系统抽样技术，估计农作物种植面积、森林覆盖面积、城市用地面积等地物面积要素变化的方法，主要应用于遥感工程、农林业工程、城市规划、国土资源调查和统计领域。

背景技术

不同尺度遥感数据，是指具有不同空间分辨率的遥感影像。对称系统抽样技术，是指对辅助变量进行排序，用对称的方法获取样本点，估计总体变量。我国是一个农林产品大国，同时正处于城市化迅速发展的时期，农作物种植面积、森林覆盖面积、城市用地面积等都是国家社会经济的基础信息，关系国计民生，对于制定国家和区域社会经济发展规划，制定农林产品进出口计划，确保国家粮食与生态安全，加强生态环境保护，掌控农林产品期货市场信息与城市用地信息，提高国家宏观调控能力和政府决策能力具有重要价值。

国家统计局系统对于农作物种植面积和森林覆盖面积通常采用按国家、省、地区、县四级行政单元逐级上报的方法，或者利用直属农业调查队的调查数据进行统计抽样。由于上报数据的质量无法控制，而多阶段抽样方法会使误差增大，造成最终结果精度不高，整个流程耗时长。对于城市用地面积而言，国土资源部每年采取顺序挑选部分城市进行遥感城市用地面积变化监测，完成一次全覆盖需要几年时间，很难及时全面掌握全国城市用地面积的变化情况。

国外普遍采用遥感数据与抽样技术结合的方法，如美国的大面积农作物估产计划(LACIE)、农业和资源的空间遥感调查计划(AGRISTARS)等采用了面积抽样框(Area Sampling Frame)方法，上述抽样方法及框架布设分两种情况：一种是有国家级详细历史数据的国家，另一种是低于国家级数据(省、州)的地区。前者采用二阶随机分层抽样布点，后者采用标准分层抽样布点。欧盟农业遥感研究项目(MARS)通过抽样设置作物面积遥感监测样区(60个)，采用分层抽样方法，实现了欧盟17种作物的面积遥感监测。国内农业遥感抽样研究始于80年代初，先后开展了“京津冀冬小麦遥感估产研究”、“北方11省市冬小麦气象卫星遥感综合估产技术研究”、“全国主要农作物遥感调查”等。但是，上述研究采用现有国家统计局的辅助统计数据和遥感高分辨率影像数据，未能充分利用低分辨率遥感影像资源，而且多采用随机抽样或分层抽样方法，因此，无论在遥感影像资源运用上，还是在抽样技术运用的科学性和精确性上都有很大的局限性。对于城市用地的面积变化估计尚未建立一套系统的方法。

本项发明首次将对称系统抽样技术应用到遥感领域，用于解决遥感应用于农作物种植面积、森林覆盖面积、城市用地面积等面积要素的变化估计上，在国内外同类技术中属于首创。将不同尺度遥感数据相结合，采用对称系统抽样技术，建立一套基于网格的、适用于农林用地和城市用地面积变化估算的方法体系，抽样统计估计量的精度和估算方法的时效性都得到了不同程度的提高，具有巨大的社会效益和经济效益。

主要引证文件：1、水稻遥感估产，王人潮等著，中国农业出版社，2002；2、中国水稻遥感动态监测与估产，赵锐等编，中国科学技术出版社，1996；3、中国农作物遥感动态监测与估产总论，孙九林等编，中国科学技术出版社，1996；4、中国小麦遥感动态监测与估产，王乃斌编，中国科学技术出版社，1996；5、中国玉米遥感动态监测与估产，万恩璞等编，中国科学技术出版社，1996；6、中国农作物遥感动态监测与估产集成系统，熊利亚编，中国科学技术出版社，1996；7、农业信息科学与农业信息技术，王人潮等著，中国农业出版社，2003；8、基于3S技术的土地利用/土地覆盖变化野外采样框架设计一以东北地区黑龙江省为例，徐新良等，遥感技术与应用，2002；9、遥感和空间分析支持下的耕地、城镇动态变化空间采样方法研究，《空间信息多重采样设计的空间统计学应用研究》课题组，统计研究，1999；10、土地资源遥感调查中的空间信息多重采样框架设计与GIS实现，《空间信息多重采样设计的空间统计学应用研究》课题组，统计研究，1999年第1期。

发明内容：

本项发明的目的，是以遥感数据源代替传统地面实地丈量，以对称系统抽样技术代替传统的随机抽样或分层随机抽样方法，从而获得代表性强、估计量效率高的采样点，用于估计农作物种植面积、森林覆盖面积、城市用地面积等面积要素的变化。本项发明是一项遥感工程、农林业工程、城市规划、国土资源调查和统计交叉领域的创新技术，也是一项具有宽广开发前景的应用技术，对我国农作物面积、森林覆盖面积和城市用地面积的动态监测及其技术手段的提高、更新和换代具有指导价值。

本项发明的主要内容是，利用中低分辨率遥感影像，如MODIS、TM，构造抽样基本单元；用面积变化大小的初略值为辅助变量进行排序；采用对称系统抽样技术选取代表性强的样本点；以此样本信息推算研究对象在整个区域的面积变化情况。

本项发明与现有技术相比，具有如下优点和积极效果：首次将对称系统抽样技术用于遥感面积抽样，极大提高了估计量的精度；将不同尺度遥感数据相结合，建立一套适用于农林用地和城市用地面积变化估算的多用途方法体系，适用面广；采用网格抽样技术，可以和元数据标准相结合，有利于数据共享；使用高时间分辨率的遥感影像源，与统计部门的估计方法相比，时效性大大提高；用遥感数据构造的抽样框，样本代表性强，样本量较传统估计方法大大减少，降低成本。

具体实施方式

本项发明是一项统计学、遥感工程、农林业工程、城市规划和国土资源调查交叉领域的创新技术，实现本项发明的基本技术方案及最佳方式，可以归纳为如下步骤：

1、抽样框的确定

以每个网格内研究对象的变化数为依据进行排序，选择恰当的变量作为排序标志，具体做法如下：

(1)根据研究区域和研究对象的性质，选取1∶10万(或者1∶5万，1∶2.5万等)地形图国际标准分幅作为抽样的基本单元；

(2)收集研究区域过去原有的遥感影像或其他能反映研究对象原来分布的信息(如土地利用图等)；

(3)获得研究区域目前的低分辨率遥感影像；

(4)对(2)、(3)获得的遥感影像进行几何校正、匹配等预处理工作；

(5)收集影响研究区域、研究对象有关的先验信息，包括自然条件、人文条件、社会经济条件以及以往的历史信息；

(6)利用土地利用与土地覆盖分类识别技术，以及先验信息获得研究对象原分布图和现分布图；

(7)将这两幅图叠加，通过土地利用与土地覆盖变化检测技术获得研究对象面积变化的初略值；

(8)按每个网格内面积变化的初略值进行排序，并编制抽样框。

2、计算样本量

因为面积变化初略值和面积变化精确值间存在很强的相关性，可以用面积变化初略值的均值和方差来估算所要抽取的样本量。

(1)计算面积变化初略值X的均值和方差：

$\stackrel{\‾}{X}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{X}_i$

$S^2=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(X_i-\stackrel{\‾}{X})}^2$

(2)计算n₀

$n_0=\frac{t^2{S}^2}{d^2}$ 或 $n_0=\frac{t^2{S}^2}{r^2{\stackrel{\‾}{X}}^2}$

其中：t＝1.96，d＝5％，r＝10％

如果n₀≥N，则取n₀＝N。

(3)如果 $\frac{n_0}{N}<0.05,$ ，则就取n₀，否则对n₀进行修正：

$n=\frac{n_0}{1+\frac{n_0}{N}}$

3、用对称系统抽样法确定抽样单元

以一般的等距抽样为例，假设抽样间距为k，总体单元数为N＝nk，将总体的N个单元排列成k行n列，如下表：

	1     2    … j        … n
		12_r_k	Y1  Yk+1 … Y(j-1)k+1 … Y(n-1)k+1Y2  Yk+2 … Y(j-1)k+2 … Y(n-1)k+2_   _      _          _Yr  Yk+r … Y(j-1)k+r … Y(n-1)k+r_   _      _          _Yk  Y2k … Y(j-1)k+r  … Ynk

对这个线性趋势总体来说，可使用Singn提出的对称系统抽样法提高系统抽样的精度，使系统抽样法有可能达到比分层随机抽样更高的效果。

当N不是n的整数倍，即抽样间距 $k=\frac{N}{n}$ 不是整数时，k可取一个与 最为接近的整数。注意到：当n大于50时，N≠nk造成的干扰很可能是微不足道的。

Singn提出的对称系统抽样具体做法如下：

当N＝nk，n为偶数时，假设起始随机数为r(1≤r≤k)，入样的 对样本单元为：

$[r+\mathrm{jk},N-r-\mathrm{jk}+1],j=\mathrm{0,1,2},\&CenterDot;\&CenterDot;\&CenterDot;,\frac{n}{2}-1$

即y₁、y₂、…、y_n依次为：

$y_r,y_{r+k},\&CenterDot;\&CenterDot;\&CenterDot;{,y}_{r+(\frac{n}{2}-1)k},y_{N-r-(\frac{n}{2}-1)k+1},\&CenterDot;\&CenterDot;\&CenterDot;,y_{N-r-k+1},y_{N-r+1}$

当n为奇数时，仍按以上步骤进行，但到 $j=\frac{n-1}{2}-1$ 后，增加靠近中间的一个单元 $r+\frac{1}{2}(n-1)k.$ 。y₁、y₂、…、y_n依次为：

$y_r,y_{r+k},\&CenterDot;\&CenterDot;\&CenterDot;,y_{r+(\frac{n}{2}-2)k},y_{r+\frac{1}{2}(n-1)k},y_{N-r-(\frac{n}{2}-2)k+1},\&CenterDot;\&CenterDot;\&CenterDot;,y_{N-r-k+1},y_{N-r+1}$

4、获取样本点更高分辨率的遥感影像，读取面积变化的精确值

因为高分辨率的遥感影像受到云量等因素的影响，在抽取样本点时，会存在缺失数据现象。这时可使用与和缺失数据大小最接近的样本点的高分辨率影像资料来代替。如果最接近的样本点的高分辨率影像资料也不好获得，可用初略值来代替。

5、估计量及估计量方差的估计

假设起始值为r，当n为偶数时，相应系统样本的平均值为：

$\hat{Y}=N{\stackrel{\&OverBar;}{y}}_{\mathrm{sy}}=\frac{N}{n}(\underset{j=0}{\overset{\frac{n}{2}-1}{\mathrm{\&Sigma;}}}{y}_{r+\mathrm{jk}}+\underset{j=0}{\overset{\frac{n}{2}-1}{\mathrm{\&Sigma;}}}{y}_{N-r-\mathrm{jk}+1})$

$=\frac{N}{n}(y_r+y_{r+k}+\&CenterDot;\&CenterDot;\&CenterDot;+y_{r+(\frac{n}{2}-1)k}+y_{N-r-(\frac{n}{2}-1)k+1}+\&CenterDot;\&CenterDot;\&CenterDot;+y_{N-r-k+1}+y_{N-r+1})$

当n为奇数时，相应系统样本的平均值为：

$\hat{Y}=N{\stackrel{\&OverBar;}{y}}_{\mathrm{sy}}=\frac{N}{n}(\underset{j=0}{\overset{\frac{n}{2}-2}{\mathrm{\&Sigma;}}}{y}_{r+\mathrm{jk}}+y_{r+\frac{1}{2}(n-1)k}+\underset{j=0}{\overset{\frac{n}{2}-2}{\mathrm{\&Sigma;}}}{y}_{N-r-\mathrm{jk}+1})$

$=\frac{N}{n}(y_r+y_{r+k}+\&CenterDot;\&CenterDot;\&CenterDot;+y_{r+(\frac{n}{2}-2)k}+y_{r+\frac{1}{2}(n-1)k}+y_{N-r-(\frac{n}{2}-2)k+1}+\&CenterDot;\&CenterDot;\&CenterDot;+y_{N-r-k+1}+y_{N-r+1})$

如果把系统样本y₁、y₂、…、y_n看成从各层抽取两个单位的分层随机抽样，可采用以下办法进行方差估计。从第二个样本单元开始，每个样本单元与前一个样本单元组成一对，共n-1对，

的估计为：

$v\left(\hat{Y}\right)=N^{2}v\left({\stackrel{\&OverBar;}{y}}_{\mathrm{sy}}\right)=\frac{N(N-n)}{2n(n-1)}\underset{i=1}{\overset{n-1}{\mathrm{\&Sigma;}}}{(y_{i+1}-y_i)}^2$

6、检察是否需要再抽样本，如果需要，则按前述方法将补充的样本分配下去比较r的值是否达到了10％， $V={\left(\frac{d}{t}\right)}^2={\left(\frac{r\stackrel{\&OverBar;}{Y}}{t}\right)}^2.$ 。如果没有，则需补充样本量，使得总的样本量为：

$n=\frac{s_1^2}{V}(1+\frac{2}{n_1})$

其中： $s_1^2=\frac{1}{2(n-1)}\underset{i=1}{\overset{n-1}{\mathrm{\&Sigma;}}}{(y_{i+1}-y_i)}^2,$ ，n₁是第2步所确定的样本量。

7、推断

总体总值的置信限为：

$\hat{Y}\&PlusMinus;u_{\frac{\mathrm{\&alpha;}}{2}}v\left(\hat{Y}\right)$

Claims (1)

1、一种以不同尺度遥感数据为数据源，采用对称系统抽样技术，估计农作物种植面积、森林覆盖面积、城市用地面积等地物面积要素变化的方法，其特征在于：

(1)抽样框的编制

以每个网格内研究对象的变化数为依据进行排序，选择恰当的变量作为排序标志，具体做法如下：

a、根据研究区域和研究对象的性质，选取1∶10万(或者1∶5万，1∶2.5万等)地形图国际标准分幅作为抽样的基本单元；

b、收集研究区域过去原有的遥感影像或其他能反映研究对象原来分布的信息(如土地利用图等)；

c、获得研究区域目前的低分辨率遥感影像；

d、对(2)、(3)获得的遥感影像进行几何校正、匹配等预处理工作；

e、收集影响研究区域、研究对象有关的先验信息，包括自然条件、人文条件、社会经济条件以及以往的历史信息；

f、利用土地利用与土地覆盖分类识别技术，以及先验信息获得研究对象原分布图和现分布图；

g、将这两幅图叠加，通过土地利用与土地覆盖变化检测技术获得研究对象面积变化的初略值；

h、按每个网格内面积变化的初略值进行排序，并编制抽样框；

(2)计算样本量

因为面积变化初略值和面积变化精确值间存在很强的相关性，可以用面积变化初略值的均值和方差来估算所要抽取的样本量；

a、计算面积变化初略值X的均值和方差：

$\stackrel{\&OverBar;}{X}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{X}_i$

$S^2=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{(X_i-\stackrel{\&OverBar;}{X})}^2$

b、计算n₀

$n_0=\frac{t^2{S}^2}{d^2}$ 或 $n_0=\frac{t^2{S}^2}{r^2{\stackrel{\&OverBar;}{X}}^2}$

其中：t＝1.96，d＝5％，r＝10％

如果n₀≥N，则取n₀＝N；

c、如果 $\frac{n_0}{N}<0.05,$ 则就取n₀，否则对n₀进行修正：

$n=\frac{n_0}{1+\frac{n_0}{N}}$

3、用对称系统抽样法确定抽样单元

以一般的等距抽样为例，假设抽样间距为k，总体单元数为N＝nk，将总体的N个单元排列成k行n列，如下表：

   1     2      …  j        …  n     12_r_k Y₁   Y_k+1  …  Y_(j-1)k+1 …  Y_(n-1)k+1Y₂   Y_k+2  …  Y_(j-1)k+2 …  Y_(n-1)k+2_    _        _          _Y_r   Y_k+r  …  Y_(j-1)k+r …  Y_(n-1)k+r_    _        _          _Y_k   Y_2k   …  Y_jk     …  Y_nk

对这个线性趋势总体来说，可使用Singn提出的对称系统抽样法提高系统抽样的精度，使系统抽样法有可能达到比分层随机抽样更高的效果；

当N不是n的整数倍，即抽样间距 $k=\frac{N}{n}$ 不是整数时，k可取一个与

最为接近的整数；注意到：当n大于50时，N≠nk造成的干扰很可能是微不足道的；

Singn提出的对称系统抽样具体做法如下：

当N＝nk，n为偶数时，假设起始随机数为r(1≤r≤k)，入样的

对样本单元为：

$[r+\mathrm{jk},N-r-\mathrm{jk}+1],j=\mathrm{0,1,2},\&CenterDot;\&CenterDot;\&CenterDot;,\frac{n}{2}-1$

即y₁、y₂、…、y_n依次为：

y_r、y_r+k、…、 $y_{r+(\frac{n}{2}-1)k},y_{N-r-(\frac{n}{2}-1)k+1},$ …、y_N-r-k+1、y_N-r+1

当n为奇数时，仍按以上步骤进行，但到 $j=\frac{n-1}{2}-1$ 后，增加靠近中间的一个单元 y₁、y₂、…、y_n依次为：

y_r、y_r+k、…、 $y_{r+(\frac{n}{2}-2)k},y_{r+\frac{1}{2}(n-1)k},y_{N-r-(\frac{n}{2}-2)k+1},$ …、y_N-r-k+1、y_N-r+1

4、获取样本点更高分辨率的遥感影像，读取面积变化的精确值

因为高分辨率的遥感影像受到云量等因素的影响，在抽取样本点时，会存在缺失数据现象；这时可使用与和缺失数据大小最接近的样本点的高分辨率影像资料来代替；如果最接近的样本点的高分辨率影像资料也不好获得，可用初略值来代替；

5、估计量及估计量方差的估计

假设起始值为r，当n为偶数时，相应系统样本的平均值为：

$\hat{Y}=N{\stackrel{\&OverBar;}{y}}_{\mathrm{sy}}=\frac{N}{n}(\underset{j=0}{\overset{\frac{n}{2}-1}{\mathrm{\&Sigma;}}}{y}_{r+\mathrm{jk}}+\underset{j=0}{\overset{\frac{n}{2}-1}{\mathrm{\&Sigma;}}}{y}_{N-r-\mathrm{jk}+1})$

$=\frac{N}{n}(y_r+y_{r+k}+\&CenterDot;\&CenterDot;\&CenterDot;+y_{r+(\frac{n}{2}-1)k}+y_{N-r-(\frac{n}{2}-1)k+1}+\&CenterDot;\&CenterDot;\&CenterDot;+y_{N-r-k+1}+y_{N-r+1})$

当n为奇数时，相应系统样本的平均值为：

$\hat{Y}=N{\stackrel{\&OverBar;}{y}}_{\mathrm{sy}}=\frac{N}{n}(\underset{j=0}{\overset{\frac{n}{2}-2}{\mathrm{\&Sigma;}}}{y}_{r+\mathrm{jk}}+y_{r+\frac{1}{2}(n-1)k}+\underset{j=0}{\overset{\frac{n}{2}-2}{\mathrm{\&Sigma;}}}{y}_{N-r-\mathrm{jk}+1})$

$=\frac{N}{n}(y_r+y_{r+k}+\&CenterDot;\&CenterDot;\&CenterDot;+y_{r+(\frac{n}{2}-2)k}+y_{r+\frac{1}{2}(n-1)k}+y_{N-r-(\frac{n}{2}-2)k+1}+\&CenterDot;\&CenterDot;\&CenterDot;+y_{N-r-k+1}+y_{N-r+1})$

如果把系统样本y₁、y₂、…、y_n看成从各层抽取两个单位的分层随机抽样，可采用以下办法进行方差估计；从第二个样本单元开始，每个样本单元与前一个样本单元组成一对，共n-1对，

的估计为：

$v\left(\hat{Y}\right)=N^{2}v\left({\stackrel{\&OverBar;}{y}}_{\mathrm{sy}}\right)=\frac{N(N-n)}{2n(n-1)}\underset{i=1}{\overset{n-1}{\mathrm{\&Sigma;}}}{(y_{i+1}-y_i)}^2$

6、检察是否需要再抽样本，如果需要，则按前述方法将补充的样本分配下去比较r的值是否达到了10％， $V={\left(\frac{d}{t}\right)}^2={\left(\frac{r\stackrel{\&OverBar;}{Y}}{t}\right)}^2;$ 如果没有，则需补充样本量，使得总的样本量为：

$n=\frac{s_1^2}{V}(1+\frac{2}{n_1})$

其中： $s_1^2=\frac{1}{2(n-1)}\underset{i=1}{\overset{n-1}{\mathrm{\&Sigma;}}}{(y_{i+1}-y_i)}^2,$ n₁是第2步所确定的样本量；

7、推断

总体总值的置信限为：

$\hat{Y}\&PlusMinus;u_{\frac{\mathrm{\&alpha;}}{2}}v\left(\hat{Y}\right)$