一种基于时序遥感影像的森林火烧迹地自动提取方法

摘要

本发明公开了一种新的基于时序遥感影像的森林火烧迹地自动提取方法，该方法包含三个步骤：第一步骤是选择森林火灾发生前和森林火灾发生后的无云遥感影像，计算两个时相遥感影像的各波段变化量和植被指数差值图像，将两者进行加权融合，得到差异图像；第二个步骤是对火灾发生前和火灾发生后的近红外波段进行线性拟合，求得拟合误差图像后对其进行二值化分类，将分类结果设为C-V模型的初始曲线；第三个步骤是将获取的差异图像和初始曲线作为C-V模型的输入，应用C-V模型对第一步获得的差异图像进行分割，实现火烧迹地信息的自动化提取。利用本发明能够很好实现应用遥感影像进行火烧迹地的快速与自动提取，具有很强的可操作性所提取的火烧迹地轮廓清晰，具有较高的检测精度，能够用于各种陆地卫星影像数据的森林区火烧迹地提取，实用性强。

说明书

技术领域

本发明涉及一种应用时序遥感影像进行森林火烧迹地自动提取的方法，适用于森林火灾监测领域。

背景技术

森林火灾是对森林破坏性最大的一种自然灾害，它不但破坏森林生态结构，甚至对人类的生存环境也 会造成重要影响。火烧迹地位置与覆盖范围的检测不仅能够有效得到森林的过火面积，从而评估火灾损失， 而且有助于林业部门有效开展火灾过后的植被恢复工作。遥感数据获取具有广视域、多波段、多时相和快 速等优势，为森林火烧迹地检测提供了一种有效的技术手段。目前利用遥感数据进行森林火烧迹地提取的 常用方法有监督分类法、基于火灾前后的归一化指数(Normalized Difference Vegetation，NDVI)差值法 和基于单时相遥感影像的图像分割法等。监督分类法是通过监督分类的方式实现火烧迹地的检测，由于选 择的训练样本很难普遍适合不同季节或不同空间覆盖范围的火烧迹地检测，因此训练样本往往需要不断增 加或更新调整，所以该方法比较依赖人为经验。基于火灾前后的NDVI差值法是通过计算火灾前遥感影像 与火灾后遥感影像的NDVI差值后通过阈值分析来提取火烧迹地，合适的阈值选取是提取精度的重要前提 条件，从而使得该方法在普适性应用方面具有在一定的局限性。基于单时相遥感影像的图像分割法，可以 对火灾的蔓延范围能够做出判断，但是却需要一定程度的人工干预，并不能实现火烧迹地提取流程的自动 化。所以当前这些方法普遍存在一定程度上的人工干预与人为经验依赖及自动化程度不够高等缺点，并不 能实现火烧迹地的自动与快速提取，影响了遥感数据在森林火灾信息提取方面的应用效率。因此为了减少 森林火灾的危害，如何通过应用多时相遥感数据快速而准确提取火烧迹地是亟待解决的问题。多时相遥感 影像的变化检测方法为火烧迹地的检测与提取了提供一种有效的途径。由于监督方法需要地面真实样本作 为训练数据，而适合不同地区与不同季节的真实地面样本数据存在很难获取的问题，因此应用非监督变化 检测方法来实现火烧迹地自动提取更能符合实际应用。在应用非监督变化检测实现火烧迹地自动地检测 时，根据不同时相的遥感影像获取差异图像是一个重要步骤。而当前的差异图像获取方法多数利用单一特 征，会导致检测的结果存在不够准确和不完整的缺陷。在获得差异图像后如何通过有效的方法分离出火烧 迹地是另一个重要的步骤，通常是通过分割差异图像的方法来实现，目前分割方法主要有阈值法、边缘提 取和分水岭分割法等。阈值法在实际应用中具有一定的局限性，主要体现在合适的阈值很难选择，而不同 的阈值会产生不同的分割效果。边缘提取算法，对边缘模糊的图像往往不能得到很好的结果。分水岭分割 算法是利用数学形态理论来实现图像分割，但容易出现过分割的问题。而森林火灾发生时由于受到风向或 地表覆盖条件的影响，在遥感图像上显示的火烧迹地往往存在边界模糊、边界凹陷等现象，因此在选择分 割方法时需要能够充分考虑到图像的全局信息，也能有效处理边缘模糊的问题。水平集Chan-Vese模型(简 称C-V模型)能够充分考虑图像的全局信息，也能有效处理边缘模糊的问题，然而该方法的初始曲线设定 对它的曲线演化速度和分割准确度存在重要的影响，不适当的初始曲线设置往往会降低火烧迹地的提取速 度与精度。

因此，为了提高应用遥感影像进行森林火烧迹地提取的速度与精度，实现森林火烧迹地的自动化提取， 需要一种实用与有效的火烧迹地自动提取方法。本发明针对以往方法的不足，针对稳定可靠的真实地面样 本数据难以获取的实际问题，和火烧迹地图像多数存在着边缘模糊和边界凹陷的特点，同时考虑到火烧迹 地提取的快速与自动化的实用性需求，提出了一种新的基于时序遥感影像的森林火烧迹地自动提取方法。 该方法通过将变化矢量分析法(Change Vector Analysis，CVA)和NDVI差值进行加权融合后产生差异图像， 并利用不同时相近红外波段影像进行线性拟合后，将拟合误差的二值化分类图作为新的初始曲线对传统水 平集C-V模型进行了改进，能够有效提高火烧迹地的提取速度和准确度，具有很强的可操作性。

发明内容

本发明的目的在于提供一种应用多时相遥感影像进行森林火烧迹地自动提取的方法，以克服现有方法 的缺点。利用本发明能够实现应用遥感影像进行森林火烧迹地的快速与准确提取，从而解决现有方法中存 在的前述问题，并提升遥感影像在森林火灾监测领域的应用水平。为了实现上述目的，本发明采用的技术 方案如下：

整个技术流程主要分为三个步骤，首先是通过对时相A(火灾发生前)和时相B(火灾发生后)的光 谱特征CVA和特征因子NDVI差值进行加权融合后得到差异图像；然后通过分析时相A的近红外波段Y和 时相B的近红外波段X之间的线性拟合误差，并通过K-means方法对误差图像进行二值化分类后得到改进 的C-V模型的初始曲线；最后将获取的差异图像和初始曲线作为C-V模型的输入，进行火烧迹地信息的自 动化提取，输出检测出的火烧迹地提取结果。

第一步、获得差异图像

具体为融合火灾发生前和火灾发生后的遥感影像的各波段光谱变化矢量CVA和NDVI差值从而获得差 异图像。

1-1)选择森林火灾发生前和森林火灾发生后的无云遥感图像，计算两个时相遥感图像的各波段的变 化量CVA；

1-2)分布计算火灾前和火灾后两个时相遥感影像的植被指数NDVI，然后计算两者的差值得到NDVI差 值图像；

1-3)将获得的变化矢量CVA和NDVI差异图像进行加权融合，得到差异图像。

第二步、C-V模型的初始曲线设置

具体为对火灾发生前与发生后遥感影像的近红外波段进行线性拟合，并对拟合误差图像进行二值化分 类，然后得到改进的C-V模型的初始曲线；

2-1)对火灾发生前遥感影像的近红外波段Y和火灾发生后遥感影像的近红外波段X进行线性拟合， 求得拟合误差图像；

2-2)应用K均值方法对拟合误差图像进行二值化分类，得到C-V模型的初始曲线。

第三步、火烧迹地提取

具体为将获取的差异图像和初始曲线作为C-V模型的输入，对获取的差异图像进行分割，实现火烧迹 地信息的自动化提取，最终输出检测出的火烧迹地提取结果。

附图说明

图1是本发明的主流程示意图

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。 应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。

如图1所示本发明公开了一种基于时序遥感影像的森林火烧迹地自动提取方法，包括以下步骤：

第一步、获得差异图像；

第二步、C-V模型的初始曲线设置；

第三步、火烧迹地提取。

以下对具体步骤进行详述：

差异图像获取是应用变化检测方法实现森林火烧迹地自动提取的首要步骤，目的是根据森林覆盖区域 火灾发生前后遥感影像的差异获得火烧迹地的初始信息。为了充分利用遥感数据的光谱特征，并考虑到森 林火灾导致了森林树木受到损害，叶绿素的数量减少从而导致植被指数下降的特点，因此本发明采用了火灾前 后遥感图像的光谱特征变化量CVA和指标指数NDVI差值进行加权融合来获取差异图像，在一定程度上增强了 火烧迹地识别的有效性。

光谱特征变化量CVA是通过对火灾发生前后遥感影像的各个波段数据进行差值运算，求得每个像素在 各个波段的变化向量。在每个变化向量中，变化的强度用变化向量的欧氏距离表示。设BVdate1和BVdate2 分别为图像大小为M×N个像元、具有n个波段的两幅经过几何配准的火灾发生前和发生后的遥感影像， CMpixel表示应用变化矢量分析方法得到的光谱特征变化量CVA：

${\mathrm{CM}}_{\mathrm{pixel}}=\underset{k=1}{\overset{n}{\Sigma }}{[{\mathrm{BV}}_{\mathrm{ijk}\left(\mathrm{datr}2\right)}-{\mathrm{BV}}_{\mathrm{ijk}\left(\mathrm{date}1\right)}]}^2$    (式1)

式(1)中CMpixel值越大表示像元值发生变化的可能性越大。

归一化植被指数NVDI被定义为近红外波段与红光波段反射率值之差和两者之和的比值。即由式(2) 计算得到。

$\mathrm{NDVI}=\frac{(\mathrm{NIR}-R)}{(\mathrm{NIR}+R)}$    (式2)

式(2)中NIR表示近红外波段的反射率值，R表示红光波段的反射率值。

将获得的光谱特征差异图像和归一化植被指数差异图像归一化处理后再进行加权融合，本文根据距离 运算中的卡方变换CST(Chi Square Transformation)[D’Addabbo A.，Satalino G.，Pasquariello G.， et al.2004.Three different unsupervised methods for change detection：an application. Proceedings of Geoscience and Remote Sensing Symposium2004，3：1980-1983]提出一种差异图像的加 权算法，是根据两种差异特征的标准差来获得光谱特征CVA和NDVI差异图像的加权系数。

$\mathrm{DI}=\frac{D_{\mathrm{CVA}}}{{\sigma }_C}+\frac{D_{\mathrm{NDVI}}}{{\sigma }_N}$    (式3)

式(3)中DI代表获得的最终差异图像(Different Image，DI)，D_CVA代表火灾发生前后的光谱特征差 异图像，D_NDVI是归一化植被指数NDVI差异图像。σc表示DCVA的标准差，σ_N表示D_NDVI的标准差。通过融 合两种差异图像得到最终的差异图像DI，并对DI进行归一化处理。

应用C-V模型[Chan T.，Vese L.，2001.Active contours without edges，IEEE Transactions on  Image Process，10：266-277]火烧迹地分割时，初始曲线的设置关系到火烧迹地提取的速度和准确性。本 发明采用通过火灾发生前后的遥感影像的近红外波段进行线性拟合后，并应用k均值二值化拟合误差的方 式来构造出水平集C-V模型的初始曲线。

假设Y为火灾发生前遥感图像的第k波段m维向量，X为t2时刻的遥感图像的第k波段m维向量。

Y^*=KX    (式4)

通过最小二乘法拟合来获得K的值，即要满足由上述约束条件可得

K=Σ₂₁Σ₁₁^-1    (式5)

式(5)中∑₂₁表示Y和X之间的协方差，∑₁₁表示X的方差。

通过马氏距离计算Y与拟合结果Y^*之间的误差，记为Error，如式(6)所示：

Error=(Y-Y^*)^T∑_ε(Y-Y^*)    (式6)

式(6)中∑ε表示的Y-Y^*方差，将误差Error中的较大值选定为变化区域。然后采用非监督的k-means 聚类方法对误差进行分析后，将误差分为变化与非变化两类。由此得到的二值化分类图作为水平集C-V模 型的初始曲线。

火烧迹地提取步骤是通过应用C-V模型分割差异图像实现。Mumford-Shah方法通过最小化能量函数找 到一条曲线C将图像X分为两类。能量函数如式(7)所示：

$F^{\mathrm{MS}}(I,C)=\underset{\Omega }{\int }{|X-I|}^2\mathrm{dxdy}+\lambda \underset{\Omega C}{\int }{|▿I|}^2\mathrm{dxdy}+\mu \left|C\right|$    (式7)

式(7)中I为由曲线C包围区域的灰度平均值，λ和μ都是正参数，X：Ω→R²。由于式(7)并不 是一个凸函数，所以在实际计算中很难求得最小值。Chan和Vese提出了针对这个问题的解决方法，即采 用水平集算法对Mumford-Shah方法进行简化，并引入了Heaviside函数。它将全部图像区域根据曲线C 划分成内部区Ω₁和外部区和Ω₂两个区域，使得在Ω₁内和Ω₂内的平均灰度值恰好能够反映出对象与背 景之间灰度平均值的差异。简化后的最小能量泛函如下：

$\left(\begin{array}{c}{F}^{\mathrm{MS}}(\Phi ,c_1,c_2)=\underset{\Omega }{\int }{|X-c_1|}^{2}H\left(\Phi (x,y)\right)\mathrm{dxdy}\\ +\underset{\Omega C}{\int }{|X-c_2|}^2(1-H\left(\Phi (x,y)\right))\mathrm{dxdy}\\ +\mu \underset{\Omega }{\int }|▿H\left(\Phi (x,y)\right)|\mathrm{dxdy}\end{array}\right)$    (式8)

式(8)中H(Φ(x，y))是Heaviside函数，

$H\left(z\right)=\left(\begin{array}{cc}1,& z>0\\ 0,& z\le 0\end{array}\right)$    (式9)

在应用中，Heaviside通常通过下式求得：

$H_{ϵ}\left(z\right)=\frac{1}{2}(1+\frac{2}{\pi }\arctan \left(\frac{z}{ϵ}\right))$    (式10)

式(8)中的c1和c2通过下式求得：

$c_1\left(\Phi \right)=\frac{{\int }_{\Omega }X(x,y)H\left(\Phi (x,y)\right)\mathrm{dxdy}}{{\int }_{\Omega }H\left(\Phi (x,y)\right)\mathrm{dxdy}}$    (式11)

$c_2\left(\Phi \right)=\frac{{\int }_{\Omega }X(x,y)(1-H\left(\Phi (x,y)\right))\mathrm{dxdy}}{{\int }_{\Omega }(1-H\left(\Phi (x,y)\right))\mathrm{dxdy}}$    (式12)

从式(11)和式(12)可以看出c1和c2分别对应曲线C内部和外部两个区域图像的灰度平均值。当 c1和c2确定后，利用变分法，式(8)可通过欧拉--拉格朗日方程得到如下偏微分方程：

$\frac{\mathrm{d\Phi }}{\mathrm{dt}}={\delta }_{ϵ}[\mu \mathrm{div}\left(\frac{▿\Phi }{|▿\Phi |}\right)-{(X-c_1)}^2+{(X-c_2)}^2]$    (式13)

式(13)中δε(Φ)是狄拉克函数

${\delta }_{ϵ}\left(\Phi \right)=\frac{d\left(H\left(\Phi \right)\right)}{\mathrm{d\Phi }}$    (式14)

根据C-V模型，将式(11)-(14)联立，求得稳定解即得到所要提取的分割曲线方程。

利用本发明的方法对美国陆地卫星Landsat遥感影像进了火烧迹地提取，并对火烧迹地检测精度与误 差进行统计，结果如表1所示。

从表1可以得知，本发明提出的方法将光谱信息CVA和归一化植被指数差值通过加权融合后得到差异 图像，从而有效地突出了火烧迹地的特征，在提取火烧迹地轮廓的过程中，根据火烧迹地导致植被损害的 原因，提出利用不同时相对应近红外波段的最小二乘法拟合误差二分类算法对传统C-V模型的初始曲线进 行改进，减少了C-V模型的曲线迭代次数，加快了检测速度，提高了检测结果的精度。实验结果表明本发 明提出的算法具有较高的检测精度，能够很好实现森林覆盖区域的火烧迹地的快速与自动提取。

表1火烧迹地检测精度及误差统计

遥感数据来源 漏检率(％) 虚警率(％) 正检率(％) Kappa系数 迭代次数 Landsat5TM 0.90 1.40 97.71 0.9684 349 Landsat5TM 1.06 0.32 98.63 0.9817 196 Landsat8OLI 0.65 2.47 96.88 0.9589 890