基于模糊集理論的單時相跨區域森林過火區遙感制圖

朱 曦，覃先林，廖 靖

(中國林業科學研究院資源信息研究所，北京 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 研究方法

本文研究的方法是基于一個區域生長的過程，該過程需要一個初始“種子”像元和一個限制生長準則?？傮w技術流程如圖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]，因為只有最大化地突出非過火像元、進而生成負面信息，通過負面信息修正剔除水體和陰影等易與過火像元混淆的類別，才能有效地減少誤判?！癕ax”算子加權向量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))對于水體陰影等都具有很好的可分性?；诳梢姽夂徒t外的光譜指數對過火區和植被具有一定的可分性[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.