基于時(shí)空融合技術(shù)的森林火災(zāi)遙感動態(tài)監(jiān)測

黃武彪，欒海軍，李大成

（1.廈門理工學(xué)院計(jì)算機(jī)與信息工程學(xué)院，福建廈門 361024；2.長安大學(xué)地質(zhì)工程與測繪學(xué)院，陜西西安 710064；3.廈門理工學(xué)院數(shù)字福建自然災(zāi)害監(jiān)測大數(shù)據(jù)研究所，福建廈門 361024；4.太原理工大學(xué)礦業(yè)工程學(xué)院，山西太原 030024）

引言

森林火災(zāi)是一種突發(fā)性強(qiáng)、破壞性大、救援困難的自然災(zāi)害，全球每年平均發(fā)生森林火災(zāi)數(shù)十萬次，森林受災(zāi)面積達(dá)數(shù)百公頃［1］，對生態(tài)環(huán)境和經(jīng)濟(jì)發(fā)展造成了極大威脅。截止2019年，全國有記錄森林火災(zāi)2 345次，其中重大火災(zāi)為8次，特別重大火災(zāi)為1次，火場總面積為39 705 hm2，受害森林面積為13 505 hm2，傷亡人數(shù)為76人，其它損失折款16 219.9萬元［2］。森林火災(zāi)監(jiān)測，要求及時(shí)發(fā)現(xiàn)著火點(diǎn)的位置及其變化、過火面積，并準(zhǔn)確評估出火災(zāi)損失及影響，因此需要具有高時(shí)間分辨率和高空間分辨率的遙感影像來進(jìn)行分析判讀。但在現(xiàn)有傳感器硬件條件限制下，衛(wèi)星遙感數(shù)據(jù)無法同時(shí)滿足高空間分辨率和高時(shí)間分辨率的要求［3-4］。因此，很多學(xué)者提出時(shí)空融合的技術(shù)方法來解決遙感傳感器時(shí)間分辨率和空間分辨率的矛盾［5-8］，將該方法應(yīng)用于森林火災(zāi)場景中，便于更加準(zhǔn)確、快速得獲得受災(zāi)區(qū)域高空間和高時(shí)間分辨率的遙感影像，為森林火災(zāi)演化的監(jiān)測和災(zāi)損評估提供有力支撐。Gao等［9］提出的時(shí)空自適應(yīng)反射率融合模型（Spatial and Temporal Adaptive Reflectance Fusion Model，STARFM），用于融合Landsat影像和MODIS影像得到高時(shí)空分辨率的數(shù)據(jù)，取得了較好的效果，該模型適用于同質(zhì)地表季節(jié)變化情形。Hilker等［10］在提出了一種針對反射率變化的時(shí)空自適應(yīng)融合模型算法（Spatial TemporalAdaptive Algorithm for mapping Reflectance Change，STA?ARCH），從低分辨率影像的密集時(shí)間序列中檢測出變化點(diǎn)，以提高土地覆蓋類型變化時(shí)的STARFM性能，適用于地表反射率突發(fā)擾動事件情形。鄔明權(quán)等［11］提出基于混合像元分解的方法（Spatial and Temporal Data Fusion Model，STDFM）來融合MODIS和Landsat影像數(shù)據(jù)，進(jìn)一步，Zhang等［12］對STDFM方法進(jìn)行了改進(jìn)，提出增強(qiáng)型基于混合像元分解的方法（Enhanced spatial and temporal data fusion model，ESTDFM），該模型適用于異質(zhì)地表季節(jié)變化情形。Zhu等［13］基于STARFM提出了一種增強(qiáng)型時(shí)空自適應(yīng)反射率融合模型（En?hanced Spatial and Temporal Adaptive Reflectance Fusion Model，ESTARFM），引入了一個轉(zhuǎn)換系數(shù)，可以更好地預(yù)測異質(zhì)性地表的反射率的變化。Huang等［14］提出的基于稀疏表示的時(shí)空反射率融合模型（SParse?repre?sentation?based SpatioTem?poral reflectance FusionModel，SPSTFM），將稀疏表達(dá)理論引入時(shí)空融合算法，該模型適用于地表季節(jié)和類別變化情形。除地表反射率參數(shù)外，學(xué)者們［15-17］對地表溫度參數(shù)也進(jìn)行了時(shí)空融合研究。學(xué)者們致力于研發(fā)出通用性與魯棒性更優(yōu)的時(shí)空融合算法。如，Cheng等［18］提出了一種時(shí)空非局部濾波融合模型（Spatial and Temporal Nonlocal Filter?Based Data FusionMethod，STNLFFM），對于異質(zhì)性地表區(qū)域具有更高的預(yù)測精度。Zhao等［19］提出了一種針對復(fù)雜地表變化的魯棒性自適應(yīng)時(shí)空數(shù)據(jù)融合模型（Ro?bust Adaptive Spatial and Temporal Fusion Model，RASTFM），在捕捉地表變化現(xiàn)象時(shí)具有更高的準(zhǔn)確度和魯棒性。黃波和姜曉璐［20］提出一種增強(qiáng)型空間像元分解時(shí)空遙感影像融合算法（Unmixing Enhanced model for Spatial and temporal image fusion，EUSTFM），能夠?qū)崿F(xiàn)對季節(jié)性變化及復(fù)雜的地物類型變化的穩(wěn)定預(yù)測，生成具有更高精度的融合影像。

鑒于森林屬于地表相對均質(zhì)的場景，故本次研究將選用經(jīng)典的STARFM算法與具有組合優(yōu)勢的基于地物內(nèi)組分時(shí)相變化模型的地表反射率時(shí)空融合算法，并綜合使用多種空間分辨率更優(yōu)（≤30m）的國內(nèi)外傳感器影像，對兩種算法所得時(shí)空融合結(jié)果優(yōu)化組合，進(jìn)而應(yīng)用于2019年3月30日四川涼山州木里縣森林火災(zāi)遙感動態(tài)監(jiān)測中。

1 研究區(qū)與研究數(shù)據(jù)

1.1 研究區(qū)概況

本次研究區(qū)域（如圖1所示）為四川涼山木里縣雅礱江鎮(zhèn)立爾村附近的雅礱江邊海拔約3 800 m的一處森林中，區(qū)域經(jīng)緯度介于28°32′10″N～28°33′27″N，101°15′02″E～101°16′49″E之間。該區(qū)域位于典型的高山峽谷區(qū)，氣候特點(diǎn)為冷熱兩季交替、干濕分明，天氣干燥，日照強(qiáng)，大面積被植被覆蓋，有少部分水及巖石。地形復(fù)雜、坡陡谷深，交通、通訊不便。火災(zāi)發(fā)生時(shí)間為2019年3月30日-4月4日。此次火災(zāi)造成30名消防員犧牲，據(jù)中國新聞網(wǎng)2019年4月5日報(bào)道，本次火災(zāi)起火點(diǎn)為一顆位于山脊上的樹齡約80年的云南松，起火原因?yàn)槔讚艋穑?1-22］。

圖1 研究區(qū)位置Fig.1 Location of the study area

1.2 研究數(shù)據(jù)

本次研究所用實(shí)驗(yàn)數(shù)據(jù)源為MOD09GA、Landsat8 OLI、Sentinel-2及GF-1 WFV影像，具體實(shí)驗(yàn)影像說明見表1，表中所列影像均為本次實(shí)驗(yàn)中所用到的數(shù)據(jù)。各遙感數(shù)據(jù)的特征與使用波段如表2所示。

表1 實(shí)驗(yàn)數(shù)據(jù)說明Table 1 Experimental data

表2 各遙感數(shù)據(jù)的特征及使用波段Table 2 The characteristics and use bands of various types of remote sensing data

Landsat8 OLI數(shù)據(jù)和Sentinel-2數(shù)據(jù)可從USGS（https：//earthexplorer.usgs.gov/）免費(fèi)獲取，MOD09GA數(shù)據(jù)可從NASA（https：//ladsweb.modaps.eosdis.nasa.gov/search/）免費(fèi)獲取，GF-1 WFV數(shù)據(jù)可從中國資源衛(wèi)星應(yīng)用中心陸地觀測衛(wèi)星數(shù)據(jù)服務(wù)平臺（http：//218.247.138.119：7777/DSSPlatform/productSearch.html）免費(fèi)獲取。

在進(jìn)行時(shí)空融合之前，需要對獲取的影像數(shù)據(jù)使用遙感圖像處理軟件進(jìn)行預(yù)處理。使用ENVI對Land?sat8 OLI影像進(jìn)行輻射定標(biāo)和大氣校正，生成反射率數(shù)值影像。Sentinel-2影像使用Sen2Cor和SNAP軟件進(jìn)行大氣校正并轉(zhuǎn)換成ENVI所能打開的img格式。GF-1WFV影像除進(jìn)行與Landsat8 OLI相同的預(yù)處理外，還需先使用ENVI遙感圖像處理軟件進(jìn)行幾何糾正，其輻射定標(biāo)時(shí)絕對輻射定標(biāo)系數(shù)在中國資源衛(wèi)星應(yīng)用中心網(wǎng)站上下載。對于MOD09GA影像需先使用MODIS Reprojection Tool（MRT）重投影至UTM/WGS84坐標(biāo)系、GeoTIFF格式，并重采樣為10 m、16 m和30 m分辨率。除上述預(yù)處理外，還應(yīng)使用ENVI遙感圖像處理軟件進(jìn)行各傳感器影像裁剪，確保其研究區(qū)域范圍相同。在本次研究中，均針對各影像的單波段進(jìn)行研究。

2 研究方法

具體研究技術(shù)路線圖如圖2所示。主要技術(shù)流程闡述如下：

圖2 本研究的技術(shù)路線圖Fig.2 Flowchart of the research

（1）對獲取到的MOD09GA、Landsat8 OLI、Sentinel-2、GF-1WFV影像數(shù)據(jù)進(jìn)行預(yù)處理；

（2）使用STARFM算法與基于地物內(nèi)組分時(shí)相變化模型的地表反射率時(shí)空融合算法分別對MOD09GA與Landsat8 OLI、Sentinel-2、GF-1 WFV影像進(jìn)行融合，以達(dá)到結(jié)合多源影像的目的，從而生成待預(yù)測時(shí)刻的中等空間分辨率影像；

（3）基于預(yù)測影像計(jì)算火災(zāi)指標(biāo)因子，進(jìn)行火災(zāi)演化趨勢分析。

2.1 面向多源數(shù)據(jù)的時(shí)空自適應(yīng)反射率融合模型

時(shí)空自適應(yīng)反射率融合模型（STARFM）在使用前需要對來自不同平臺的觀測數(shù)據(jù)進(jìn)行幾何配準(zhǔn)，并進(jìn)行大氣校正轉(zhuǎn)化至地表反射率。STARFM算法的總體流程描述如下［9］。STARFM算法的基本原理是：假定在tk日低分辨率影像數(shù)據(jù)地表反射率值M（xi，yj，t k）與中高分辨率影像數(shù)據(jù)地表反射率值L（xi，yj，tk）之間的關(guān)系可以表示為

則在t0日的低分辨率影像數(shù)據(jù)地表反射率值M（xi，yj，t0）與在t0日的中高分辨率影像數(shù)據(jù)地表反射率預(yù)測值L（xi，yj，t0）的關(guān)系可以表示為

εk與ε0分別表示觀測到的低分辨率影像數(shù)據(jù)和中高分辨率影像數(shù)據(jù)地表反射率之間由于不同的波段寬度和太陽幾何造成的差異。假設(shè)在預(yù)測日期t0和日期tk時(shí)像素（xi，y j）的地面覆蓋類型和系統(tǒng)誤差不變，則ε0=εk，因此

但這僅僅是一種理想情況，它們之間的關(guān)系受到以下3方面的影響：1）觀測的低分辨率影像數(shù)據(jù)在用與高分辨率影像數(shù)據(jù)相同的空間分辨率考慮時(shí)可能包括混合的土地覆蓋類型；2）在預(yù)測期間，土地覆蓋可能從一種類型變?yōu)榱硪环N類型；3）土地覆蓋狀態(tài)和太陽幾何雙向反射率分布函數(shù)的變化將改變從預(yù)測日期t0到日期tk的反射率。

因此，通過引入鄰近像素的附加信息，使用加權(quán)函數(shù)來計(jì)算日期t0時(shí)中心像素的地表反射率：

其中ω是搜索窗口的大小是該移動窗口的中心像素。為確保使用來自鄰近像素的正確信息，僅使用來自同一光譜類別的和來自移動窗口內(nèi)中高分辨率影像數(shù)據(jù)地表反射率的無云像素來計(jì)算反射率。權(quán)重Wijk代表每個鄰近像素對中心像素的預(yù)測反射率的貢獻(xiàn)程度，從光譜、時(shí)間、距離3個方面來確定每個光譜相似像素的最終權(quán)重。

本次研究所做改進(jìn)：上述內(nèi)容通常針對單一類型中等空間分辨率傳感器，當(dāng)面向多源數(shù)據(jù)時(shí)，需要對STARFM時(shí)空融合模型進(jìn)行優(yōu)化組合應(yīng)用。具體優(yōu)化組合策略設(shè)計(jì)如下：在經(jīng)典的STARFM算法基礎(chǔ)上，在30 m Landsat衛(wèi)星16天或更長的預(yù)測周期內(nèi)利用多種國內(nèi)外傳感器影像（空間分辨率≤30 m）進(jìn)行分段，分段原則為：時(shí)間最鄰近及空間分辨率優(yōu)先。上述原則以案例描述為：（1）若引入10 m Sentinel-2可見光影像，該影像在Landsat基期與末期時(shí)間段之間（如第5天），則將預(yù)測區(qū)間分為兩段（Landsat基期-第5天）、（第6天-第16天）；（2）對于每一時(shí)間段（包含新的基期和末期兩天），若天數(shù)為偶數(shù)，對其均分，前半段影像預(yù)測由本時(shí)間段基期影像完成，后半段影像預(yù)測由本時(shí)間段末期影像完成，若天數(shù)為奇數(shù)，前半段（如Landsat基期-第2天）的影像預(yù)測由“最鄰近”的Landsat基期影像完成，后半段（如第4天-第5天）的影像預(yù)測由“最鄰近”的Sentinel-2影像完成，中間一天（如Landsat基期-第5天時(shí)間段，中間一天為第3天）的影像預(yù)測依照空間分辨率優(yōu)先原則，由空間分辨率高的傳感器影像（也就是Sentinel-2影像）完成；（3）當(dāng)引入的傳感器影像更多時(shí)，參照上述案例進(jìn)行擴(kuò)展即可。各時(shí)間分段內(nèi)，當(dāng)新的基期或末期影像中高空間分辨率優(yōu)于30 m，需要通過現(xiàn)有商業(yè)遙感軟件進(jìn)行空間重采樣至30 m，然后進(jìn)行預(yù)測影像生成。

2.2 面向多源數(shù)據(jù)的基于地物內(nèi)組分時(shí)相變化模型的地表反射率時(shí)空融合算法

對于相同的土地覆蓋類型的均質(zhì)像素都具有相同的反射率，并且這些像素的季節(jié)性和雙向反射率變化對于不同的像元大小也應(yīng)該相同［23］。對于不同分辨率的影像，每種土地覆蓋類型在獲取日期和預(yù)測日期之間的變化關(guān)系大致相同，即在不同組分之間具有空間尺度不變性。因此，可以將基于獲取日期的低分辨率影像數(shù)據(jù)和中高分辨率影像數(shù)據(jù)建立的變化關(guān)系映射到中高分辨率遙感影像中，從而預(yù)測出待預(yù)測時(shí)刻的中高分辨率影像。

基于上述前提，提出了基于地物內(nèi)組分時(shí)相變化模型的地表反射率時(shí)空融合算法，具體流程如圖3所示。先對中高空間分辨率遙感影像采用迭代自組織數(shù)據(jù)分析算法（ISODATA）進(jìn)行非監(jiān)督分類，將其分為多個地物類別，得到不同的地物組分。然后分別對各個地物組分的變化，建立相應(yīng)時(shí)相變化模型進(jìn)行預(yù)測。

圖3 基于地物內(nèi)組分時(shí)相變化模型的地表反射率時(shí)空融合算法流程圖Fig.3 Flow chart of Spatiotemporal fusion algorithm of land surface reflectance basedon temporal phase change model of components in surface features

尋找“純像元”的基礎(chǔ)為將非監(jiān)督分類得到的結(jié)果聚合到低分辨率影像像元大小，其聚合過程可理解為通過計(jì)算聚合前后的分辨率像元的比例n，將非監(jiān)督分類影像自左上開始，組成若干n*n大小的矩陣，并記錄聚合后的行列號。分別計(jì)算在每一個n*n矩陣中各組分所占比例，取比例最大的組分作為該像元的類別，并判斷比例值是否大于20%（20%為經(jīng)驗(yàn)值，針對不同的分辨率設(shè)置不同的閾值），若大于該值，則將該聚合后的像元稱為該地類的一個“純像元”。依次得到各類地物的“純像元”，并找出兩個時(shí)間的低分辨率影像中對應(yīng)組分內(nèi)部的無云像元是如何變化的，使用差值、比值或變化率建立不同地物組分的時(shí)相變化模型。整個過程基于MATLAB軟件實(shí)現(xiàn)。

基于地物內(nèi)組分時(shí)相變化模型的地表反射率時(shí)空融合算法三類模型依次描述如下。

差值模型：

比值模型：

變化率模型：

式中，M OD09G A_T0表示在T0時(shí)刻低分辨率遙感影像的地表反射率值，M O D09 G A_T1表示在T1時(shí)刻低分辨率遙感影像的地表反射率值，S entinel_T0表示在T0時(shí)刻的中高分辨率遙感影像地表反射率值，pre_Sen?tinel_T1表示預(yù)測得到的T1時(shí)刻的中高分辨率遙感影像地表反射率值，mean表示對所有的值取平均。

本次研究所做改進(jìn)：上文所述算法通常針對單一類型中等空間分辨率傳感器，對于多源數(shù)據(jù)的實(shí)施原則可參照2.1節(jié)“優(yōu)化組合策略”，其實(shí)施過程相似，但無需“重采樣至30 m”步驟。

2.3 火災(zāi)指標(biāo)因子介紹

在森林火災(zāi)監(jiān)測中應(yīng)用遙感影像，需要對使用兩種時(shí)空融合算法得到的影像數(shù)據(jù)計(jì)算燃燒面積指數(shù)和歸一化燃燒指數(shù)。

2.3.1 燃燒面積指數(shù)

燃燒面積指數(shù)（Burn Area Index（BAI））是采用影像的紅色波段（Red）和近紅外波段（N I R）來增強(qiáng)火燒以后的地表信息，也就是增強(qiáng)火災(zāi)過火以后圖像上的木炭信號，燃燒區(qū)域的BAI值比較大。計(jì)算公式［24］如下：

在ENVI5.3中利用光譜指數(shù)計(jì)算工具選擇“Burn Area Index”即可計(jì)算研究區(qū)燃燒面積指數(shù)。

2.3.2 歸一化燃燒指數(shù)

歸一化燃燒指數(shù)（Normalized Burn Ratio（N B R））是基于近紅外波段（N IR）和短波紅外波段2（SWI R2）來增強(qiáng)較大范圍的火災(zāi)區(qū)域，燃燒區(qū)域的NB R值比較小，計(jì)算公式［25-26］如下：

在ENVI5.3中利用光譜指數(shù)計(jì)算工具選擇“Normalized Burn Ratio”即可計(jì)算研究區(qū)歸一化燃燒指數(shù)。

3 分析與討論

3.1 融合結(jié)果分析

利用STARFM算法和基于地物內(nèi)組分時(shí)相變化模型的地表反射率時(shí)空融合算法分別對MOD09GA與Landsat8 OLI、Sentinel-2及GF-1 WFV進(jìn)行時(shí)空融合。根據(jù)融合結(jié)果（見表3），將預(yù)測后的影像與可從其他類型傳感器獲得的同時(shí)刻的影像進(jìn)行目視對比，分析本次實(shí)驗(yàn)所采用的時(shí)空融合算法的精度，并對火災(zāi)發(fā)生前后的火災(zāi)指標(biāo)因子進(jìn)行提取分析。如圖4是3月22日GF-1影像與MOD09GA影像預(yù)測3月21日GF-1影像，并與Sentinel-2影像對比，圖5是3月18日Landsat8 OLI影像與MOD09GA影像預(yù)測3月26日Landsat8 OLI影像，并與Sentinel-2影像對比，圖6是3月26日Sentinel-2影像與MOD09GA影像預(yù)測3月30日Senti?nel-2影像，并與GF-1影像對比。

圖4 使用3月22日GF-1與MODIS影像預(yù)測3月21日GF-1影像，并與Sentinel-2影像對比Fig.4 Use GF-1 and MODIS images on March 22 to predict GF-1 image on March 21 and compare with Sentinel-2 images

圖5 3月18日Landsat與MODIS影像預(yù)測3月26日Landsat影像，并與Sentinel-2影像對比Fig.5 Use Landsat and MODIS images on March 18 to predict Landsat image on March 26 and compare with Sentinel-2 images

圖6 3月26日Sentinel-2與MODIS影像預(yù)測3月30日Sentinel-2影像，并與GF-1影像對比Fig.6 Use Sentinel-2 and MODIS images on March 26 to predict Sentinel-2 image on March 30 and compare with GF-1 images

表3 使用兩種算法的融合結(jié)果Table 3 Fusion results based on the two algorithms

兩種算法在不同類型遙感數(shù)據(jù)融合應(yīng)用中各有其局限性。將時(shí)空融合后的單波段影像進(jìn)行疊加，得到RGB彩色圖像，如圖4～6所示。基于圖4～6分析：使用STARFM算法對地表變化明顯區(qū)域預(yù)測精度差，從圖4（e）可發(fā)現(xiàn)影像左側(cè)存在黑色區(qū)域，與實(shí)際地表存在明顯差距，說明STARFM算法對GF-1 WFV的預(yù)測周邊裸露巖石的效果不是很好；從圖5（e）、圖6（e）與圖5（d）、圖6（d）的山脈、植被、巖石等地區(qū)進(jìn)行對比可以發(fā)現(xiàn)STARFM對于Landsat8 OLI和Sentinel-2的預(yù)測效果較好，但Landsat8 OLI的結(jié)果更優(yōu)。基于地物內(nèi)組分時(shí)相變化模型的地表反射率時(shí)空融合算法包括3種變化模型，將圖4（f）、圖5（f）和圖6（f）與對應(yīng)的（d）影像選取山脈、植被、巖石等地區(qū)對比，發(fā)現(xiàn)使用地物組分時(shí)相變化差值模型來進(jìn)行預(yù)測時(shí)，對Landsat8 OLI、Sentinel-2和GF-1預(yù)測效果都比較好；將圖4（g）、圖5（g）和圖6（g）與對應(yīng)的（d）影像對比，發(fā)現(xiàn)使用地物組分時(shí)相變化比值模型進(jìn)行預(yù)測時(shí)，圖4（g）、圖5（g）影像圖左側(cè)巖石地區(qū)呈現(xiàn)黃色，與真實(shí)地表存在明顯差異，因此地物組分時(shí)相變化比值模型對GF-1和Landsat8 OLI的周邊預(yù)測效果較差，對Sentinel-2的預(yù)測效果較好；將圖4（h）、圖5（h）和圖6（h）與對應(yīng)的（d）影像對比，發(fā)現(xiàn)使用地物組分時(shí)相變化的變化率模型進(jìn)行預(yù)測時(shí)，圖4（h）影像圖左側(cè)裸露巖石呈現(xiàn)明顯差異，因此地物組分時(shí)相變化的變化率模型對GF-1預(yù)測效果較差，Sentinel-2和Landsat8 OLI預(yù)測效果較好。綜合分析，對于Landsat8 OLI的預(yù)測可以選用STARFM和地物組分時(shí)相變化差值模型，對于Sentinel-2的預(yù)測，地物組分時(shí)相變化差值模型和變化率模型要好于STARFM和比值模型，對于GF-1的預(yù)測，盡量選擇地物組分時(shí)相變化差值模型。因此，通過影像對比與分析，基于地物內(nèi)組分時(shí)相變化的差值模型對于上述3種影像均適用，而使用比值模型進(jìn)行預(yù)測的效果都不是很好。

對不同傳感器中等空間分辨率影像與相同時(shí)間不同算法預(yù)測所得結(jié)果進(jìn)行比對分析，計(jì)算相同波段上的影像相關(guān)系數(shù)，結(jié)果如表4所示。表4中，相關(guān)系數(shù)越大表明算法在該類型傳感器該波段上的預(yù)測效果越好。可以看出，對于GF-1 WFV和Sentinel-2，在預(yù)測的每個波段，基于地物組分時(shí)相變化的差值模型的相關(guān)系數(shù)最高；對于Landsat8 OLI，在近紅外和短波紅外1兩個波段上，基于地物組分時(shí)相變化的變化率模型更優(yōu)。同時(shí)，不同的算法在不同傳感器影像的同一波段上的預(yù)測效果也不同，具體為：1）對于藍(lán)光波段，STARFM、基于地物組分時(shí)相變化的比值模型和變化率模型對Sentinel-2的預(yù)測效果最好，差值模型對GF-1 WFV的預(yù)測效果最優(yōu)；2）對于綠光波段，STARFM和基于地物組分時(shí)相變化的比值模型對Sentinel-2的預(yù)測效果最好，差值模型對GF-1 WFV的預(yù)測效果最優(yōu)，變化率模型對Landsat8 OLI的預(yù)測效果更好；3）對于紅光波段，STARFM、基于地物組分時(shí)相變化的差值模型、比值模型和變化率模型分別對Landsat8 OLI、GF-1 WFV、Sentinel-2和Landsat8 OLI的預(yù)測效果最好；4）對于近紅外波段，四種算法的預(yù)測效果最好的傳感器影像分別為Landsat8 OLI、GF-1 WFV、Landsat8 OLI、Landsat8 OLI；5）對于短波紅外1和2，4種算法在Landsat OLI上的預(yù)測效果均最好。對于提取2種火災(zāi)指標(biāo)因子所用波段，可選用差值模型進(jìn)行預(yù)測，其對幾個波段均適用，其次是變化率模型，最后是比值模型和STARFM模型。

表4 各波段不同時(shí)空融合方法的相關(guān)系數(shù)Table 4 Correlation coefficients of different fusion methods in each band

3.2 火災(zāi)變化時(shí)序分析

對2種時(shí)空融合算法所得結(jié)果進(jìn)行優(yōu)化組合，獲取質(zhì)量更優(yōu)的逐日中等空間分辨率預(yù)測影像，進(jìn)而對預(yù)測影像計(jì)算火災(zāi)指標(biāo)因子，由于GF-1 WFV影像無短波紅外波段，故無法計(jì)算歸一化燃燒指數(shù)，以“時(shí)間最鄰近及空間分辨率優(yōu)先”為原則選取符合條件的Landsat8 OLI/Sentinel-2預(yù)測影像作為替代計(jì)算該指數(shù)。經(jīng)過計(jì)算，如圖7所示是火災(zāi)發(fā)生前期、中期、后期3個時(shí)間的燃燒面積指數(shù)變化，圖8所示是火災(zāi)發(fā)生前期、中期、后期3個時(shí)間的歸一化燃燒指數(shù)對比。

圖7 火災(zāi)發(fā)生前期、中期、后期的燃燒面積指數(shù)Fig.7 Burn Area Index before，during and after the fire

圖8 火災(zāi)發(fā)生前期、中期、后期的歸一化燃燒指數(shù)Fig.8 Normalized Burn Ratio before，during and after the fire

圖7（a）中植被覆蓋區(qū)域呈現(xiàn)黑色或暗灰色，圖7（c）中植被覆蓋區(qū)域（未發(fā)生火災(zāi)區(qū)域）呈現(xiàn)黑色，而火災(zāi)發(fā)生區(qū)域呈現(xiàn)白色，顯著的對比度差異利于確定火災(zāi)發(fā)生以后的位置和面積；但圖7（b）研究區(qū)呈現(xiàn)黑色，無法看到燃燒區(qū)域，這是由于在火災(zāi)發(fā)生時(shí)有植被燃燒產(chǎn)生的濃煙，而燃燒面積指數(shù)計(jì)算所使用的紅波段（Red）和近紅外波段（NIR）均無法穿透煙霧，所以對于燃燒面積指數(shù)會產(chǎn)生一定的影響。圖8（a）中植被覆蓋區(qū)域呈現(xiàn)白色或淺灰色，圖8（c）中植被覆蓋區(qū)域（未發(fā)生火災(zāi)區(qū)域）呈現(xiàn)白色或淺灰色，而火災(zāi)發(fā)生區(qū)域呈現(xiàn)黑色，明顯的黑白對比可確定火災(zāi)發(fā)生以后的位置和面積；圖8（b）中可以看到火災(zāi)發(fā)生時(shí)零星的燃燒區(qū)域，這是由于歸一化燃燒指數(shù)計(jì)算所使用的短波紅外波段（SWIR）可以穿透煙霧，GF-1 WFV影像沒有短波紅外波段，這也造成了GF-1在火災(zāi)指標(biāo)因子分析中的局限性。從圖7整體來看，未發(fā)生火災(zāi)區(qū)域仍然保持黑色或暗灰色，而發(fā)生火災(zāi)區(qū)域變?yōu)榘咨粓D8則與圖7相反，從整體來看未發(fā)生火災(zāi)區(qū)域仍然保持白色或淺灰色，而發(fā)生火災(zāi)區(qū)域變?yōu)楹谏Ｇ铱傮w而言，火災(zāi)演化進(jìn)程中NBR的變化更加顯著，BAI效果略差。通過兩種火災(zāi)指標(biāo)因子的對比研究，結(jié)果表明：基于時(shí)空融合影像分析火災(zāi)演化態(tài)勢時(shí)，歸一化燃燒指數(shù)計(jì)算結(jié)果更敏感、更有效。

火災(zāi)發(fā)生時(shí)間為2019年3月30日-2019年4月4日。基于預(yù)測得到的影像并計(jì)算NBR，統(tǒng)計(jì)過火面積，我們發(fā)現(xiàn)：由于2019年3月30日影像獲取時(shí)火災(zāi)還未發(fā)生，故NBR結(jié)果中沒有火災(zāi)區(qū)域；2019年3月31日過火面積約為72 087.96m2；2019年4月1日過火面積約為187 616.86m2；2019年4月2日過火面積約為249 256.77m2；2019年4月3-5日MODIS數(shù)據(jù)云量太大無法使用，故無預(yù)測結(jié)果；2019年4月6日過火面積約為335 834.01m2。可以看到從2019年3月31日開始，火災(zāi)面積逐步擴(kuò)大，最終過火面積為335 834.01m2，這與官方所公布的結(jié)果總體一致。

3.3 討論

通過分析發(fā)現(xiàn)，結(jié)合多類型中高空間分辨率影像（如Landsat8 OLI、Sentinel-2和GF-1影像）、MODIS影像，并針對兩種時(shí)空融合算法結(jié)果優(yōu)化組合，可彌補(bǔ)單一時(shí)空融合方法或使用單一中等空間分辨率影像（如Landsat影像）和MODIS影像時(shí)空融合的不足，可獲取更精確、更好的時(shí)間和空間分辨率預(yù)測影像。但對于實(shí)驗(yàn)存在的問題進(jìn)行討論：在應(yīng)用時(shí)空融合時(shí)，我們發(fā)現(xiàn)在使用較差的預(yù)處理后的影像進(jìn)行融合時(shí)，效果比較差，因此，時(shí)空融合的2種算法對數(shù)據(jù)的預(yù)處理工作有很高的要求，包括輻射定標(biāo)、大氣校正等，對于高分系列衛(wèi)星還受到輻射校正和幾何校正精度的影響，因此時(shí)空融合結(jié)果的精度一定程度上受影像預(yù)處理精度的影響；其次，該算法對于輸入的影像數(shù)據(jù)有很高的要求，要求獲取的影像必須是無云的純凈像素，MODIS影像的時(shí)間分辨率雖然比較好，但不可避免的也會出現(xiàn)研究區(qū)被云層覆蓋的情況，如果是薄霧，可使用去云處理，但如果云層較厚，關(guān)于去云研究目前還未成熟。因此，在后續(xù)研究中，可以考慮加入雷達(dá)影像，利用雷達(dá)影像能夠穿透云層，可以看到火災(zāi)的邊界等信息的特點(diǎn)，探索雷達(dá)強(qiáng)度圖像與光學(xué)數(shù)據(jù)的之間的關(guān)系，實(shí)現(xiàn)兩者的融合，以期得到更好的應(yīng)用。

4 結(jié)論

本文綜合利用經(jīng)典的STARFM算法與具有組合優(yōu)勢的基于地物內(nèi)組分時(shí)相變化模型的地表反射率時(shí)空融合算法，聯(lián)合使用多種空間分辨率更優(yōu)（≤30m）的傳感器影像（Landsat8 OLI、Sentinel-2、GF-1 WFV）與MODIS影像，以“時(shí)間最鄰近及空間分辨率優(yōu)先”為原則對傳統(tǒng)單一中等空間分辨率影像預(yù)測周期（如Land?sat影像為16天）進(jìn)行分段獨(dú)立預(yù)測，并優(yōu)化組合兩種預(yù)測方法的預(yù)測結(jié)果，進(jìn)而對四川涼山木里縣3·30森林大火預(yù)測結(jié)果計(jì)算BAI和NBR指數(shù)，提取森林火災(zāi)區(qū)域動態(tài)變化信息。研究結(jié)果表明：本文所設(shè)計(jì)的時(shí)空融合策略可彌補(bǔ)單一時(shí)空融合方法或使用單一中等空間分辨率影像（如Landsat影像）和MODIS影像時(shí)空融合的不足，同時(shí)，使用NBR提取火災(zāi)變化信息效果更優(yōu)。本文所提出的方法在森林火災(zāi)遙感動態(tài)監(jiān)測場景中具有可行性，具有進(jìn)一步深入研究的價(jià)值與意義。