2018, Vol. 43引用本文 |
| 面向地理国情监测的地表覆盖变化信息提取方法 |  [PDF全文] |
利用遥感、地理信息系统等技术, 全面掌握地表覆盖情况, 并监测自然与人类社会经济活动引起的地表覆盖动态变化, 是地理国情监测重大专项的任务之一[1]。目前已提出的基于遥感技术的地表覆盖变化信息提取方法主要有基于像元和面向对象方法两类。其中, 基于像元方法又可分为分类后比较法和光谱差异法[2-3]两种。
高空间分辨率遥感影像的广泛应用, 促进了基于遥感影像的变化检测从传统像元级的检测向面向对象的检测转变。于信芳等以内蒙古草原区为例, 对这两类变化检测方法进行了比较, 得出面向对象的变化检测方法在总体精度、Kappa系数上都有明显的优越性[4]; 陆苗等研究了利用多尺度几何特征向量的面向对象变化检测方法, 并验证了方法的检测效果[5]。
在这些方法中, 光谱差异法能够确定变化发生区域、缩小分类范围、提高检测速度, 但不能同时获得具体的变化类型[6]; 分类后比较法能够确定变化类型, 但精度依赖于对不同时相影像分类的精度; 面向对象的变化检测技术目前尚不成熟, 尚未被广泛应用于工程项目中。
国家地理国情监测重大专项项目在普查阶段, 将建成普查成果本底数据库, 为监测提供公共、基础和标准的参考数据。本底数据库包含覆盖全国范围的地表覆盖数据、遥感影像数据, 分别为矢量和栅格格式。
基于丰富的本底数据库资源, 针对地理国情监测的需要, 本文从实用性、有效性角度, 提出了一种基于矢栅协同的地表覆盖变化信息提取方法, 结合矢量和栅格数据结构, 研究了试验区模型和数据处理过程, 验证了其可行性与应用潜力。
1 方法模型在地表覆盖发生变化的区域, 其遥感影像特征会发生变化, 通过不同时相的遥感影像与本底数据中的遥感影像进行对比计算, 便可获取变化发生区域。此时, 将变化发生区域矢量化, 与本底数据中的地表覆盖数据进行叠置分析, 便可获取变化类型信息, 得到变化类型分布范围以及面积等统计指标, 从而掌握土地空间利用和开发状态的变化。本文将这种采用矢量、栅格数据交互分析的方法称为矢栅协同方法。方法实现过程如图 1所示。
 |
| 图 1 矢栅协同方法的流程图 Figure 1 Flow Chart of Raster and Vector Data Synergy |
1.1 本底数据分析
地理国情监测中, 可用于地表覆盖变化信息提取的本底数据主要包括地表覆盖数据和遥感影像数据。地表覆盖数据为ShapeFile矢量格式, 几何类型为面, 属性项包含地理国情信息分类码(class code, CC)等; 遥感影像数据为Tiff或Erdas Img栅格格式, 使用TFW(tiff world file)格式记录影像坐标信息, 可扩展的标识语言(extensible markup language, XML)记录影像投影信息和元数据信息, 元数据信息包含影像获取时间、多光谱波段数量、多光谱波段名称等; 其中, 分幅遥感影像数据经过了融合、增强、彩色变换处理, 整景遥感影像数据未进行这些辐射处理。变化信息提取在自动化实现过程中, 适宜使用未作辐射处理的影像。
1.2 影像预处理影像预处理包含几何、辐射两个方面。
1) 正射纠正。采用经典数字微分纠正算法, 如共线方程法、RFM模型(rational function model)、RPC模型(rational polynomial coefficients)等, 对影像进行正射纠正。此处理须保证输出DOM (digital orthophoto map)的数学基础、分辨率与本底数据中的遥感影像一致。
2) 辐射量化统一。在不同时相的遥感影像与本底数据中的遥感影像辐射量化比特数(如WorldView-2多光谱影像为16 bits)不同时, 需进行辐射分辨率处理, 使经过辐射处理后的比特数与本底数据中的遥感影像统一, 以保证变化检测指标各因子的计算在同一尺度、同一量级上进行。
3) 彩色变换。为了目视判读方便, 可对不同时相的遥感影像进行彩色变换, 生成一景新影像, 使其接近自然真彩色, 或突出部分地表覆盖信息, 以增强判读效果。
1.3 遥感指数计算遥感指数本身便是对特定的、感兴趣的地表覆盖类型信息的增强, 因此, 在变化指标设计时, 可考虑顾及遥感指数因子的算法, 将遥感指数差值作为变化检测指标因子之一[7]。
经典的遥感指数, 如归一化植被指数(normalized difference vegetation index, NDVI)[8-9]、归一化水体指数(normalized difference water index, NDWI)等, 均进行了归一化处理, 其值域范围为[-1, 1], 由此得到的遥感指数差值也在较小的区间范围。但较小的遥感指数差值可能表征着较大的变化信息。因此, 在计算过程中, 需将遥感指数差值进行如下变换:
1) 遥感指数差值正值化。将遥感指数差值拉伸到非负区间, 使用正值化后的遥感指数差值表征变化量。
2) 量化范围拉伸。由于不同时相的遥感影像与本底数据中遥感影像的光谱差异值跟辐射量化比特数以及色彩饱和度相关, 与遥感指数差值不在同一量化级。因此, 需对遥感指数差值进行拉伸处理, 使其量化范围与影像的光谱差异值处于相同量化级。其具体方法是:对影像的光谱差异值进行直方图统计; 然后采用直方图匹配或线性拉伸算法, 对遥感指数差值进行拉伸。
3) 因子权重设置。根据地表覆盖类型, 可对遥感指数因子设置不同的权重, 使其发挥不同权重的作用[10]。
1.4 变化检测指标设计变化发生区域范围与变化检测指标的选择以及阈值的设置有关。现有的地表覆盖变化检测方法主要使用单一的距离指数或相似性指数等指标来比较不同时相遥感影像间的特征差异, 然后根据确定的域值判断像元对应区域地表覆盖类型是否发生变化。常见的指数有欧氏距离、绝对距离、夹角余弦、相关系数[11-12]。本文选择经典欧氏距离作为变化检测指标, 其计算公式为:
| $ {d_i} = \sqrt {{\rm{ }}\sum\limits_{p = 1}^n {{{({x_{{t_1}}}-{x_{{t_2}}})}^2}} {\rm{ }}} $ | (1) |
式中, di表示第i个像元的变化检测指标值; xt1、xt2分别表示第i个像元第p波段在时相t1、t2的特征值; n表示波段数量。
1.5 变化发生区域与变化类型确定计算出变化检测指标值后, 需设置阈值和判定条件[13], 确定像元是否发生变化, 从而检测出发生变化的可疑像元。经典算法是对变化检测指标值进行直方图统计, 利用加权平均值与标准方差, 定义阈值的计算公式[14]为:
| $ c = \bar d + {a_i}{\sigma _i} $ | (2) |
式中, c表示阈值; d表示变化检测指标值的加权平均值; σi表示标准方差; ai表示调整参数, 其值域范围为[0, 1.5][15-18]。通过调整参数, 使确定的阈值适于研究区域内不同地表覆盖类型变化信息的确定。
变化检测指标值大于等于阈值的像元, 为发生变化的可疑像元。在计算变化检测指标值的过程中, 为防止计算值溢出, 需根据计算值的值域范围, 定义比遥感影像辐射量化比特数更高阶的变量, 或求取均方根, 使其值域范围在定义的变量存储范围内, 式(1)采用的是求取均方根的方法。
2 试验与结果分析试验区位于中国南部, 面积为8.0 km2, 属于平地地形, 高程范围为0~95.0 m, 地表覆盖类型多样, 分布有耕地、林地、草地、水域、房屋建筑、道路等类型, 人工地表覆盖比较密集, 属于人类活动比较频繁的地区。试验数据为两期RapidEye遥感影像, 包含蓝、绿、红、红边、近红外等5个波段, 获取时间分别为2011-12-24和2012-11-06, 时相较一致, 影像质量较好, 无云覆盖; 将2011-12-24影像模拟为遥感影像本底数据, 将2012-11-06影像作为不同时相遥感影像。同时, 收集了该地区30 m分辨率DEM数据, 作为影像纠正、分类辅助数据。
2.1 数据处理采用RPC模型对两期影像进行了正射纠正, 对两期影像进行了精配准处理; 按照地理国情监测地表覆盖分类体系, 对2011-12-24影像进行了分类, 并形成了地表覆盖本底数据。
地理国情监测地表覆盖分类体系包含10个一级类(耕地、园地、林地、草地、房屋建筑(区)、道路、构筑物、人工堆掘地、荒漠与裸露地表、水域), 87个三级类, 本试验区共分类出9个一级类, 20个三级类, 分类结果如图 2(c), 一级类面积统计如表 1。
 |
| 图 2 遥感影像与地表覆盖分类图 Figure 2 Classification Result of Images and Corresponding Land Cover |
| 表 1 地表覆盖分类面积统计 Table 1 Area Statistics of Land Cover for Each Class |
 |
2.2 结果与分析
利用式(1)计算出变化检测指标值, 得到变化强度图及其统计直方图, 如图 3(a)、3(b)。从图 3(b)可以得到, 加权平均值为370.5, 标准方差为345.6, 这里设置调整参数为1.1, 计算阈值为750.7, 从而提取到变化发生的可疑像元分布范围, 如图 3(c) (白色区域), 也即提取到地表覆盖变化发生可疑区域。
 |
| 图 3 变化信息提取 Figure 3 Land Cover Change Information Extraction |
将变化发生可疑区域与地表覆盖本底数据进行叠置分析, 主要包括空间叠加、逻辑交、属性信息投射等, 并进行噪声去除处理, 最终提取到变化发生区域、变化类型信息。变化发生的具体情况以及面积统计信息如表 2。
| 表 2 地表覆盖变化信息统计/m2 Table 2 Area Statistics of Land Cover Change for Each Class/m2 |
 |
1) 图 2、图 3反映了试验区在监测周期内(2011-12-24至2012-11-06)的土地开发与利用状态及变化情况。
2) 从表 2可以看出, 试验区地表覆盖变化范围面积合计649 359.2 m2, 占试验区总面积8.1%。
结果表明, 9个一级类中, 有7个一级类发生了变化。其中, 变化发生面积最大的是林地, 达到了459 417.8 m2; 大部分变化成了构筑物(三级类为其他硬化地表), 占林地变化总面积的97.9%。通过对影像的对比分析, 这种变化是由于试验区在开发道路、房屋建筑(区)等过程中, 占用了较多的林地区域所致。
变化发生面积较大的是人工堆掘地, 大部分变化成了道路、房屋建筑(区), 分别占人工堆掘地变化总面积的46.9%、46.2%。部分荒漠与裸露地表变化成了水域, 是由于水域发生了季节性变化, 部分水域较浅地区的泥土地表被水域覆盖。在发生变化的可疑区域中, 有部分构筑物变成了新的构筑物, 通过对影像的对比分析, 判断是这些区域同物异谱导致的。因此, 没有计入最终的变化信息统计中。地表覆盖本底数据中的耕地、水域没有发生变化。
3) 图 4标注了部分变化发生显著区域(带标注号的多边形范围), 叠加影像为2011-12-24影像。从标注1所在的区域可以看到, 为了建设道路, 在2011-12-24属于林地的地块, 被用作道路建筑工地, 从自然地表覆盖变化为人工地表覆盖。标注2处的路段在2012-11-06建设完成, 而在2011-12-24还处于分段建设阶段, 整条道路的贯通, 只剩中间部分路段。从标注3可以看到, 一块在2011-12-24属于土质硬化地表的地块, 在2012-11-06建成一块标准露天体育场, 由跑道和草坪构成。标注4反映了水域在不同季节的变化情况, 由于不同位置水量、水深的差异, 水域在遥感影像上的光谱特征有所差异, 裸露出来的泥土地表的范围也不同。
 |
| 图 4 部分变化发生显著区域 Figure 4 Some Regions of Land Cover Change Obviously |
3 结束语
本文提出的地表覆盖变化信息提取方法, 有效利用了地理国情普查构建的地表覆盖本底数据库资源, 通过空间叠置分析对检测出的变化发生区域进行变化类型的确定, 减少了分类工作量, 克服了单一方法不能同时获得变化类型且过度依赖不同时相影像分类精度的缺点, 能够为地理国情监测提供方法参考。
| [1] | 陈俊勇. 地理国情监测的学习札记[J]. 测绘学报, 2012, 41(5): 633–635 |
| [2] | Lu D, Mausel P, Brondízio E, et al. Change Detection Techniques[J]. International Journal of Remote Sensing, 2004, 25: 2 365–2 407 DOI: 10.1080/0143116031000139863 |
| [3] | Civco D L, Hurd J D, Wilson E H, et al. A Comparison of Land Use and Land Cover Change Detection Methods[C]. ASPRS-ACSM Annual Conference and FIG ⅩⅩⅡ Congress, Washington, D C, USA, 2002 |
| [4] | 于信芳, 罗一英, 庄大方, 等. 土地覆盖变化检测方法比较—以内蒙古草原区为例[J]. 生态学报, 2014, 34(24): 7 192–7 201 |
| [5] | 陆苗, 梅洋, 赵勇, 等. 利用多尺度几何特征向量的变化检测方法[J]. 武汉大学学报·信息科学版, 2015, 40(5): 623–627 |
| [6] | 张振龙, 曾志远, 李硕, 等. 遥感变化检测方法研究综述[J]. 遥感信息, 2005, (5): 64–66 |
| [7] | Cakir H I, Khorram S, Nelson S A C. Correspondence Analysis for Detecting Land Cover Change[J]. Remote Sensing of Environment, 2006, 102(3): 306–317 |
| [8] | 黄昌狄, 徐芳. 基于真彩色高分辨遥感影像稀疏植被覆盖检测[J]. 测绘地理信息, 2016, 41(5): 42–46 |
| [9] | 徐芳, 陈联, 唐雅梅, 等. 独立坐标下卫星影像立体测图技术研究[J]. 测绘地理信息, 2013, 38(4): 27–30 |
| [10] | Yue T X, Chen S P, Xu B, et al. A Curve-Theorem Based Approach for Change Detection and Its Application to Yellow River Delta[J]. International Journal of Remote Sensing, 2002, 23(11): 2 283–2 292 DOI: 10.1080/01431160110106041 |
| [11] | 李月臣, 陈晋, 宫鹏, 等. 基于NDVI时间序列数据的土地覆盖变化检测指标设计[J]. 应用基础与工程科学学报, 2005, 13(3): 261–275 |
| [12] | 刘春霞, 李月臣, 杨华, 等. 三峡库区重庆段生态与环境敏感性综合评价[J]. 地理学报, 2016, 66(5): 631–642 |
| [13] | Rosin P. Thresholding for Change Detection[J]. Computer Vision and Image Understanding, 2002, 86(2): 79–95 DOI: 10.1006/cviu.2002.0960 |
| [14] | Rosin P, Ioannidis E. Evaluation of Global Image Thresholding for Change Detection[J]. Pattern Recognition Letters, 2003, 24(14): 2 345–2 356 DOI: 10.1016/S0167-8655(03)00060-6 |
| [15] | Morisette J T, Khorram S. Accuracy Assessment Curves for Satellite-Based Change Detection[J]. Photogrammetric Engineering & Remote Sensing, 2000, 66(7): 875–880 |
| [16] | Tian J, Nielsen A A, Reinartz P. Improving Change Detection in Forest Areas Based on Stereo Panchro matic Imagery Using Kernel MNF[J]. IEEE Transactions on Geoscience & Remote Sensing, 2014, 52(11): 7 130–7 139 |
| [17] | Castle C, Finomore V, Simpson B, et al. Evaluation of Spatial Audio for Improving Change Detection on Large Screen Displays[DB/OL]. Proceedings of the Human Factors & Ergonomics Society Annual Meeting. [2016-06-18]. http://journals.sagepub.com/doi/abs/10.1177/1071181312561330 |
| [18] | Nielsen A A. Kernel Maximum Autocorrelation Factor and Minimum Noise Fraction Transformations[J]. IEEE Transactions on Image Processing, 2011, 20(3): 612–624 DOI: 10.1109/TIP.2010.2076296 |