CN106840409A - 一种基于modis的森林火灾火点检测方法 - Google Patents

Abstract

本发明公开了一种基于MODIS的森林火灾火点检测方法，先对来自MODIS的若干天的MOD02及MOD03数据构建时间序列数据，得到该若干天的NDVI和亮温的时间序列日值数据，然后采用时间序列断点识别算法识别NDVI和亮温的时间序列日值数据中的断点，将识别出的断点作为潜在火点；采用空间上下文分析法从潜在火点中判别出相对火点；从自主模块化传感器数据中提取真实火点，将相对火点与真实火点进行对比，验证相对火点的正确性和遗漏率。本发明利用BFAST算法检测时间序列数据断点的方法判别潜在火点，而不是固定阈值确定潜在火点，因此本发明算法对于火点判别更具有自适应性和灵活性，克服了固定阈值判别潜在火点无法满足区域性差异的不足，提高了火点判对率。

Description

一种基于MODIS的森林火灾火点检测方法

技术领域

本发明涉及一种基于MODIS的森林火灾火点检测方法。

背景技术

森林是一种主要的自然资源，在维持环境平衡中发挥着重要作用。而频发的森林火灾会直接导致森林退化，严重威胁地表植被和野生动物的生命健康。鉴于森林火灾的严重性，及时、准确地检测森林火点对减少灾害损失、保护森林资源和野生动物、掌握气候变化规律、维护生态平衡等具有重要意义。近年来，覆盖范围广、时空分辨率高的遥感手段在监测森林火灾中得到广泛的运用。而MODIS传感器扫描范围宽、属于太阳同步近极地轨道，空间分辨率与光谱分辨率适宜，在火灾监测方面表现良好的适用性和独特的优势。

随着MODIS火点检测产品的发展，相应算法不断丰富。目前林火火点检测算法主要有光谱法、阈值法、空间上下文分析法、多时相法等。光谱法借助火点与其他地物的光谱差异，增强图像信息，目视定性地了解火点大致位置。阈值法是根据区域和季节的先验信息，设定火点判别物理量的阈值，从而确定火点；Kaufman、Kennedy、Arino、Melinotte和覃先林等人分别引入亮温、反射率、反照率、NDVI阈值、亮温-植被指数作为火点判别物理量，Lee和Tag技术依据像元背景亮温和地表温度的差异对通道亮温进行了校正；阈值方法简单省时高效，但很难适用于所有的气候和环境状况。针对固定阈值法的不足，Giglio提出了上下文算法，并发展成为了MODIS V4算法；基于此，Wang引入烟羽掩膜确定潜在火点区域、Wang调整了潜在火点判别阈值中MIR和TIR通道的亮温差、彭光雄等采用了烟羽掩膜和自适应背景窗口技术、HEL等提出了消除太阳污染的方法、Ying Li提出特征空间马氏距离法、肖霞提出类间方差和烟羽掩膜法，使得MODIS V4算法日益成熟。与此同时，多时相法得到迅速发展，Keiji和Kushida基于多时相MODIS V4产品构建双谱阈值的随机模型比MODIS V4精度更高；Movaghati S利用自动Agent算法检测区域火点，对面积小、温度低火点检测比MODIS V4更有优势。总的来说，上下文分析法设定自适应背景窗口内目标与背景的差异度判识火点，提高了方法的适用范围和准确度，但通过固定阈值确定潜在火点的方法灵活性差，区域性差异考虑不足；多时相火点检测方法克服了不同区域环境阈值设置问题，更具有自适应性，能在早期发现火点，但要对历史数据进行统计分析和变化阈值设定。

发明内容

本发明的目的是解决目前森林火点检测方法中，上下文分析法需通过固定阈值确定潜在火点，灵活性差，区域性差异考虑不足；多时相火点检测方法要对历史数据进行统计分析和变化阈值设定，使用不便的技术问题。

为实现以上发明目的，本发明提供一种基于MODIS的森林火灾火点检测方法，包括如下步骤：

(1)识别潜在火点：先对来自MODIS的若干天的MOD02及MOD03数据构建时间序列数据，得到该若干天的NDVI和亮温的时间序列日值数据，然后采用时间序列断点识别算法识别所述NDVI和亮温的时间序列日值数据中的断点，将识别出的断点作为潜在火点；

(2)判别相对火点：基于所述潜在火点及其背景像元，进行空间上下文分析，将识别出的火点作为相对火点，所述相对火点为真实火点。

进一步地，所述步骤(1)中，在对所述MOD02及MOD03数据构建时间序列数据之前，先对所述MOD02及MOD03数据进行预处理。

进一步地，所述预处理包括几何校正、数据有效性检查、亮温反演和NDVI反演。

进一步地，所述几何校正包括去除“蝴蝶结”效应、投影转换和图像规则裁剪。

进一步地，所述数据有效性检查包括对经几何校正后的数据进行重新赋值以及数据拼接处理。

进一步地，步骤(1)中所述采用时间序列断点识别算法识别所述NDVI和亮温的时间序列日值数据中的断点的步骤如下：

(a1)利用可加性分解模型迭代拟合趋势组分和季节组分，表示为：

Y_t＝T_t+S_t+e_t (t＝1,...,n)；

其中：t为观测时间，Y_t为t时刻观测值；T_t为趋势组分；S_t为季节组分；e_t为噪声组分，是观测值中除去趋势组分和季节成分后的剩余部分；

设断点为 则趋势组分T_t的线性模型表示为：

其中：j为趋势断点所在位置，j＝1,...,m；α_j为截距，β_j为斜率，用来反演突变的幅度和断点前后渐变的斜率，对于断点t，若满足则断点t突变的幅度表示为：

Magnitude＝(α_j-1-α_j)+(β_j-1-β_j)^t；

其中，β_j-1和β_j分别为断点前、后的渐变斜率；

设断点为 则季节组分S_t的谐波模型表示为

其中，j为季节断点所在位置，j＝1,...,p；k为谐波项数，f为频率；γ_j,k＝a_j,kcos(δ_j,k)，θ_j,k＝a_j,ksin(δ_j,k)，a_j,k和δ_j,k分别为频率f/k的振幅和相位：和δ_j,k＝tan^-1(θ_j,k/γ_j,k)；

(a2)断点监测：

利用基于最小二乘法的移动求和检验是否存在断点，利用最小贝叶斯信息准则确定断点的最优数量，利用最小残差平方和估算断点的位置。

进一步地，所述步骤(a2)的步骤如下：

(a21)如果最小二乘法的移动求和检验表明趋势组分T_t中有断点，运用最小二乘法拟合季节性因素调整后的数据得到趋势断点的数量和位置；

(a22)由基于M-估计的鲁棒回归估算得到趋势组分系数α_j和β_j，基于步骤(a1)中所述公式Magnitude＝(α_j-1-α_j)+(β_j-1-β_j)t，趋势组分估计值表示为：

(a23)如果最小二乘法的移动求和检验表明季节组分S_t中有断点，运用最小二乘法拟合趋势性因素调整后的数据得到季节断点的数量和位置；

(a24)由基于M-估计的鲁棒回归估算得到季节组分系数θ_j,k和γ_j,k，基于谐波模型，季节组分估计值表示为：

(a25)以上步骤不断地迭代，直到断点的数量和位置不再变化。

进一步地，步骤(2)中所述采用空间上下文分析法从所述潜在火点中判别出相对火点的步骤如下：

(b1)云体掩膜：

利用高级超高分辨率辐射计数据监测体，白天满足下式的像元为云体像元：

ρ₁+ρ₂＞0.9orT₃₂＜265Kor(ρ₁+ρ₂＞0.7and T₃₂＜285K)；

其中，ρ₁、ρ₂、T₃₂分别为MODIS波段1的反射率、波段2的反射率、波段32的亮温；

(b2)确定有效背景像元：

以潜在火点为中心，向左、右、上、下分别延长1-5个像元，构建5×5、7×7……21×21大小的窗口，直到有效背景像元个数占窗口总像元数量的25％以上停止搜索，来确定窗口大小；

(b3)火点判别：

在获取所述有效背景像元之后，采用上下文空间统计分析方法对潜在火点进行判别，判别条件如下式，满足该条件则所述有效背景像元为火点像元，

其中，为有效背景像元22和31波段亮温差的平均值，为有效背景像元22和31波段亮温差的平均绝对偏差；和分别为有效背景像元22和31波段的亮温差；和分别为有效背景像元22和31波段亮温的平均绝对偏差；为背景火点像元22波段亮温的平均绝对偏差。

进一步地，所述步骤(2)后面还包括：

从自主模块化传感器数据中提取真实火点，将步骤(2)中所述相对火点与所述真实火点进行对比，验证所述相对火点的正确性和遗漏率。

本发明综合运用时间序列分析法和空间上下文分析法，利用时间序列断点识别算法(BFAST)检测归一化植被指数(NDVI)和亮温的时间序列日值数据断点，通过断点识别潜在火点，进而利用空间上下文统计分析确定相对火点；将检测结果与MODIS V4结果进行对比分析，并采用自主模块化传感器(AMS)数据的真实火点验证结果的准确性。相对于MODISV4火点提取算法，本发明提高了火点判对率，这是因为本发明利用BFAST算法检测时间序列数据断点的方法判别潜在火点，而不是固定阈值确定潜在火点，因此本发明算法对于火点判别更具有自适应性和灵活性，克服了固定阈值判别潜在火点无法满足区域性差异的不足。

附图说明

图1是圣迭戈林区分布图。

图2a是美国加利福尼亚州地理位置图。

图2b是加利福尼亚州破坏性最大的林火分布图。

图2c是加利福尼亚州伤亡最大的林火分布图。

图2d是加利福尼亚州规模最大的林火分布图。

图3是样本点位置图。

图4a是1号样本点NDVI时间序列图。

图4b是1号样本点亮温时间序列图。

图4c是2号样本点NDVI时间序列图。

图4d是2号样本点亮温时间序列图。

图5a是1号样本点NDVI断点检测结果。

图5b是1号样本点亮温断点检测结果。

图5c是2号样本点NDVI断点检测结果。

图5d是2号样本点亮温断点检测结果。

图6a是AMS的火点分布图。

图6b是AMS的3.7μm波段亮温分布图。

图6c是AMS的3.7μm和10.76μm波段亮温散点图。

图6d是AMS的10.76μm波段亮温分布图。

图7a是1号样本点的AMS火点分布图。

图7b是2号样本点的AMS火点分布图。

图8a是MODIS V4算法的火点提取结果。

图8b是本发明的火点提取结果。

图9a是2007年10月24日圣迭戈火点采用MODIS V4算法的火点提取结果。

图9b是2007年10月24日圣迭戈火点采用本发明的火点提取结果。

图10a是确定有效背景像元时的5×5大小的窗口示意图。

图10b是确定有效背景像元时的7×7大小的窗口示意图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明。

如图1和2a所示，本实施例选取美国圣迭戈作为样本进行研究。圣迭戈(SanDiego)是美国加利福尼亚州的一个太平洋沿岸城市，位于加利福尼亚州的西南端，南距墨西哥边境20公里，西经117.6°-116.0°，北纬32.5°-33.6°。人口134万，面积963km2，是加利福尼亚州的第二大城市。迭戈市年平均气温约为华氏70度，以温暖的气候和多处的沙滩著名，属温带大陆性干旱气候，阳光充足，降雨较少，年均降水量约为25cm。

从图2b-2d可以看出，圣迭戈发生破坏严重、严重伤亡、大规模火灾的次数分别为4、5、3，属于火灾较严重的地区。参考MCD12Q1的第一种分类体系的分类结果：1常绿针叶林，2常绿阔叶林，3落叶针叶林，4落叶阔叶林，5混交林，6封闭灌丛，7开放灌丛，8多树草原，共8类像元作为本发明的林区范围。

本发明主要使用的数据包括实验数据和验证数据。

(1)实验数据为2002年1月1日到2009年9月30日共2828天的MOD02和MOD03数据。MOD02即MODIS原始数据MODIS L1B经过定位和辐射定标处理的数据之一，分别率为1KM。包含：反射率(Reflectance)、辐亮度值(Radiance)和发射率(Emissive)三种类型数据，不同分辨率与波段组合的四个波段子集。主要的预处理包括几何校正、数据有效性检查、亮温反演和NDVI反演；

具体地，几何校正包括去除“蝴蝶结”效应、投影转换和图像规则裁剪，基于java语言，使用MODIS Swath格式数据的简易软件平台MRTSwath去除“蝴蝶结”效应，进行投影转换、图像规则裁剪等。需要注意的是，该软件需要MODIS L2数据对应的几何定位数据MOD03。

数据有效性检查：为保证每个像元每天有且仅有一个时间序列数据，需对hdf文件经MRTSwath几何校正生成的TIF文件进行重新赋值以及数据拼接处理。

亮温反演，根据反解Planck函数求解亮温的原理，由辐射率得出MOD021KM/MYD021KM定标辐射热红外波段(20-25，27-36)的亮温T_λ，计算公式为：

该公式用来反演亮温，直接下载的影像MODIS02数据提供波长、波段辐射亮度、热辐射偏移量、拉伸量，已知常数参量即可求出时间序列数据的亮温。其中：π＝3.14159；C₁＝3.7418×10^-16W·m²；C₂＝14388μm·K；λ为波长(m)；L_λ为波段辐射亮度(W·m^-2·sr^-1)；L_offset_λ和L_scales_λ分别为热辐射偏移量和拉伸量，均可从HDF文件属性中获取。

NDVI反演，归一化植被指数(Normalized Difference Vegetation Index,NDVI)表示为：

或

该公式基于MODIS02数据来计算时间序列的NDVI。其中：ρ_NIR和ρ_RED分别为近红外和红光波段的反射率，DN_NIR和DN_RED分别为近红外波段和红光波段的DN(Digital Number是遥感影像像元亮度值)值。MODIS的第一、二波段分别为红光和近红外波段，可以通过MOD021KM/MYD021KM数据集250m重采样到1000m的反射率波段得到。

(2)验证数据为AMS(自主模块化传感器)数据，该数据具备相关的光谱测量能力，可为火灾温度、面积反演以及高分辨率火点辐射功率(FRP)的验证提供数据。为保证评估的准确性，要求MODIS和AMS过境时间差控制在15-17分钟内，本发明选取Ikhana于2007年10月24日19:11-19:18获取的圣迭戈AMS数据作为验证数据。

本发明采用如下方法来检测火点：

1时间序列分析

时间序列断点识别算法BFAST(Breaks For Additive Seasonal and Trend)是Verbesselt等提出的时间序列变化检测模型，目前主要用于NDVI产品的断点识别来研究物候扰动。本发明利用BFAST算法识别的NDVI和亮温两个时间序列观测数据断点作为潜在火点。

1.1分解模型

BFAST算法将长时间序列影像通过迭代方式进行逐像元的曲线分离，过滤出趋势组分(即年际间变化，由气候、土地管理、人类活动引起的变化)、季节组分(即年内周期性变化)及噪声变化(主要指云、雨、雪、仪器噪声等引起的瞬时变化)。

具体来说，利用可加性分解模型迭代拟合趋势组分和季节组分，表示为：

Y_t＝T_t+S_t+e_t (t＝1,...,n)； (3)

其中：t为观测时间，Y_t为t时刻观测值；T_t为趋势组分；S_t为季节组分；e_t为噪声组分，是观测值中除去趋势组分和季节成分剩余部分。

(1)趋势组分T_t：描述为分段直线函数。设断点为 则线性模型表示为：

其中：j为趋势断点所在位置，j＝1,...,m；截距α_j和斜率β_j可用来反演突变的幅度和断点前后渐变的斜率。对于断点t，若满足则断点t突变的幅度表示为：

Magnitude＝(α_j-1-α_j)+(β_j-1-β_j)t； (5)

其中，β_j-1和β_j分别为断点前后的渐变斜率。

(2)季节组分S_t：相邻断点间的S_t是固定的，但不同断点间S_t可能不同，且季节断点与趋势断点可能发生在不同时间。设断点为 季节组分S_t为谐波模型，表示为：

其中，j为季节断点所在位置，j＝1,...,p；k为谐波项数；f为频率(如：16天时间序列的年观测值f＝23)；γ_j,k＝a_j,kcos(δ_j,k)，θ_j,k＝a_j,ksin(δ_j,k)，a_j,k和δ_j,k分别为频率f/k的振幅和相位：和δ_j,k＝tan^-1(θ_j,k/γ_j,k)。

谐波模型有如下优点：模型对短期数据的变化和固有噪声(云层和大气散射)不敏感；必需的观测值更少，提高了算法速度和效率；拟合参数(a_j,δ_j)更容易用来描述物候变化。

1.2断点监测

利用基于最小二乘法的移动求和(OLS-MOSUM)检验测试是否存在断点，利用最小贝叶斯信息准则(BIC)确定断点的最优数量，利用最小残差平方和(MRSS)估算断点的位置。具体步骤如下：

(1)如果OLS-MOSUM测试表明趋势组分T_t中有断点，运用OLS拟合季节性因素调整后的数据得到趋势断点的数量和位置。

(2)考虑到潜在异常值，趋势组分系数α_j和β_j由基于M-估计的鲁棒回归估算得到，基于公式(5)，趋势组分估计值表示为：

(3)如果OLS-MOSUM测试表明季节组分S_t中有断点，运用OLS拟合趋势性因素调整后的数据得到季节断点的数量和位置。

(4)季节组分系数θ_j,k和γ_j,k由基于M-估计的鲁棒回归估算得到，基于谐波模型，季节组分估计值表示为：

(5)以上步骤不断地迭代，直到断点的数量和位置不再变化。

2空间上下文分析

BFAST算法检测的断点不一定为林火，因此，本发明将潜在火点进行空间上下文分析判定其是否为火点，主要包括云体掩膜、确定有效背景像元和火点判别三部分。

2.1云体掩膜

云体掩膜借鉴IGBP(International Geosphere Biosphere Program)，即利用AVHRR(高级超高分辨率辐射计)数据监测全球火灾的云监测技术，白天满足式(9)的像元认为是云体像元。

ρ₁+ρ₂＞0.9or T₃₂＜265Kor(ρ₁+ρ₂＞0.7and T₃₂＜285K)； (9)

其中，ρ₁、ρ₂、T₃₂分别为MODIS波段1反射率、波段2反射率、波段32的亮温。

2.2确定有效背景像元

有效背景像元指以潜在火点为中心并且满足下列四个条件的像元：①该像元数据有效；②该像元为非云像元；③该像元为陆地像元；④该像元为非背景火点像元。背景火点像元指白天T₂₂＞325K且T₂₂-T₃₁>20K的像元(T₂₂为22波段的亮温)。以潜在火点为中心，向左右上下分别延长1-5个像元，即构建5×5、7×7……21×21大小的窗口，如图10a和10b所示，直到有效背景像元个数占窗口总像元数量的25％以上即停止搜索，来确定窗口大小。

2.3火点判别

在获取有效背景像元之后，采用上下文空间统计分析方法对潜在火点进行判别，判别条件如式(10)，满足条件则该像元为火点像元。

其中，为有效背景像元22和31波段亮温差的平均值，为有效背景像元22和31波段亮温差的平均绝对偏差；和分别为有效背景像元22和31波段的亮温差；和分别为有效背景像元22和31波段亮温的平均绝对偏差；为背景火点像元22波段亮温的平均绝对偏差。

结果分析

1样本的火点判别

选取两个MODIS像元作为样本点对本发明算法进行分析，如图3为1、2号样本点。

1.1构建时间序列数据

如图4a-4d所示，利用2002年1月1日到2009年9月30日共2828天的MOD02和MOD03数据，经过几何校正、转投影、拼接，提取了从Band_250m_RefSB和Band_1KM_Emissive两个波段子集，分别提取MODIS 1、2和22波段的Scaled Integer数据，利用参数offset和scalefactor得到时间序列影像1、2波段的反射率和22波段的辐亮度。在matlab中逐像元读取1、2波段的反射率和22波段的辐亮度的时间序列数据，对有效范围之外的数据进行时间序列线性插值，保证像元时间序列数据质量上的有效性，则每个像元将得到2828天的NDVI和22波段亮温时间序列日值数据，分别保存在以各自行列号命名的txt文件中。

1.2潜在火点的识别

本发明基于R语言，使用Verbesselt团队编写BFAST软件包，根据bfast函数语法输入潜在断点的最小间隔比例、季节性模型、最大断点个数(本发明设为4)等参数来进行断点检测。图5a和5c分别为1号样本点和2号样本点NDVI时间序列断点提取结果，图5b和5d分别为1号样本点和2号样本点22通道亮温时间序列断点提取结果，表1为BFAST算法潜在火点的识别结果提取信息统计表，从图表中可以得知：(1)1号点在2004年第270天(降)、2006年第64天(降)，2007年第297天(降)NDVI数据趋势组分有断点，2号点2005第168天(升)、2006第290天(升)、2008第36天(升)NDVI数据趋势组分存在断点。如果NDVI断点为“升”表明NDVI突然增大，可能因为造林绿化；如果NDVI断点为“降”表明NDVI突然减小，可能原因有火灾、森林砍伐、洪灾等。(2)1号点在2004年第290天(降)、2007年第297天(升)亮温数据趋势组分有断点，在2007年第298天亮温数据季节组分有断点；2号点在2003第100天(降)、2004第291天(降)、2006第58天(降)、2008年第37天(升)亮温数据趋势组分有断点。如果亮温断点为“升”表明亮温突然增大，可能因为林火灾害、新增加的人工热源；如果亮温断点为“降”表明亮温突然下降，可能因为植被增加。

表1 BFAST潜在火点统计数据

注：2004(270)表示2004年第270天；带*者为季节组分断点，其他为趋势组分断点

1.3相对火点的判别

将样本点NDVI和亮温的所有断点作为潜在火点，进行空间上下文分析，表2表明1号点2007年297天为林火火点，2号点没有林火火点。图6a和6b的火点内部分布图，可以发现1号点检测结果是正确的，2号点为火点像元，但本发明算法发生漏检，原因是2号点火点面积很小，且火点造成的NDVI和亮温时间序列数据趋势几乎没有变化，而BFAST算法通过结构性变化显著性检验判断时间序列数据中是否存在断点，在结构性变化不显著时会发生断点漏检的情况。

表2样本点的火点判别结果

2结果验证

根据样本点的火点判别结果：1号点2007年297天为林火火点。将2007年10月24日19：11-19：18 UTC Ikhana获取的圣迭戈AMS数据作为本发明的验证数据。

2.1AMS真实火点提取

利用AMS 11通道(3.7μm)和12通道(10.76μm)的亮温直方图进行火点提取。第一步，确定背景温度阈值。(图6b和6d中的黑色虚线Tb分别为两波段的背景温度阈值)；第二步，基于通道间比较方法(ICT)移除假火点(低于背景温度和10.76μm温度低3.7μm温度高的非火点。图6c中的白色三角点为ICT值，低于ICT的深色圆点亮温散点为非火点)；ICT值即为判断潜在火点与非火点的阈值。基本原理：将4μm通道观测值划分为10个柱状区间，计算每一个柱状区间11μm亮温的25％四分位数作为ICT值，如果11μm亮温的25％四分位数高于350K，ICT值设为350K。排除非火点的依据是：高于ICT值的点作为潜在火点，低于ICT值的点不是潜在火点；第三步，确定火点判定阈值，判别火点。(图6d灰色竖线Tf为10.76μm波段火点判定阈值；图6b灰色竖线实线Tf为3.7μm波段火点判定阈值。图6c黑色菱形散点为火点，图6a为AMS火点分布图，)。图7a和7b为AMS真实火点内部分布图。

2.2火点提取结果对比分析

将AMS火点提取结果作为火点真实数据，将本发明算法与MODIS V4算法火点提取结果进行对比分析。可以发现1号点检测结果是正确的，2号点为火点像元，但本发明算法发生漏检。

图8a-b底图为AMS数据中红外、近红外、绿光波段辐亮度对应RBG通道的合成影像，火点在中红外波段辐亮度高，在图像中显示为红色，植被在近红外高反射，在图像中呈绿色；黑框区域为火点提取算法验证区域，条纹像元为AMS提取的MODIS火点像元，作为MODIS的“真实”火点；图8a深色像元为MODIS V4算法的火点检测结果中判对的像元，置信度都为9，浅色像元为MODIS V4算法误判为火点的像元，条纹背景的像元为MODIS V4漏判的火点像元；图8b深色像元为本发明算法的火点提取结果中判对的像元，条纹背景的像元为本发明算法漏判的火点像元。

表3为MODIS火点统计信息，验证区共有46个火点像元，MODIS V4算法提取出44个火点像元，其中，判对火点41个火点像元，判错火点3个火点像元，另外，漏判5个；而本发明算法提取火点39个，判对火点39个，判错0个，漏判7个。表4为本发明算法漏检的火点像元统计信息，漏检火点对应的AMS火点像元数分别为4、2、58、27、17、2、7，火点面积都比较小。

表3火点验证结果

表4漏检火点的统计信息

分析可得：

(1)相对于MODIS V4火点提取算法，本发明提出的基于时间序列分析和空间分析的火点提取算法提高了火点判对率，这是因为本发明利用BFAST算法检测时间序列数据断点的方法判别潜在火点，而不是固定阈值确定潜在火点，因此本发明算法对于火点判别更具有自适应性和灵活性，克服了固定阈值判别潜在火点无法满足区域性差异的不足；

(2)本发明算法的漏检火点像元个数高于MODIS V4火点检测算法，这是因为火点面积很小，且火点造成的NDVI和亮温时间序列数据趋势几乎没有变化，而BFAST算法通过显著性检验判断时间序列数据中是否存在断点，在结构性变化不显著时发生断点漏检的情况。

3研究区火点检测

将本发明火点检测算法应用于整个研究区，得到2007年10月24日的火点检测结果图9a-b，底图为MOD02 1km真彩色合成图像，图9a为MODIS V4算法的火点检测结果，按照置信度从低到高分为7、8、9三个级别；图9b为本发明算法的火点检测结果，火点标记为红色深黑色，不包含置信度信息。从图中可以看出2007年10月24日圣迭戈火点分布比较分散，西北部、中北部和中南部都有火点存在。表5为两种算法的火点提取结果对比统计数据，可以得出本发明算法检测的火点总数(133)低于MODIS V4算法检测的火点总数(253)，这一方面由于MODIS V4火点检测精度依赖植被情况，森林砍伐活跃区误判率约为35％；另一方面，由于对于面积很小且火点造成的NDVI和亮温时间序列数据趋势几乎没有变化的火点，BFAST可能出现漏判的情况。

表5火点像元提取结果统计表

本发明针对MODIS数据的时间分辨率优势，提出时间序列分析和上下文空间分析相结合的MODIS火点检测方法；引入时间序列断点识别算法BFAST提取潜在火点，克服固定阈值区域性限制，从而对MODIS V4火点检测算法进行了改进；并将该算法应用于圣迭戈(San Diego)的火点检测，成功识别出研究区森林火灾情况。无论从方法论还是实证研究方面，对林火检测都具有重要意义。主要结论如下：

运用本发明算法，成功检测出1号样本点2007年10月24日的火点；同时识别出整个研究区当天的133个火点，这些火点分散于西北部、中北部和中南部。相对于MODIS V4火点检测算法，本发明的火点检测算法提高了火点检测的正确率，对于面积很小和未造成NDVI和亮温时间序列数据趋势明显变化的火点可能出现漏判情况。

空间上下文算法的原理是基于固定阈值提取潜在火点，再依据潜在火点与背景像元的空间关系进行相对火点判别，得到真实火点；本发明将时间序列算法引入来自动判别潜在火点，提高MODIS火点检测的灵活性，其中时间序列算法BFSAT能排除噪声干扰、提取断点无需设定阈值，相比固定阈值，可以适应区域性的差异，提高火点检测的精度。

除上述实施例外，本发明还可以有其他实施方式，凡采用等同替换或等效变换形成的技术方案，均落在本发明的保护范围内。

Claims (9)

1.一种基于MODIS的森林火灾火点检测方法，其特征在于，包括如下步骤：

(1)识别潜在火点：先对来自MODIS的若干天的MOD02及MOD03数据构建时间序列数据，得到该若干天的NDVI和亮温的时间序列日值数据，然后采用时间序列断点识别算法识别所述NDVI和亮温的时间序列日值数据中的断点，将识别出的断点作为潜在火点；

(2)判别相对火点：基于所述潜在火点及其背景像元，进行空间上下文分析，将识别出的火点作为相对火点，所述相对火点为真实火点。

2.如权利要求1所述的方法，其特征在于，所述步骤(1)中，在对所述MOD02及MOD03数据构建时间序列数据之前，先对所述MOD02及MOD03数据进行预处理。

3.如权利要求2所述的方法，其特征在于，所述预处理包括几何校正、数据有效性检查、亮温反演和NDVI反演。

4.如权利要求3所述的方法，其特征在于，所述几何校正包括去除“蝴蝶结”效应、投影转换和图像规则裁剪。

5.如权利要求3或4所述的方法，其特征在于，所述数据有效性检查包括对经几何校正后的数据进行重新赋值以及数据拼接处理。

6.如权利要求1所述的方法，其特征在于，步骤(1)中所述采用时间序列断点识别算法识别所述NDVI和亮温的时间序列日值数据中的断点的步骤如下：

(a1)利用可加性分解模型迭代拟合趋势组分和季节组分，表示为：

Y_t＝T_t+S_t+e_t (t＝1,...,n)；

其中：t为观测时间，Y_t为t时刻观测值；T_t为趋势组分；S_t为季节组分；e_t为噪声组分，是观测值中除去趋势组分和季节成分后的剩余部分；

设断点为则趋势组分T_t的线性模型表示为：

$T_t={\mathrm{\&alpha;}}_j+{\mathrm{\&beta;}}_{j}t,(t_{j-1}^{*}<t\&le;t_j^{*});$

其中：j为趋势断点所在位置，j＝1,...,m；α_j为截距，β_j为斜率，用来反演突变的幅度和断点前后渐变的斜率，对于断点t，若满足则断点t突变的幅度表示为：

Magnitude＝(α_j-1-α_j)+(β_j-1-β_j)t；

其中，β_j-1和β_j分别为断点前、后的渐变斜率；

设断点为则季节组分S_t的谐波模型表示为

$S_t=\underset{k=1}{\overset{K}{\&Sigma;}}\&lsqb;{\mathrm{\&gamma;}}_{j,k}sin\left(\frac{2\mathrm{\&pi;}kt}{f}\right)+{\mathrm{\&theta;}}_{j,k}\cos \left(\frac{2\mathrm{\&pi;}kt}{f}\right)\&rsqb;;$

其中，j为季节断点所在位置，j＝1,...,p；k为谐波项数，f为频率；γ_j,k＝a_j,k cos(δ_j,k)，θ_j,k＝a_j,k sin(δ_j,k)，a_j,k和δ_j,k分别为频率f/k的振幅和相位：和δ_j,k＝tan^-1(θ_j,k/γ_j,k)；

(a2)断点监测：

利用基于最小二乘法的移动求和检验是否存在断点，利用最小贝叶斯信息准则确定断点的最优数量，利用最小残差平方和估算断点的位置。

7.如权利要求6所述的方法，其特征在于，所述步骤(a2)的步骤如下：

(a21)如果最小二乘法的移动求和检验表明趋势组分T_t中有断点，运用最小二乘法拟合季节性因素调整后的数据得到趋势断点的数量和位置；

(a22)由基于M-估计的鲁棒回归估算得到趋势组分系数α_j和β_j，基于步骤(a1)中所述公式Magnitude＝(α_j-1-α_j)+(β_j-1-β_j)t，趋势组分估计值表示为：

(a23)如果最小二乘法的移动求和检验表明季节组分S_t中有断点，运用最小二乘法拟合趋势性因素调整后的数据得到季节断点的数量和位置；

(a24)由基于M-估计的鲁棒回归估算得到季节组分系数θ_j,k和γ_j,k，基于谐波模型，季节组分估计值表示为：

${\hat{S}}_t={\&Sigma;}_{k=1}^3{a}_{j,k}sin(\frac{2\mathrm{\&pi;}kt}{f}+{\mathrm{\&delta;}}_{j,k});$

(a25)以上步骤不断地迭代，直到断点的数量和位置不再变化。

8.如权利要求1所述的方法，其特征在于，步骤(2)中所述采用空间上下文分析法从所述潜在火点中判别出相对火点的步骤如下：

(b1)云体掩膜：

利用高级超高分辨率辐射计数据监测体，白天满足下式的像元为云体像元：

ρ₁+ρ₂＞0.9 or T₃₂＜265K or(ρ₁+ρ₂＞0.7 and T₃₂＜285K)；

其中，ρ₁、ρ₂、T₃₂分别为MODIS波段1的反射率、波段2的反射率、波段32的亮温；

(b2)确定有效背景像元：

以潜在火点为中心，向左、右、上、下分别延长1-5个像元，构建5×5、7×7……21×21大小的窗口，直到有效背景像元个数占窗口总像元数量的25％以上停止搜索，来确定窗口大小；

(b3)火点判别：

在获取所述有效背景像元之后，采用上下文空间统计分析方法对潜在火点进行判别，判别条件如下式，满足该条件则所述有效背景像元为火点像元，

$\begin{array}{c}(T_{22}-T_{31}>\stackrel{\&OverBar;}{T_{22}-T_{31}}+3.5{\mathrm{\&delta;}}_{T_{22}-t_{31}})and(T_{22}-T_{31}>\stackrel{\&OverBar;}{T_{22}-T_{31}}+6K)and(T_{22}>\stackrel{\&OverBar;}{T_{22}}+3{\mathrm{\&delta;}}_{T_{22}})\\ and\left\{(T_{31}>\stackrel{\&OverBar;}{T_{31}}+{\mathrm{\&delta;}}_{T_{31}}-4K)or({\mathrm{\&delta;}}_{T_{22}}^{\&prime;}>5K)\right\}\end{array};$

其中，为有效背景像元22和31波段亮温差的平均值，为有效背景像元22和31波段亮温差的平均绝对偏差；和分别为有效背景像元22和31波段的亮温差；和分别为有效背景像元22和31波段亮温的平均绝对偏差；为背景火点像元22波段亮温的平均绝对偏差。

9.如权利要求1所述的方法，其特征在于，所述步骤(2)后面还包括：

从自主模块化传感器数据中提取真实火点，将步骤(2)中所述相对火点与所述真实火点进行对比，验证所述相对火点的正确性和遗漏率。