一种森林火灾亚像元火点特性反演方法

技术领域

本发明涉及利用遥感信息监测火灾领域，尤其涉及一种森林火灾亚像元火点特性反演方法。

背景技术

森林火灾属于世界性的自然灾害，越来越受到各国政府的高度重视。频发的森林火灾不仅烧毁大量的森林资源，破坏自然环境和生态平衡，而且会严重影响人类的正常生活，甚至给人类带来灾难性的损失。卫星遥感具有监测范围广，时、空分辨率高等特点，在森林火灾监测中具有特殊的功能和应用前景。

MODIS是搭载在terra和aqua卫星上的一个重要的传感器，是卫星上唯一将实时观测数据通过x波段向全世界直接广播，并可以免费接收数据并无偿使用的星载仪器，全球许多国家和地区都在接收和使用MODIS数据。

MODIS传感器凭借其特殊的仪器特征参数广泛应用于森林火灾监测，但 MODIS火点仅仅代表含有火点的像元，火点像元实际为混合像元，包括背景组分、火点组分等。仅仅提取火点像元，并不能得到该像元实际火点面积、温度等特性，因此这样做存在很大的不确定性。目前对亚像元火点面积和温度的研究，主要有基于亚像元温度场理论的Dozier模型和基于多端元光谱混合分析法。 Dozier模型利用两个不同波段对亚像元高温物体的不同响应来检测亚像元火点，此方法假设一个火点像元由火点和背景两个组分构成，背景均质且温度相同。 Dozier、Prins、Giglio等分别将Dozier模型应用到AVHRR数据、GOESVAS数据以及MODIS数据中；Dozier亚像元求解使用牛顿迭代法，求解时需估算背景温度，背景温度的估算精度直接决定了亚像元火点特性的计算精度，张顺谦等利用加速遗传算法对林火亚像元辐射率平衡方程求解，结果表明加速遗传算法比牛顿迭代法具有更高的计算精度和可信度；Dozier算法不同波段的选择对结果精度影像较大，Giglio利用MODIS不同波段组合得到的亚像元火点面积与温度差异很大，刘诚利用AVHRR数据不同通道组合，提出了亚像元火点特性估算的3种方法，精度可达10-4以上；Dozier算法对波段见配准精度很敏感，Shephard等证明，MODIS 1km传感器10％的配准误差会导致火点温度和面积误差分别为150K、210％。Matson、Prins、Giglio、Peterson等对Dozier模型进行了改进。为解决混合像元端元在光谱和空间上的变化问题，Roberts等首次提出了多端元光谱混合分析(MESMA)方法，即通过将像元光谱信息和空间信息相结合，逐像元调整端元的数量与类型动态进行混合像元分解，此方法广泛应用在土地覆盖/利用分类、监测城市环境等方面。Dennison等首次将MESMA用于林火火点特性研究。通过建立发射辐射端元光谱库(500-1500K)和反射辐射端元光谱库，拟合并得到亚像元火点面积和温度，但是地表端元的光谱相似性和烟雾的存在影响模型精度。Eckmann等将此方法用于MODIS数据亚像元火点特性研究，并验证了MESMA对于研究亚像元火点面积和温度的实用性和优越性。Eckmann等用 MESMA对ASTER数据亚像元火点温度和面积进行研究，并进一步探讨了火点的季节变化。Dennison研究了不同波段选择对MESMA反演结果的影响。

针对Dozier方法假设背景热特性相同，且只将像元分为火点和背景两种组分、验证困难且存在很多误差源的不足，MESMA方法考虑到端元数目和光谱的变化，为亚像元火点特性研究提供了新的方法和可能性，具有更好的应用前景，但是MESMA建立的发射辐射端元光谱库仅有植被背景光谱信息，不包含背景为烟羽的光谱信息，使得MESMA无法适用于所有火点像元的亚像元火点特性研究。

发明内容

为了克服现有的缺陷，本发明提出一种森林火灾亚像元火点特性反演方法。

该方法包括以下步骤：步骤A：接收设定区域的MODIS数据并探测火点；步骤B：建立火点发射辐射端元库；步骤C：建立反射太阳辐射端元库；步骤D: 建立阴影端元库；步骤E：基于MESMA的亚像元火点特性反演，计算出亚像元火点面积、温度、背景类型信息。

进一步的，步骤B为：步骤B01：结合步骤A得到数据，设置模型参数，包括大气模型、光谱反射率、消光模型、传感器高度、二氧化碳柱、地面高程、传感器天顶角、可降水量；步骤B02：设置模拟的温度范围和温度间隔，将模型参数输入MODTRAN软件得到参数文件；步骤B03：根据参数文件运行MODTRAN软件得到光谱响应文件，包括不同温度下每个波长的辐亮度值；步骤B04：将辐亮度值分别重采样到MODIS每个响应波段，得到MODTRAN模拟的各波长的MODIS 波段响应率。

优选的，在步骤B02中，模拟的温度范围是500K-1500K，温度间隔为10K。

进一步的，步骤C包括：步骤C01：将端元分为灰、丛林、农田、草地、橡树林、岩土和烟羽7种组分类型；步骤C02：在google earth中勾画每类地物类别的矢量边界，为每类端元生成kml文件；步骤C03：kml文件转成shp文件，定义投影，并在遥感处理软件中提取反射太阳辐射端元光谱，作为初始光谱；步骤C04：选取最优单元。

优选的，步骤C02中，边界为多边形，边长大于2km。

进一步的，步骤C04包括：步骤C0401：建立方阵，记录步骤C03中初始光谱模拟反射辐射端元中不同背景像元光谱所产生的均方根误差、光谱角、端元丰度和约束限制编码结果；其中，约束值编码有3种取值：0、1、2，满足丰度和均方根误差阈值时，约束值编码取0；部分约束满足均方根误差约束但超越丰度阈值，约束值编码取1；完全约束超过均方根误差约束时，约束值编码取2；步骤C0402：根据端元均方根误差法、最小平均光谱角法和计数端元选择法选取最优端元。

其中，部分约束模型具有最大丰度、最小丰度及均方根误差阈值的约束；当丰度值越界时，将根据阈值内的丰度值计算均方根误差，如果得到的均方根误差在其阈值内，则表示该模型满足部分约束；如果得到的均方根误差不在其阈值内，则表示该模型不满足部分约束。

进一步的，步骤D为：根据步骤B建立阴影端元库，并且在步骤B02中，温度为10K。

进一步的，步骤E包括：步骤E01：对每个像元进行端元组合；步骤E02：计算每个像元的均方根误差，按所有波段B的均方根误差最小的原则选取每个像元的模型；步骤E03：计算得到亚像元的面积、温度、背景类型。

优选的，步骤E02中，波段为21、29、31、32。

本申请基于能量平衡方程，将火点像元传感器接收的辐亮度分为火点发射辐射、地表反射太阳辐射和阴影辐射三部分，其中，火点发射辐射加入烟羽背景，并基于MESMA模型反演得到亚像元火点的面积、温度和背景信息。通过与 AMS火点结果、MOD14热异常产品FRP进行相关性分析后，结论如下：相比FRP 与AMS火点的相关性(r＝0.293)，MESMA火点面积与AMS火点结果具有更高的相关性(r＝0.743)，MESMA火点温度与AMS火点同样具有显著的相关性(r＝0.339)；且加入烟羽的MESMA模型丰富了火点背景信息，使MESMA方法适用于植被、烟羽背景的亚像元火点特性研究，扩大了MESMA的使用范围。

附图说明

图1是根据本发明一个实施例的森林火灾亚像元火点特性反演方法的流程图；

图2是根据本发明一个实施例的MODTRAN模拟的火点端元传感器辐亮度示意图；

图3是根据本发明一个实施例的背景端元辐亮度的均值和标准差示意图；

图4是根据本发明一个实施例的阴影端元辐亮度模拟结果示意图；

图5是根据本发明一个实施例的样本点分布示意图；

图6是根据本发明一个实施例的美国加利福尼亚州概况图；

图7是根据本发明一个实施例的利用MODIS V4算法的火点提取结果；

图8是根据本发明一个实施例的AMS的火点提取流程示意图；

图9是根据本发明一个实施例的样本点的真实火点分布图；

图10是根据本发明一个实施例的FRP和MESMA火点面积关系图；

图11是根据本发明一个实施例的FRP和AMS火点像元数关系图；

图12是根据本发明一个实施例的AMS火点像元数和MESMA火点面积关系图。

具体实施方式

下面结合附图和具体实施例对本发明提供的一种森林火灾亚像元火点特性反演方法进行详细描述。

本发明提供一种森林火灾亚像元火点特性反演方法。图1所示为本发明的具体步骤，包括：步骤A：接收设定区域的MODIS数据并探测火点；步骤B：建立火点发射辐射端元库；步骤C：建立反射太阳辐射端元库；步骤D:建立阴影端元库；步骤E：基于MESMA的亚像元火点丰度反演，计算出亚像元火点面积、温度、背景类型等信息。

1.端元光谱库建立

端元光谱是火点像元分解的参考波谱，其质量和精度直接影响火点像元混合分解精度。本发明端元光谱库分为火点发射辐射端元库、反射太阳辐射端元库和阴影端元库。

1.1火点发射辐射端元库的建立

MODTRAN(中等光谱分辨率大气透过率及辐射传输算法软件)是目前对遥感影像进行大气纠正的主要大气辐射传输模型之一。针对MODIS数据特点，本发明利用MODTRAN4辐射传输模型建立发射辐射端元库模拟传感器端辐射值。MODTRAN将火点视为黑体，建立500-1500K(227-1227℃，间隔 10K)共101个发射辐射端元的辐亮度数据。为MODTRAN建立发射辐射端元的重要参数设置，是研究区各像元参数的中位数，所有输入数据都可以从 MODIS同步产品中得到。

表1发射辐射模型的MODTRAN参数

MODTRAN模拟发射辐射的参数文件：MODTRAN tape5文件用于存放所有参数，包括控制运行参数(CARD 1、CARD 5)、传感器参数(CARD 1A3、CARD 4)、大气参数(CARD 1、CARD2)、观测几何参数(CARD 1、CARD 3)、地表参数(CARD 1、CARD 4)和输出参数六类。

本发明使用的MODTRAN版本为Mod4v1r1(MODTRAN4Version 1 Revision 1)，需要模拟的温度间于[500-1500K]，温度间隔10K，共101 组参数，将所有的参数写入tape5文件，运行软件。本发明利用输出的.plt 文件，得到101组数据，每组数据由N行组成，每行一个波长对应一个辐亮度；将MODIS光谱响应文件内插分别得到MODTRAN模拟的各波长的MODIS 波段响应率，根据公式(1)将辐亮度值分别重采样到MODIS每个响应波段。图2中(a)(b)(c)分别为采样前500K、1000K和1500K在各波长的辐亮度； (d)(e)(f)分别为采样后500-900K、1000-1500K、500-1500K(间隔100K) 在MODIS各波段的辐亮度。

其中，L(λ_i)和f(λ_i)分别为波长λ_i对应的辐亮度值和光谱响应率，每隔Δλ就计算一个辐亮度值，Δλ为2nm。

1.2反射太阳辐射端元库的建立

1.2.1光谱选取

多端元光谱混合分析中的分解端元获取途径有两种：直接从已知的端元库获取“参考端元”；直接从遥感影像中提取“图像端元”。后者的端元光谱与需要分解的遥感影像具有相同的测量比例和相同状态的噪声，且容易获取，但对影像空间分辨率要求较高。本发明需要获取反射太阳辐射数据，因此选择直接从遥感影像获取不同地物的辐亮度数据的方法。首先根据研究区遥感影像，通过图像特征空间几何顶点将研究区组分类型，将端元分为灰(ash)、丛林(chaparral)、农田(farmland)、草地(grass)、橡树林(oak forest)、岩土(soilrock)、烟羽(smoke)等7类；在google earth中勾画每类地物类别的矢量边界(为保证每个多边形至少覆盖一个 MODIS像元，要求所有的多边形边长大于2km)，将生成的7个kml文件转成shp文件，定义投影，并在ENVI(The Environment for Visualizing Images，遥感处理软件)中提取反射太阳辐射端元光谱，作为初始光谱。 ENVI遥感处理软件具有输入、输出图像、图像定标、增强、校正、融合以及各种变换、信息提取、图像分类等功能。

1.2.2最优端元的选取

步骤一：建立方阵

方阵记录所有光谱组合模拟反射辐射端元中不同背景像元光谱所产生的均方根误差、光谱角、端元丰度和约束限制编码结果。部分约束模型具有最大丰度、最小丰度及均方根误差(RMSE)阈值的约束，但是当丰度值越界时，将根据阈值内的丰度值计算均方根误差，如果得到的均方根误差在其阈值内，则表示该模型满足部分约束；如果得到的均方根误差不在其阈值内，则表示该模型不满足部分约束。

(1)均方根误差RMSE

RMSE矩阵第i行第j列元素RMSE_i,j表示光谱j模拟光谱i的均方根误差，>

其中，ε_λ为光谱j模拟光谱i各波段亮度差。

(2)光谱角SA

SA矩阵第i行第j列元素SA_i,j表示光谱j模拟光谱i的角距，表达式为：

SA矩阵为对角线对称矩阵，即光谱A模拟光谱B的光谱角等同于光谱 B模拟光谱A的光谱角。

(3)端元丰度EM fraction

端元是非阴影端元，EM Fraction矩阵的i行j列为利用端元i光谱混合分析端元j得到的丰度，由于亮度差异导致SMA丰度和RMSE模型的差异， EM Fraction矩阵不是对角线对称矩阵。当亮光谱i模拟暗光谱j时，SMA 丰度间于[0，100％]且RMSE为暗光谱j中低于光谱i的亮度差异；而当暗光谱j模拟两光谱i时，SMA丰度将大于100％且RMSE为亮光谱i中高于暗光谱j的亮度差异，因此后者的RMSE会比前者的RMSE大。

(4)约束值编码Constraint Code

约束值编码有3种取值：0、1、2。满足丰度和RMSE阈值时，约束值编码取0；部分约束满足RMSE约束但超越丰度阈值，约束值编码取1；完全约束超过RMSE约束时，约束值编码取2。

步骤二：选取最优端元

在实际多端元光谱混合分析中，每个像元仅用一个最优模型来模拟，因此需要用所有类别的所有端元逐个模拟影像中的像元，确定该端元对所在类别库中端元光谱的代表性。选取最优端元的指标有：

(1)端元均方根误差法EAR

EAR(Endmember Average RMSE)作为类内最优端元光谱定位的方法，通过计算端元库中每个光谱模拟同类别其他所有光谱的RMSE均值。每类中 RMSE均值越小，代表此端元对于该类别所有端元的信息表达度越高，越能作为最佳端元。

其中，EAR_i为光谱i的端元均方根误差，N为该类别中光谱数量，RMSE_i,j为光谱i模拟光谱j的均方根误差。

(2)最小平均光谱角法(MASA)

MASA(Minimum Average Spectral Angle)类似于EAR，但MASA衡量标准是光谱角而非RMSE。通过计算每个光谱与端元库同类别其他所有光谱的光谱角均值，最优MASA候选端元为平均光谱角最小的光谱。

其中，MASA_i为光谱i模拟类内其他光谱的光谱角平均值，SA_i,j为光谱i模拟类内光谱j的光谱角。

(3)计数端元选择法(CoB)

CoB(Count based Endmember selection)是为模拟类别中最大数量光谱的光谱库选取最优端元的方法。即通过候选光谱组合在分解类别中其他光谱时是否满足丰度、RMSE、残差约束三种指标来评估次该光谱组合的优度，这些信息存放在约束值编码矩阵。计算得到每个光谱组合被模拟为类内光谱的次数(in_CoB)、被模拟为内外光谱的次数(out_CoB)，类别最优端元应该是in_CoB最高和out_CoB最低的端元，代表该端元对类内光谱信息表达度高，且有较低的类间混淆度，最理想的端元是只能模拟类内光谱，不能模拟类外光谱。

前两种方法是为了选择反射辐射端元库中最优的模拟同类别光谱的端元，后一种方法是综合RMSE、丰度等指标选择的类内表达度高、类间混淆度低的端元。通过上述方法，本发明建立了13种反射辐射端元光谱，图3 为背景端元辐亮度均值和标准差图。

1.3阴影端元光谱的建立

考虑到因地形和气象因素生成的阴影图像中，阴影会影响像元整体灰度，因此在混合像元分解时必须要考虑阴影端元的提取。本发明阴影端元光谱的建立同发射辐射端元光谱的建立方法一样，但MODTRAN模拟的边界温度为10K，将生成的pltout文件的辐亮度借助MODIS光谱响应率文件转换为MODIS各个波段的辐亮度数据，如图4所示。

2.基于MESMA的亚像元火点特性反演

MESMA混合像元分解原理是将混合像元辐亮度看作各端元辐亮度的线性组合，允许每个像元有不同的端元组合，表示为：

其中，N＝3，L_iλ表示位置i的混合像元，即3(N＝3)个端元L_kλ的总和，每个端元的端元丰度为f_ki，剩余部分ε_iλ表示不在端元库的端元，每个像元的模型选择标准是对于所有波段B的RMSE最小(式(7))。

考虑水汽影响，本发明基于MODIS的21、29、31、32波段辐亮度进行多端元光谱混合分析，因此B＝4，分解得到亚像元火点面积、温度、背景类型等信息。

实施例1：

(一)研究区及数据处理

1.研究区概况

圣迭戈(San Diego)是美国加利福尼亚州的一个太平洋沿岸城市，位于加利福尼亚州的西南端，南距墨西哥边境20公里，西经117.6°- 116.0°，北纬32.5°-33.6°。人口134万，面积963km2，是加利福尼亚州的第二大城市。迭戈市年平均气温约为华氏70度，以温暖的气候和多处的沙滩著名，属温带大陆性干旱气候，阳光充足，降雨较少，年均降水量约为25cm。图6为加利福尼亚州概况，从(b)-(d)中我们可以看出圣迭戈发生破坏严重、严重伤亡、大规模火灾的次数分别为4、5、3，属于火灾较严重的地区。参考MCD12Q1的第一种分类体系的分类结果：1常绿针叶林，2常绿阔叶林，3落叶针叶林，4落叶阔叶林，5混交林，6封闭灌丛，7开放灌丛，8多树草原，共8类像元作为本发明的林区范围，本发明选取林区内的两个样本点：1号点、2号点作为研究对象，如图5所示。

2.实验数据

本实施例使用的数据包括实验数据和验证数据。

(1)实验数据为样本点2007年10月24日的火点检测结果，如图7 所示，火点提取算法为MODIS V4火点检测方法。经验证，以该算法提取的 41个正确火点作为研究火点。

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

(二)结果验证与分析

1.AMS真实火点提取

利用AMS 11通道(3.7μm)和12通道(10.76μm)的亮温直方图进行火点提取。第一步，确定背景温度阈值。(图8(b)(d)中的黑色虚线Tb 分别为两波段的背景温度阈值)；第二步，基于通道间比较方法(ICT)移除假火点(低于背景温度和10.76μm温度低3.7μm温度高的非火点。图 8(c)中的白色三角点为ICT值，低于ICT的黑色圆点亮温散点为非火点)；第三步，确定火点判定阈值，判别火点。(图8(d)黑色实线Tf为10.76μm 波段火点判定阈值；图8(b)黑色实线Tf为3.7μm波段火点判定阈值。图8(c)黑色菱形点为火点，(a)为AMS火点分布图)。图9为AMS真实火点内部分布图。

2.亚像元火点特性验证与分析

将MODIS的41个火点进行混合像元分解，来研究亚像元火点特性。可以得出(详细结果见表2)，大部分火点面积小于1％，大部分为背景组分，充分说明了亚像元火点特性研究的意义。

表2 2007年10月24日火点详细数据

41个火点的小样本量限制了本发明结论的显著性，但通过简单相关性分析和线性回归分析，我们可以得到一些结论。本发明将MOD14热异常产品中的火点辐射功率FRP与MESMA反演的火点面积(简称MESMA火点面积)、 AMS火点数进行相关性分析和线性回归分析；本发明将MESMA火点温度和 AMS火点数进行相关性分析来验证MESMA反演的火点温度(简称MESMA火点温度)的精度。

2.1 MESMA火点面积分析

MESMA火点面积与FRP、AMS火点数量三者的相关性统计分析结果见表 2，具体分析结果如下：

表2 FRP、MESMA火点面积、AMS火点数量相关性统计分析

注：^**在0.01水平(双侧)上显著相关，r表示Pearson相关系数，p表示相关系数显著性检验

FRP与MESMA火点面积Pearson相关系数双侧显著性检验p＝0.541，二者不具有相关性(如图10所示)。可知，FRP和MESMA所描述的火点信息完全不同。例如，FRP分别为152.5MW和480.3MW的两个火点，MESMA火点面积分别为3.66％和1.03％，温度分别为620K和500K，这种矛盾可能由于背景辐射不同造成。MESMA反演的背景类型为植被和烟羽时，FRP和MESMA 火点面积同样不具有相关性。

FRP和AMS火点像元数的Pearson相关系数为0.293，中等相关，相关系数双侧显著性检验p值为0.063(如图7所示)；背景为植被时二者的 Pearson相关系数r为0.381，相关系数双侧显著性检验p值分别为0.080，具有一定的统计学意义；烟羽的加入使得二者相关性下降了23％。总体来说，虽然FRP能反映AMS火点数量的变化，但是相关性很低，造成FRP和AMS火点像元数相关性很低的另一原因是：AMS火点像元为混合像元，并不是实测火点面积。

AMS火点像元数与MESMA火点面积的Pearson相关系数r为0.743^**，在0.01水平(双侧)上显著相关，具有统计学意义(如图8所示)。烟羽端元的加入对二者极其显著的相关关系(仅有植被p<0.01)没有改变。总体来说，MESMA火点面积与AMS火点像元数相关性(r＝0.743，p<0.01)高于FRP与AMS火点像元数相关性(r＝0.293)。

2.2 MESMA火点温度比较分析

MESMA火点温度绝大多数大于1000K，只有5个为500K，尚没有数据可以验证本发明火点温度的反演精度。本发明仅将MESMA火点温度与AMS 火点像元数的相关性分析统计，结果见表3，可知MESMA火点温度与AMS 火点像元数的Pearson相关系数r为-0.339*，在0.05水平(双侧)上显著相关，当背景为植被时，MESMA火点温度与AMS火点像元数的Pearson 相关系数r为-0.707**，在0.01水平(双侧)上显著相关。总之，烟羽端元的加入使MESMA火点温度与AMS火点像元数极其显著的相关性降低，但二者相关性仍是显著的。

表3 MESMA火点温度相关性分析统计结果

注：*在0.05水平(双侧)上显著相关，**在0.01水平(双侧)上显著相关；r表示Pearson 相关系数，p表示相关系数显著性检验。

结论与讨论

本发明针对像元级MODIS火点面积、燃烧组分及温度等特性的不确定性，开展基于多端元光谱混合分析的MODIS亚像元火点特性研究：(1)利用辐射传输模型和目标MODIS影像的大气参数，建立不同温度火点发射辐射端元光谱库；(2)从图像中选取烟羽区和非烟羽区(植被区)的背景端元，分别建立背景反射辐射光谱库；(3)建立阴影光谱库；通过对MODIS火点进行MESMA分析，确定火点亚像元中特性，得到亚像元火点面积、温度以及背景信息。亚像元火点特性结果表明：大部分火点面积小于混合像元面积的1％，说明亚像元火点的研究意义。本发明利用AMS火点像元数以及FRP验证MESMA反演的火点面积与温度。结果表明MESMA反演的火点面积与AMS火点像元数相关性(r＝0.743，p<0.01)高于FRP与AMS火点像元数相关性(r＝0.293)，烟羽端元的加入对二者极其显著的相关关系(仅有植被p<0.01)没有改变；烟羽端元的加入使MESMA反演的火点温度与AMS 火点像元数极其显著的相关性(r＝-0.707)降低，但二者相关性仍是显著的(r＝-0.339)。总体而言，烟羽信息的加入扩大了MESMA使用范围，使得火点亚像元分解结果信息更丰富。