劉 洋，玉 山，2，都瓦拉，楊曉穎

(1. 內(nèi)蒙古師范大學地理科學學院，內(nèi)蒙古 呼和浩特市 010022； 2. 內(nèi)蒙古自治區(qū)遙感與地理信息系統(tǒng)重點實驗室，內(nèi)蒙古 呼和浩特市 010022； 3. 中國農(nóng)業(yè)科學院草原研究所，內(nèi)蒙古 呼和浩特市 010022)

1 引言

我國草原飽受火災困擾，草原火災突發(fā)性強、蔓延速度快、危害性大。造成經(jīng)濟損失和人員傷亡的同時毀壞草場資源，生物多樣性降低，破壞草地生態(tài)系統(tǒng)的平衡。為了及時、準確地防控草原火災，有必要對草原火蔓延進行研究。

新一代靜止氣象衛(wèi)星時空分辨率得到大幅提高，其高頻次、覆蓋范圍廣的特點在監(jiān)測草原火災方面有較大的優(yōu)勢。目前，國內(nèi)外學者利用Himawari-8影像監(jiān)測草原火。在此基礎(chǔ)上，一些學者利用Himawari-8監(jiān)測成果結(jié)合王正非火蔓延模型與元胞自動機等計算機仿真算法，以柵格或矢量圖形來預測草原火蔓延，并取得了較好的結(jié)果。

GRAPES＿MESO是GRAPES(Global/Regional Assimila-tion and Prediction System)系統(tǒng)中區(qū)域中尺度數(shù)值預報模式，其在天氣預報業(yè)務(wù)和防災減災方面發(fā)揮著重要的作用。國內(nèi)學者利用其高時空分辨率的特點在氣象預警、模型構(gòu)建等方面進行研究與應(yīng)用。

綜合前人研究經(jīng)驗，本文以內(nèi)蒙古草原高火險地區(qū)為研究區(qū)域，首先利用Himawari-8與GRAPES＿MESO氣象預報數(shù)據(jù)特征，由時序線性外推法確定草原火蔓延初速度，然后耦合元胞自動機和王正非火蔓延模型對草原火蔓延速度和蔓延邊界進行預測，并應(yīng)用到實際案例驗證預測結(jié)果。

2 草原火蔓延模擬預測模型及算法

2.1 數(shù)據(jù)

2.1.1 Himawari-8數(shù)據(jù)介紹

Himawari-8是由日本氣象廳發(fā)射的新一代靜止氣象衛(wèi)星，該衛(wèi)星對地觀測頻率高達 10 min/次，空間分辨率最高達500m。其搭載的AHI(Advanced Himawari Imager)傳感器波段范圍覆蓋可見光到遠紅外，設(shè)有3個可見光通道分辨率為500m，13個近紅外和紅外通道分辨率為2km，其中第7(3.74-3.96μm)、14(11.0-11.30μm)波段的敏感對象為地表溫度、云頂溫度。能夠很好的滿足草原火蔓延過程中需要高時效、多頻次、高精確獲得火情信息的條件(https：∥www.eorc.jaxa.jp/ptree/index＿j.html)。

本文選取內(nèi)蒙古錫林郭勒盟東烏珠穆沁旗薩麥蘇木2016年3月29日2時-7時分辨率為2km的影像數(shù)據(jù)。

2.1.2 GRAPES＿MESO簡介

GRAPES＿MESO中國及周邊區(qū)域數(shù)值預報產(chǎn)品，其空間分辨率為10km，時間分辨率為3小時。預報時次為00時、12時，預報最高時效72小時，要素包括位勢高度、溫度、風的u向量、風的v向量、相對濕度、降水量等(http：∥data.cma.cn)。

本文利用風的u和v向量計算2016年3月29日00時次6小時預報時效的10m特定高度的風速和風向。

2.1.3 輔助數(shù)據(jù)

本文模型和算法主要基于柵格數(shù)據(jù)處理計算，模型的主要輸入?yún)?shù)來源于30m、500m和2000m空間分辨率的遙感產(chǎn)品和氣象數(shù)值預報產(chǎn)品。遙感數(shù)據(jù)獲取于美國航天航空局(National Aeronautics and Space Administration，NASA)的MODIS 500米分辨率地表覆蓋數(shù)據(jù)(https：∥ladsweb.modaps.eosdis.nasa.gov/search/order)；草地類型數(shù)據(jù)由中國農(nóng)科院草原研究所提供；DEM數(shù)據(jù)收集自美國航天飛機雷達地形測繪使命(Shuttle Radar Topography Mission，SRTM)30米分辨率1＿Arc＿Second SRTMDEM數(shù)據(jù)(https：∥earthexplorer.usgs.gov/)；可燃物修正系數(shù)參照文獻[10]中可燃物燃燒實驗結(jié)果獲取。

2.2 時序線性外推法

Hamawari-8衛(wèi)星對地觀測頻率高達 10 min/次，基于連續(xù)兩次觀測反演的火情得到各個方向上火場邊界的擴展距離，進而計算該間隔內(nèi)蔓延邊界上草原火蔓延初速度。假設(shè)蔓延初始時刻(設(shè)=0)、蔓延開始后時刻、蔓延開始后時刻的火情如圖1所示(坐標原點為初始火點)，則在某一方向上，到時間段內(nèi)火場邊界的蔓延速度則為

(1)

其中,為時刻方向上火場邊界與初始火點的距離。同理，在該方向上到時間段內(nèi)火場邊界的蔓延速度為：

(2)

一般地，在方向上，-1到時間段內(nèi)火場邊界的蔓延速度為：

(3)

圖1 連續(xù)三個觀測時刻的火場示意圖

得到蔓延速度后，即可預測當前時刻往后一定時間段內(nèi)的火蔓延情況。對于方向，+1時刻的火場邊界位置為

，+1=，*(+1-)

(4)

2.3 地理元胞自動機

地理元胞自動機(Geograph Cellular Automaton，GeoCA)具有許多特性，使其在地理時空建模中具有吸引力。利用簡單的局部規(guī)則，它們有可以模擬復雜的現(xiàn)象。此外，它們本質(zhì)上具有空間性，其結(jié)構(gòu)與各種數(shù)字源提供的地理空間數(shù)據(jù)集兼容。標準元胞自動機是一個四元組即

=(，，，)

(5)

式中：為正整數(shù)，為元胞的基本單元；為元胞自動機中元胞的狀態(tài)，每個元胞狀態(tài)是一個有限或無限的元素；為元胞自動機的鄰域，鄰域的狀態(tài)是決定該元胞是否進行轉(zhuǎn)換，鄰域的值通過土地覆蓋數(shù)據(jù)賦予；為元胞自動機空間中的轉(zhuǎn)換規(guī)則。當元胞自動機在二維空間上時，元胞空間結(jié)構(gòu)與柵格空間結(jié)構(gòu)極為相似，故以500m分辨率元胞自動機來模擬草原火蔓延。

根據(jù)草原火蔓延過程中各個階段，每個元胞的狀態(tài)用連續(xù)的整數(shù)值來表示()：=0：“不可燃燒”狀態(tài)是指此類元胞的屬性本身不具有可燃性，不具備向相鄰元胞擴散的能力。如沙地、水體、公路等；=1“未燃燒”狀態(tài)是指該元胞具有可燃性，由于草原火未將該元胞所表示單元點燃；=2：“燃燒”狀態(tài)是指元胞具備可燃性并且可以向其相鄰的擴散；=3：“熄滅”狀態(tài)是指該元胞具有可燃性并且該元胞已經(jīng)燃燒完畢。

草原火蔓延空間可以解釋為二維元胞空間，針對草原火蔓延提出以下元胞轉(zhuǎn)換規(guī)則：1：如果元胞(，)在時刻的狀態(tài)為2，該元胞相鄰元胞狀態(tài)不存在0或3 則其狀態(tài)根據(jù)式(6)計算。

2：如果元胞(，)在時刻的狀態(tài)為2，在+1時刻元胞(，)狀態(tài)轉(zhuǎn)變?yōu)?，則表示該元胞燒毀，元胞狀態(tài)值保持在3不變，并且不會向相鄰元胞擴散。

3：如果元胞(，)在時刻的狀態(tài)為2，但是相鄰元胞狀態(tài)值為0或3，則元胞(，)不會傳遞。

在+1時刻，元胞(，)狀態(tài)利用規(guī)則對相鄰元胞進行擴散，即

(6)

2.4 王正非模型

草原火蔓延過程中不同路徑上可燃物、氣象、地形要素的影響對草原火蔓延速度產(chǎn)生的差異，毛賢敏對王正非模型中風和地形因子的細化，結(jié)合文獻[10]可燃物修正系數(shù)的研究成果對模型進行修正。王正非火蔓延模型涉及四個參數(shù)：初始蔓延速度、可燃物配置格局修正系數(shù)、風力修正系數(shù)、地形坡度修正系數(shù)。

=

(7)

式中：為火蔓延速度m/min；為初始蔓延速度m/min；是可燃物配置格局修正系數(shù)(無因次)；為風力修正系數(shù)(無因次)；(1cos)為地形坡度修正系數(shù)(無因次)。

2.5 精度評價方法

251 蔓延范圍面積相對偏差Bias

從監(jiān)測到初始火點開始的一定時間段內(nèi)，火蔓延范圍的面積可由蔓延邊界包含在內(nèi)的像元面積相加而得。若模型模擬得到的火蔓延面積為S，參考數(shù)據(jù)得到的該面積為S′，則蔓延范圍面積的相對誤差Bias由式(8)計算得到。Bias越小，蔓延預測結(jié)果越精確，反之亦然。

Bias=(-′)′

(8)

2.5.2 蔓延邊界均方根誤差BRMSE(Boundary Root Mean Square Error)

(9)

式中BRMSE單位為像元(Pixel)，Res為像元分辨率。

2.5.3 像元平均絕對誤差CMAPE(Cell Mean Absolute Percentage Error)

模型模擬的像元總數(shù)為N，模型模擬區(qū)域與實際火情區(qū)域相交區(qū)域像元數(shù)為N，CMAPE為N與N的差的比值。CMAPE0%表示完美模型，CMAPE大于 100 %則表示劣質(zhì)模型。

(10)

3 模型求解與精度驗證

3.1 輸入數(shù)據(jù)預處理

草原火蔓延過程較為復雜，在算法實現(xiàn)過程中輸入數(shù)據(jù)的空間分辨率不同，草原火蔓延預測計算統(tǒng)一采用500m像元分辨率和WGS＿1984＿Albers投影系統(tǒng)。

首先，使用ArcGIS提取DEM數(shù)據(jù)中包含的坡度數(shù)據(jù)、利用IDL平臺將GRAPES＿MESO中國及周邊區(qū)域3小時數(shù)值預報產(chǎn)品中表示經(jīng)度方向上的風向量u和緯度方向上的風向量v計算得出風速和風向數(shù)據(jù)，將提取的坡度數(shù)據(jù)和風速數(shù)據(jù)分級賦值(見表1和表2)。將Himawari-8數(shù)據(jù)經(jīng)過幾何校正、重投影、圖像配準得到更為精準的分析數(shù)據(jù)，發(fā)生火災時中紅外第7(3.74-3.96μm)波段在時相上的亮溫變化明顯，亮溫會迅速升高，而第14(11.0-11.30μm)通道亮溫變化較小，產(chǎn)生背景亮溫差。因此根據(jù)時相變化設(shè)置閾值檢測火點，通過運算獲得影像過火區(qū)以及火點，作為草原火蔓延的現(xiàn)勢火情數(shù)據(jù)。

保證草原火蔓延模型與輔助數(shù)據(jù)之間的兼容性，算法統(tǒng)一采用IDL語言編寫，包括草原火蔓延主程序和地形數(shù)據(jù)、氣象數(shù)據(jù)標準化函數(shù)。程序首先讀取Himawari-8數(shù)據(jù)，之后進行數(shù)據(jù)標準化和分級賦值處理。地形數(shù)據(jù)根據(jù)在相關(guān)研究中將坡度分劃為6個方向的修正公式，根據(jù)研究區(qū)范圍內(nèi)DEM數(shù)據(jù)，通過修正公式計算獲得分級數(shù)值(見表1)，風力等級標準根據(jù)國家風力等級劃分標準賦值(見表2)，可燃物數(shù)據(jù)與草地類型數(shù)據(jù)一一對應(yīng)賦值分級。

表1 地形坡度修正系數(shù)(值)

表2 風力修正系數(shù)(值)

我國北方草原當年10月至來年4月容易發(fā)生火災。在此基礎(chǔ)上，通過內(nèi)蒙古草原區(qū)可燃物樣本采集進行燃燒實驗，獲得適用于草原區(qū)的王正非模型可燃物配置系數(shù)。使得王正非模型在模擬草原火蔓延的研究中有更高的精度和廣泛應(yīng)用性，由此彰顯模型的優(yōu)勢。玉山對內(nèi)蒙古草原70種主要草本可燃物燃燒特性進行測定，并分析草本間在各指標上的差異，采用主成分分析法對可燃物燃燒性的各項指標進行分析后對草本的燃燒性進行排序。利用聚類分析法對可燃物燃燒強度、速度、碳排放量、燃燒速率、產(chǎn)煙能力和點燃難易度等進行分析后，以內(nèi)蒙古草原高火險區(qū)草地類型為基礎(chǔ)將可燃物修正系數(shù)進行細分(見表3)。

表3 可燃物更正系數(shù)對應(yīng)表(Ks值)

3.2 草原火蔓延模擬預測

草原火蔓延模擬預測的實質(zhì)是草原火蔓延邊界，其計算最大的困難在于計算不同時刻草原火蔓延的速度?，F(xiàn)實情況下由于氣象、可燃物類型、地形等影響因子不同。在火蔓延的不同階段，即使在同一方向上蔓延速度也不盡相同。

3.2.1 草原火蔓延預測速度求解

草原火蔓延過程中火頭部分延伸和發(fā)展最快，所以利用火頭速度計算草原火蔓延過程中的最大速度。具體的，以Himawari-8數(shù)據(jù)為基礎(chǔ)，使用時序線性外推法求得的當前時刻的蔓延速度作為下一時段的蔓延初速度()；而下一時段的蔓延實際速度()則由、 可燃物配置格局修正系數(shù)、通過GRAPES＿MESO計算得到的預報風力修正系數(shù)、以及地形坡度修正系數(shù)相乘后得到。

欲更好的反映實際情況，將柵格空間近似為元胞空間，每個元胞即為代替的柵格。以元胞自動機為為基礎(chǔ)，將可燃物層、地形層、氣象層等因子賦值到元胞中，利用轉(zhuǎn)換規(guī)則和元胞狀態(tài)判斷火行為是否滿足持續(xù)發(fā)生條件。利用元胞自動機得到火行為發(fā)生結(jié)果與王正非模型和時序線性外推法結(jié)合起來，建立最優(yōu)速度求解模型。同時利用Hamawari-8靜止衛(wèi)星多時相的特點，采用動態(tài)預測模型，對火蔓延預測結(jié)果進行實時的更新和修正，以避免誤差積累，得到更為可靠的預測結(jié)果。具體來說，當新的火情監(jiān)測數(shù)據(jù)產(chǎn)生后，利用最近兩次(最新和上一監(jiān)測時刻)的火情數(shù)據(jù)更新各個方向上的蔓延速度，再次利用最優(yōu)速度模型對往后一定時間段內(nèi)的火蔓延速度進行預測，從而更新預測數(shù)據(jù)。

草原火從發(fā)生到發(fā)現(xiàn)的時間內(nèi)草原火已形成一定的規(guī)模，可燃物燃燒已達到一定的程度。利用Hamawari-8衛(wèi)星10min分辨率為觀測周期的特點，假設(shè)建立4個以10min為時間間隔的時間段，在第一個周期內(nèi)即0-10min未發(fā)現(xiàn)火情，在第二個周期10-20min內(nèi)發(fā)現(xiàn)火情則利用第三個周期20-30min利用時序線性外推法計算在本周期內(nèi)火頭蔓延的速度并且利用大興安嶺森林防火中心提供的火頭、火翼火尾蔓延速度關(guān)系經(jīng)驗資料(見表4)，預測各個方向的火蔓延速度。在Himawari-8衛(wèi)星未觀測到火情的時間段內(nèi)，火場已經(jīng)擁有一定的規(guī)模和強度，利用王正非模型計算該時間段內(nèi)的速度。所以算法獲取該時間段前一定時間的可燃物數(shù)據(jù)、氣象數(shù)據(jù)、地形數(shù)據(jù)計算衛(wèi)星未觀測到火情前一定時間內(nèi)火蔓延的速度。在第四周期30-40min 將線性時序外推法求得的當前時刻的蔓延速度作為下一時段的蔓延初速度(V)而下一時段的蔓延實際速度(V)則以V為基礎(chǔ)，進行各項影響因子的修正后得到最優(yōu)速度模型，式中參數(shù)與式(7)相同。

=***

(11)

表4 火頭、火翼、火尾蔓延速度關(guān)系表

3.2.2 基于最優(yōu)火蔓延速度的邊界預測

如圖所示，在火蔓延過程中，火頭通過元胞空間蔓延速度是最快的，而蔓延邊界上任意一點的蔓延速度與的比例關(guān)系()可近似為二者到初始火點的距離之比，即

=*()

(12)

圖2 蔓延邊界上任意一點F的蔓延速度與火頭H速度關(guān)系示意圖

據(jù)此，可以逐像元估算草原火沿任一方向蔓延到任一位置所需的時間。在給定的時間內(nèi)，則能估算任一方向蔓延所能到達的位置。

時序線性外推法要求計算相鄰兩次監(jiān)測時間段內(nèi)，各個方向上的火蔓延速度。因此，需要確定在每個監(jiān)測時刻各方向火場邊界所處的位置。首先，以初始火點為原點建立直角坐標系，為正東方向，為正北方向；然后利用初始火點和當前時刻其它所有火點的經(jīng)緯度，結(jié)合目標分辨率，將所有火點位置坐標化(，)。坐標化方法由式(13)給出。

(13)

式中，為火點的坐標位置；和分別為點的實際經(jīng)緯度；為目標分辨率0005°。

′=*+*

(14)

′=-*+*

(15)

={′|′=0′≥0}

(16)

圖3 通過坐標旋轉(zhuǎn)計算火場邊界位置

3.3 案例應(yīng)用與精度驗證

算法使用IDL語言編寫并實現(xiàn)可視化，針對實際火情，算法將產(chǎn)生以監(jiān)測到的初始火點時刻向后一定時間段內(nèi)(150min)的火情蔓延結(jié)果。蔓延結(jié)果每10分鐘更新一次，以500米空間分辨率的二值化圖像形式呈現(xiàn)。另外，算法將按時間順序疊合預測時間段內(nèi)的多個火情圖。算法的精度驗證是使用閾值法提取Himawari-8火點作為參考數(shù)據(jù)，與蔓延預測結(jié)果進行比較，從而計算出蔓延面積相對偏差、蔓延邊界均方差和像元平均絕對誤差3個精度指標對蔓延模型進行評價。

3.3.1 應(yīng)用案例

研究案例選取2016年3月29日位于內(nèi)蒙古錫林郭勒盟東烏珠穆沁旗薩麥蘇木發(fā)生的草原大火。薩麥蘇木位于內(nèi)蒙古高原東部，地勢總體上呈北高南低，由東向西傾斜，海拔在800-1956m之間，東烏珠穆沁旗全旗地形坡度在20°以下。地理位置介于北緯45°70′到46°58′，東經(jīng)116°25′到118°09′。薩麥蘇木春秋季多西北風，月平均風速3.6m/s，太陽輻射強烈，濕度低等因素致使烏珠穆沁草原火災頻繁并且薩麥草原北接蒙古國，蒙古國鄰接草原火災頻繁，過境大火經(jīng)常性侵襲，烏珠穆沁草原成為內(nèi)蒙古乃至全國草原火最為嚴重的區(qū)域之一，如圖4所示。

圖4 案例區(qū)示意圖

3.3.2 精度驗證

算法使用IDL編寫并實現(xiàn)可視化，對薩麥蘇木草原火做案例分析。算法從3月29日3時20分開始到5時50分共歷時150min。每10min更新一次模擬結(jié)果，蔓延結(jié)果見圖5。

圖5 模型模擬結(jié)果圖

根據(jù)模型模擬結(jié)果與實際過火數(shù)據(jù)做驗證分析，疊加Himawari-8第7、14波段解譯數(shù)據(jù)與模型模擬結(jié)果進行比較，在驗證分析中利用統(tǒng)計學方法對模型模擬結(jié)果與實際過火結(jié)果統(tǒng)計得到模型精度。得到蔓延范圍面積相對偏差Bias各時間點變化在1.49%-42.87%(見圖6)；蔓延邊界均方根誤差BRMSE在1.43-22.14個像元(見圖7)。

圖6 蔓延范圍面積相對偏差變化圖

圖7 蔓延邊界均方根誤差變化圖

如圖8(a)所示，模擬的整體形狀接近實際火場，在空間位置上與實際位置產(chǎn)生微偏差。圖8(b)中淺綠色部分為漏報區(qū)域、藍色部分為與實際重疊部分、黃色為誤報區(qū)域。像元平均絕對誤差CMAPE為18.09%；模型正確模擬區(qū)域為84.68%(見表6)。模型模擬結(jié)果整體與實際過火區(qū)域之間的偏移，是由于氣象數(shù)據(jù)不全且衛(wèi)星數(shù)據(jù)分辨率較低所引起。

圖8 模型精度對比圖

表5 蔓延精度像元數(shù)統(tǒng)計表

4 結(jié)論與展望

研究方法應(yīng)用于2016年3月29日3點20分位于內(nèi)蒙古錫林郭勒盟東烏珠穆沁旗薩麥蘇木發(fā)生的草原大火。分析了本文提出的方法在實際案例中的可用性，得出以下結(jié)果：

1)以Himawari-8靜止氣象衛(wèi)星結(jié)合時序線性外推法計算代替實驗獲取草原火蔓延初速度，并將其運用到王正非草原火模型中存在一定可行性。通過驗證得到蔓延范圍面積相對偏差Bias在1.49%-42.87%；蔓延邊界均方根誤差BRMSE在1.43-22.14個像元；像元平均絕對誤差CMAPE為18.09%；模型正確模擬區(qū)域為84.68%。

2)GRAPES＿MESO區(qū)域氣象數(shù)值預報產(chǎn)品可用于草原火蔓延模擬預測。彌補草原火蔓延模型在氣象數(shù)據(jù)缺失地區(qū)的所存在的不足。

3)本文在研究參閱大量相關(guān)國內(nèi)外文獻后，以薩麥蘇木草原大火為研究對象，以計算機技術(shù)為基礎(chǔ)，IDL系統(tǒng)平臺，進行草原火蔓延模擬可視化研究，得到了較好的可視化結(jié)果。

模擬草原火蔓延時著火區(qū)可燃物量的變化，風力的變化等都是面臨的問題。因而，建立火蔓延模型變的尤為重要。在建立草原火蔓延模型的過程中，模型的選擇，模型參數(shù)的細化都需要結(jié)合區(qū)域因素綜合考慮。同時結(jié)合多源、高時空分辨率的影像數(shù)據(jù)、增加更為詳細的氣象數(shù)據(jù)，能夠更好的為草原火蔓延預測研究提供數(shù)據(jù)基礎(chǔ)。