基于模糊集理论的单时相跨区域森林过火区遥感制图

朱 曦，覃先林，廖 靖

(中国林业科学研究院资源信息研究所，北京 100091)

0 引言

森林过火区监测对灾后重建和温室气体定量评估都非常重要。卫星遥感技术已成为森林过火区识别的一种重要技术手段。环境减灾一号小卫星B星(HJ－1B)中分辨率多光谱CCD相机数据具有3个可见光波段和1个近红外波段，正常植被在近红外波段反射率较高，在可见光部分主要吸收红光[1];而在燃烧过后，由于叶片组织的破坏，植被被裸露的木炭和土壤代替，近红外波段的反射率会降低，而可见光波段反射率会上升[2]。因此，在多时相的森林过火区识别中，近红外和可见光是非常有效的波段。然而，多时相卫星数据获取难度大，物候、时相上的差异需要更精确的大气校正和几何配准，而且计算量较单时相也会更大[3]。

多种研究方法曾广泛应用于森林过火区遥感制图中。Henry[4]将佛罗里达作为试验区进行过火区制图研究，并对最大似然分类法与分类回归树法进行了对比，结果最大似然分类法的精度比分类回归树法更高;Cassidy[5]在湿地的过火区制图研究中将ISODATA的非监督分类法加入了整个分类体系，为后期的过火区提取提供了基础;Petropoulos等[3]则对基于人工神经网络(artificial neural network，ANN)方法的过火区制图进行了评估，该方法总体精度高达90.29%，其研究结果表明ANN在地中海地区过火区识别中具有很好的潜力。上述方法中都利用了各种基于不同波段的光谱指数[6]，然而光谱指数的优选以及阈值的确定仍然存在一定的困难;至今，仍没有被一致认为是过火区识别最优的一组或一个光谱指数，而阈值的设定又受到很多因素(火前植被状况、火烧严重程度等)的影响。因此，基于硬性阈值的识别方法不够稳定强健，很可能在不同区域的精度差异非常大[7]。

针对上述问题，本文主要围绕2个方面进行研究:①尝试利用单时相遥感数据进行过火区识别。单时相较多时相而言，遥感图像中的水体、阴影等与过火区易混淆的区域难以区分;为解决这个问题，先将HJ－1B IRS(红外多光谱相机)数据重采样成30 m空间分辨率，然后对过火区进行识别。HJ－1B IRS数据具有短波红外和热红外波段，短波红外对植被水分含量非常敏感，由于火后植被水分大量减少，木炭和土壤的部分裸露会导致短波红外反射率上升[8]，因此短波红外波段数据对水体和阴影等与过火区在可见光近红外波段易混淆的区域具有一定的识别能力[9]。另外，一些学者发现过火区的地表温度(land surface temperature，LST)会有所上升，故本文也加入了LST参数，对过火区进行评估。②采用模糊分类法进行森林过火区制图。其优点在于使用所有可用的光谱指数，自适应地突出不同指数的优势和抑制冗余信息、增强过火区信息，而不需要对光谱指数进行优选;在跨区域遥感制图中，该方法会更加强健。此外，利用软分类法，不硬性设定阈值，使用多指数隶属度函数生成“正面信息”和“负面信息”来减少误判，提高了识别精度[10]。

1 实验区与数据

1.1 实验区

本文选择的实验区位于中国西南、华北、东北和俄罗斯等曾经发生过森林火灾的地区，这些地区都是森林火灾的频发区域。其中，山西阳泉、内蒙古根河和云南安宁为训练样本区域，俄罗斯斯科沃罗季诺和黑龙江逊克为验证样本区域。

1.2 数据获取与预处理

本文所使用的数据为我国环境减灾一号卫星B星(HJ－1B)数据[11](表1)。

表1 环境减灾一号卫星B星参数Tab.1 Parameters of HJ－1B satellite

本文实验区中各实验地点的遥感数据获取情况见表2。

表2 研究区遥感数据Tab.2 Remote sensing data of research area

本文主要采用了HJ－1B CCD的4个波段和HJ－1B IRS的短波红外(1.55 ～1.75 μm)及热红外(10.5 ～12.5 μm)波段数据。IRS 数据的空间分辨率较低，而CCD数据的空间分辨率较高，因此对这2种载荷数据进行了融合处理，以达到图像空间分辨率较高而利于森林过火区识别的目的。

对HJ－1B CCD和IRS数据分别按照其图像头文件提供的定标系数进行辐射定标。定标后对可见光、近红外和短波红外波段数据进行了大气校正，以消除大气和光照等因素对地物反射的影响，获得地物的反射率。大气校正使用 ENVI 4.8提供的FLAASH大气校正模块，其中波谱响应函数由国家资源卫星中心提供[11]。

同一天的CCD和IRS数据由于成像条件一致，可以较好地完成几何配准。以CCD数据为参考图像，IRS数据为待配准图像;由于CCD数据的第4波段和IRS的第1波段均为近红外波段，图像的地物波谱响应较为一致，图像具有较强的相似性，因此采用这2个波段自动计算匹配点[12]，匹配误差小于20 m;然后对IRS数据进行重采样，将其空间分辨率重采样到30 m。

2 研究方法

本文研究的方法是基于一个区域生长的过程，该过程需要一个初始“种子”像元和一个限制生长准则。总体技术流程如图1所示。

图1 总体技术流程Fig.1 Flowchart of general technology

本文通过定义一系列隶属度函数，对所有使用的光谱指数SIi和LST进行归一化处理，生成突出过火区的“正面信息”(PEi)和“负面信息”(NEi)，每个像元的PEi和NEi都是在[0，1]范围内。对于PEi，像元值越接近1表示它是过火区的可能性越大;而对于NEi，像元值越接近1表示它是非过火区的可能性越大。利用模糊集理论中不同的算子对这些PEi和NEi进行聚合，生成聚合后的正面信息PEi′(即候选“种子”和候选边界);然后利用聚合后的负面信息NEi′对候选“种子”和候选边界进行修正，提高“种子”像元选取和限制生长准则边界的精度;最后利用最终修正过的“种子”像元和生长边界进行区域生长，得到森林过火区的提取结果。

2.1 地表温度反演

本文利用辐射传输方程法(即大气校正法)对热红外波段数据进行LST反演。其基本思路为:首先，利用与卫星过顶时间同步的实测大气数据来估计大气对地表辐射的影响;其次，把这部分大气影响从卫星高度上传感器所观测到的热辐射总量中减去，得到地表热辐射强度;再通过普朗克方程将其转化为LST(其具体计算方法见文献[13]);最后，对得到LST进行归一化处理。所使用的辐射传输方程为式中:Lλ为卫星高度上传感器测得的辐射强度;τ为大气在热红外波段的透过率;ε为地表辐射率;TS为地表真实温度;B(TS)为用普朗克定律推导得到的黑体在TS中的热辐射亮度;L↑和L↓分别为大气上行和下行辐射。

2.2 光谱指数计算

模糊分类法不需要对光谱指数进行优选，可以通过算法本身抑制冗余信息和突出有用信息。因此，本文使用了在森林过火区遥感制图中常用的一系列光谱指数，其计算公式和参考文献见表3。

表3 本文使用的光谱指数Tab.3 Spectral indices used in this paper

2.3 隶属度函数提取过火区信息

采用隶属度函数[19]对过火区的正面信息PEi及负面信息NEi进行有效提取。隶属度函数的定义为:设Z为一个集合，Z={z};给定一个映射μA:Z→[0，1]，使得 Z中的每一个元素 z都有一个A(z)∈[0，1]与之对应，则确定了Z中的一个模糊集合A，μA称为模糊集合A的隶属度函数，μA在z∈Z点处的值μA(z)称为z对A的隶属度。μA的值越接近1，则z对A的隶属度越高。即

通过定义一系列隶属度函数，将光谱参数转化为值为[0，1]的PEi和NEi图层，用于增强过火区的信息或非过火区的信息。PEi的定义采用了NDVI，EVI，SAVI和CSI;NEi的定义采用了NBR和LST。定义隶属度函数有多种方法，本文采用数据统计驱动的方式，对每个类别选取大量样本，根据类别的直方图进行隶属度函数的定义。训练样本的选择是综合4个训练区域完成的，而类别包括了过火区、阴影、水体和植被。本文选取了不同区域和不同燃烧程度的过火区训练样本，以增强隶属度函数的强健性。

2.4 模糊聚类

通过一些模糊聚类算法分别将一系列PEi和NEi图层进行聚类，形成所需要的候选“种子”像元、候选边界和负面修正图层。

1)对候选“种子”和候选边界的聚合。采用基于模糊的有序加权平均(ordered weighting averaging，OWA)算子[20]，其定义为

设 OWA:Rn→R，即

式中:w=(w1，w2，…，wn)为与函数 OWA 相关联的加权向量为(a1，a2，…，an)中第 j个大的元素，即对(a1，a2，…，an)从大到小进行了排序，而ai与wi没有任何关联(wi只与集合中第i个位置有关)。

有序加权向量w=(w1，w2，…，wn)可确定为

式中:μQ为模糊语义量化算子;Q为通过非递减函数μQ得出的一个语义量化模糊集，可用来对聚合策略进行一个语义上的定义。

图2 “most”语义量化函数Fig.2 Linguistic relative quantifier function“most”

图2 给出了“most”这个语义的量化函数[21]。这些函数从“most50%”到“most90%”，表示对Q中有效元素的满意程度;这些函数都是分段函数，满足3 个基点:(0，0)，(x%，0)和(1，1)。x%表示提供有效正面信息的最低百分比，x值越大，聚合的语义越严格。为了得到满意结果，需要有更多的图层都对过火区具有正面增强。

对“种子”像元选取的要求是尽量减少误判。因此，本文利用了非常严格的算子“most90%”，要求一个像元的90%信息都是增强过火区的正面信息时，才认为这个像元是过火区;而对于限制生长准则的边界，则选用了较为宽松的聚合算子“most50%”，即只要一个像元的50%信息是增强过火区的正面信息，就将这个像元归并。

2)对负面信息的聚合。为了突出非过火区像元(主要是水体和阴影)，需对正面信息进行修正从而减少误判。本文采用了另一种模糊算子“Max”[19]，因为只有最大化地突出非过火像元、进而生成负面信息，通过负面信息修正剔除水体和阴影等易与过火像元混淆的类别，才能有效地减少误判。“Max”算子加权向量W的定义为

比较式(3)和式(5)可以看出:

式中NE(negative evidence)为负面信息。

继而可通过下式计算出修正后的“种子”像元和限制生长准则边界图层，即

式中:cPESEED和cPEGROW分别为候选“种子”和候选边界;rPESEED和rPEGROW分别为修正后的“种子”和限制生长准则边界图层。

通过设定阈值在rPESEED图层中选取过火“种子”像元，在rPEGROW图层中进行区域生长。生长准则为只要rPEGROW像元的值大于0，则将其归并入过火像元，直到没有邻域像元可以加入。

3 结果和讨论

3.1 隶属度函数定义

本文中隶属度函数的定义是通过不同类别直方图获取的，图3给出了不同类别的归一化燃烧率(normalized burn ratio，NBR)和归一化地表温度(normalized LST，nLST)的直方图。

图3 本文所选类别直方图Fig.3 Histograms of different classes

从图3可以看出，过火区和水体阴影在基于可见光和近红外的NDVI(图3(a))和EVI(图3(b))的直方图中难以区分，而NBR(图3(c))和 nLST(图3(d))对于水体阴影等都具有很好的可分性。基于可见光和近红外的光谱指数对过火区和植被具有一定的可分性[2]，但水体和阴影等类别使用这些波段的数据就难以区分。因此，在对负面信息修正时，本文采用了NBR和nLST这2个参数。NBR应用于过火区识别具有很长的历史，植被燃烧后叶片的损毁导致水分减少，使得过火区的短波红外反射率比水体和阴影高[5];而近红外波段的反射率相近，因而过火区的NBR值会小于水体和阴影。很多学者也发现火后过火区温度会有所上升[6]，本文加入nLST对正面信息进行修正，由图3可以看出，过火区的nLST比除裸地外的其他类别偏高。

3.2 “种子”像元阈值和语义量化函数选择

对训练样本区域用不同语义量化函数(包括未经过负面信息修正的和经过负面信息修正的语义量化函数)和不同阈值提取的“种子”像元进行了误判评价，选择的阈值包括[0，1]区间内每间隔0.1的值。通过目视解译，在原图像中选取过火像元和非过火像元样本作为分类参考数据，与提取的“种子”像元对比，得到误差矩阵，从而找出误判像元。图4给出了3个训练区的误判评价。

图4 “种子”像元误判Fig.4 Commission error for“seed”pixels

由图4可以看出，3个训练区域在经过负面信息修正前(空心圆点)，“种子”像元的整体误判率较高;随着阈值的增大，误判率显著降低;阈值不变时，随着语义量化函数从“most60%”到“most90%”，算子越来越严格，误判率也显著降低。

定量地看，未修正的算子无法满足“种子”像元误判率要很低的要求;经过修正的算子误判率与修正前的算子相比显著降低，而阈值对经过修正的算子影响很小(阈值即使在最宽松的0.1处，误判率也小于10%)。对修正后的“most90%”而言，在阈值为0.5处的误判率已经在1%以下，因此选择语义量化函数“most90%”和阈值0.5作为选取“种子”像元的准则。修正的算子是NBR和LST隶属度函数的聚合，这也证明了短波红外和温度信息对过火区提取的重要性(这2个波段大大减少了对于水体和阴影的误判率，提高了过火区提取的精度)。

对于“种子”像元，采用的是“most90%”这个函数和0.5的阈值;而区域生长的生长准则应相对宽松，才能形成过火区的边界，减少漏判。本文的生长准则是利用量化函数“most50%”，阈值为0，即只要经过函数“most50%”聚合并经负面信息修正后的像元值大于0，就将其纳入过火区，并作为新的“种子”像元继续生长，直到生长结束。

3.3 精度验证

因为缺少现场实测的过火区边界数据，得到过火区提取结果后，在原始卫星图像上通过目视解译选取样本点进行了精度验证(分别选取了过火区和未过火区样本各大约1 700个样本点，对得到的过火区提取结果进行精度验证)。图5为过火区结果矢量边界与B4(R)，B3(G)，B2(B)假彩色合成图像的叠置显示。

图5 过火区矢量边界与B4(R)，B3(G)，B2(B)假彩色合成图像叠置显示Fig.5 Overlapping display of burned area boundary and HJ image composed of B4(R)，B3(G)，B2(B)

从图5可以看出，在所提取的过火区中，误判非常少;而由于在“种子”像元提取阶段采用了比较严格的语义量化函数和阈值，致使少量燃烧程度较轻的小面积过火区的孤立像元群被漏判。但从总体上看，大部分过火区都通过上述算法被成功地识别出来。表4给出了2个验证样本的过火区制图精度。

表4 过火区制图精度Tab.4 Accuracy of burned area mapping (%)

从表4可以看出，使用本文算法可以有效地对火烧迹地进行提取，其漏判率和误判率都较低，过火区用户精度和总体精度都在85%以上。该方法解决了光谱指数的优选问题，因为没有一个光谱指数对所有区域或一个区域的所有像元都是最优的;在“种子”像元选取时，采用模糊分类的算法，逐像元、自适应地对光谱指数进行了优选聚合，当只有90%的光谱指数都是突出过火区信息时，才认为这个像元是“种子”像元，从而提高了在不同区域的适用性和强健性;模糊分类没有采用一个简单的阈值进行硬性分类，而是通过对每个像元与类别属性的相似度进行分类，从而降低了使用硬性阈值划分类别带来的误判或漏判。

4 结论与展望

1)以环境减灾一号小卫星B星(HJ－1B)CCD和IRS图像为数据源，采用模糊分类方法对森林过火区进行提取。通过目视解译对森林过火区提取结果进行精度验证的结果表明，本文提出的方法效果较好，能够满足应用需求。

2)采用加入了温度参数的负面信息对“种子”像元和边界图像进行修正，修正结果显著降低了误判率。在单时相的森林过火区制图中，用可见光和近红外波段数据难以区分水体和阴影，而过火区在温度上高于水体和阴影，说明了温度对过火区提取的重要性。

3)本文算法逐像元地对光谱指数进行聚合，对不同区域的像元能够自适应地突出森林过火区信息并抑制冗余信息;并且模糊分类对每个像元的归属不采用硬性的阈值，不仅解决了不同区域参数优选的问题，而且在一定程度上解决了阈值设置的问题。因此，从单时相制图的角度看，本文方法能够满足快速、跨区域森林过火区遥感制图的适用性和精度要求。

4)单时相数据易于获取、计算量小，但是缺少火灾前后波段信息对比的优势。为了提高精度、补偿缺陷，本文对短波红外和热红外波段数据进行了从低分辨率向高分辨率的重采样，这样也容易带来几何配准和类别边界像元归属的误差。对于大尺度的森林过火区遥感制图，将CCD数据的分辨率从30 m重采样到300 m可能是更好的选择。

[1]Kontoes C C，Poilvé H，Florsch G，et al.A comparative analysis of a fixed thresholding vs a classification tree approach for operational burn scar detection and mapping[J].International Journal of Applied Earth Observation and Geoinformation，2009，11(5):299－316.

[2]Smith A M S，Drake N A，Wooster M J，et al.Production of Landsat ETM+reference imagery of burned areas within Southern African savannahs:Comparison of methods and application to MODIS[J].International Journal of Remote Sensing，2007，28(12):2753－2775.

[3]Petropoulos G P，Vadrevu K P，Xanthopoulos G，et al.A comparison of spectral angle mapper and artificial neural network classifiers combined with Landsat TM imagery analysis for obtaining burnt area mapping[J].Sensors，2010，10(3):1967－ 1985.

[4]Henry M C.Comparison of single－ and multi－ date Landsat data for mapping wildfire scars in Ocala National Forest，Florida[J].Photogrammetric Engineering and Remote Sensing，2008，74(7):881－891.

[5]Cassidy L.Mapping the annual area burned in the wetlands of the Okavango panhandle using a hierarchical classification approach[J].Wetlands Ecology Management，2007，15(4):253－ 268.

[6]Fraser R H，Li Z，Cihlar J.Hotspot and NDVI differencing synergy(HANDS):A new technique for burned area mapping over boreal forest[J].Remote Sensing of Environment，2000，74(3):362－376.

[7]Stroppiana D，Bordogna G，Carrara P，et al.A method for extracting burned areas from Landsat TM/ETM+images by soft aggregation of multiple spectral indices and a region growing algorithm[J].ISPRS Journal of Photogrammetry and Remote Sensing，2012，69:88－102.

[8]Veraverbeke S，Harris S，Hook S.Evaluating spectral indices for burned area discrimination using MODIS/ASTER(MASTER)airborne simulator data[J].Remote Sensing of Environment，2011，115(10):2702－2709.

[9]Bastarrika A，Chuvieco E，Martín M P.Mapping burned areas from Landsat TM/ETM+data with a two－phase algorithm:Balancing omission and commission errors[J].Remote Sensing of Environment，2011，115(4):1003－1012.

[10]Stroppiana D，Bordogna G，Carrara P，et al.Positive and negative information for assessing and revising scores of burn evidence[J].IEEE Geoscience and Remote Sensing Letters，2012，9(3):363－367.

[11]中国资源卫星应用中心.HJ卫星参数介绍[EB/OL].[2009].http://www.cresda.com/n16/index.html.China Centre for Resources Satellite Data and Application.HJCCD sensor technical specifications[EB/OL].[2009].http://www.cresda.com/n16/index.html.

[12]熊文成，魏 武，孙中平，等.基于环境一号卫星B星CCD与红外相机融合的澳洲火灾监测[J].遥感技术与应用，2010，25(2):178－182.Xiong W C，Wei W，Sun Z P，et al.Australian forest fire disaster monitoring based on CCD and IRS data of HJ－1－ B[J].Remote Sensing Technology and Application，2010，25(2):178－182.

[13]丁 凤，徐涵秋.TM热波段图像的地表温度反演算法与实验分析[J].地球信息科学，2006，8(3):125－130.Ding F，Xu H Q.Comparison of two new algorithms for retrieving land surface temperature from Landsat TM thermal band[J].Geo－ Information Science，2006，8(3):125－130.

[14]Rouse J W，Haas R H，Schell J A，et al.Monitoring vegetation systems in the Great Plains with ERTS[C].Proc.Third ERTS Symposium，1973，NASA SP－351，1:309－317.

[15]Huete A，Didan K，Miura T，et al.Overview of the radiometric and biophysical performance of the MODIS vegetation indices[J].Remote Sensing of Environment，2002，83(2):195－213.

[16]Huete A.A soil－ adjusted vegetation index(SAVI)[J].Remote Sensing of Environment，1998，25(3):295－309.

[17]Smith A M S，Wooster M J，Drake N A，et al.Testing the potential of multi－spectral remote sensing for retrospectively estimating fire severity in African savanna hs[J].Remote Sensing of Environment，2005，97(1):92－115.

[18]Key C H，Benson N C.Measuring and remote sensing of burn severity[C]//Neuenschwander L F，Ryan K C.(Eds.)Proc.Joint Fire Science Conference and Workshop，1999，2:15－17.

[19]毕 翔，韩江洪，刘征宇.基于多特征相似性融合的隶属度函数研究[J].电子测量与仪器学报，2011，25(10):835－841.Bi X，Han J H，Liu Z Y.Research of membership function based on fusion of multi－ feature similarity[J].Journal of Electronic Measurement and Instrument，2011，25(10):835－841.

[20]Yager R R.On ordered weighted averaging aggregation operators in multi－ criteria decision making[J].IEEE Transactions on Systems，Man，and Cybernetics，1988，18(1):183－190.

[21]Yager R R.Families of OWA operators[J].Fuzzy Sets and Systems，1993，59(2):125－148.