刘 洋,玉 山,2,都瓦拉,杨晓颖
(1. 内蒙古师范大学地理科学学院,内蒙古 呼和浩特市 010022; 2. 内蒙古自治区遥感与地理信息系统重点实验室,内蒙古 呼和浩特市 010022; 3. 中国农业科学院草原研究所,内蒙古 呼和浩特市 010022)
我国草原饱受火灾困扰,草原火灾突发性强、蔓延速度快、危害性大。造成经济损失和人员伤亡的同时毁坏草场资源,生物多样性降低,破坏草地生态系统的平衡。为了及时、准确地防控草原火灾,有必要对草原火蔓延进行研究。
新一代静止气象卫星时空分辨率得到大幅提高,其高频次、覆盖范围广的特点在监测草原火灾方面有较大的优势。目前,国内外学者利用Himawari-8影像监测草原火。在此基础上,一些学者利用Himawari-8监测成果结合王正非火蔓延模型与元胞自动机等计算机仿真算法,以栅格或矢量图形来预测草原火蔓延,并取得了较好的结果。
GRAPES_MESO是GRAPES(Global/Regional Assimila-tion and Prediction System)系统中区域中尺度数值预报模式,其在天气预报业务和防灾减灾方面发挥着重要的作用。国内学者利用其高时空分辨率的特点在气象预警、模型构建等方面进行研究与应用。
综合前人研究经验,本文以内蒙古草原高火险地区为研究区域,首先利用Himawari-8与GRAPES_MESO气象预报数据特征,由时序线性外推法确定草原火蔓延初速度,然后耦合元胞自动机和王正非火蔓延模型对草原火蔓延速度和蔓延边界进行预测,并应用到实际案例验证预测结果。
2.1.1 Himawari-8数据介绍
Himawari-8是由日本气象厅发射的新一代静止气象卫星,该卫星对地观测频率高达 10 min/次,空间分辨率最高达500m。其搭载的AHI(Advanced Himawari Imager)传感器波段范围覆盖可见光到远红外,设有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)。
本文选取内蒙古锡林郭勒盟东乌珠穆沁旗萨麦苏木2016年3月29日2时-7时分辨率为2km的影像数据。
2.1.2 GRAPES_MESO简介
GRAPES_MESO中国及周边区域数值预报产品,其空间分辨率为10km,时间分辨率为3小时。预报时次为00时、12时,预报最高时效72小时,要素包括位势高度、温度、风的u向量、风的v向量、相对湿度、降水量等(http:∥data.cma.cn)。
本文利用风的u和v向量计算2016年3月29日00时次6小时预报时效的10m特定高度的风速和风向。
2.1.3 辅助数据
本文模型和算法主要基于栅格数据处理计算,模型的主要输入参数来源于30m、500m和2000m空间分辨率的遥感产品和气象数值预报产品。遥感数据获取于美国航天航空局(National Aeronautics and Space Administration,NASA)的MODIS 500米分辨率地表覆盖数据(https:∥ladsweb.modaps.eosdis.nasa.gov/search/order);草地类型数据由中国农科院草原研究所提供;DEM数据收集自美国航天飞机雷达地形测绘使命(Shuttle Radar Topography Mission,SRTM)30米分辨率1_Arc_Second SRTMDEM数据(https:∥earthexplorer.usgs.gov/);可燃物修正系数参照文献[10]中可燃物燃烧实验结果获取。
Hamawari-8卫星对地观测频率高达 10 min/次,基于连续两次观测反演的火情得到各个方向上火场边界的扩展距离,进而计算该间隔内蔓延边界上草原火蔓延初速度。假设蔓延初始时刻(设=0)、蔓延开始后时刻、蔓延开始后时刻的火情如图1所示(坐标原点为初始火点),则在某一方向上,到时间段内火场边界的蔓延速度则为
(1)
其中,为时刻方向上火场边界与初始火点的距离。同理,在该方向上到时间段内火场边界的蔓延速度为:
(2)
一般地,在方向上,-1到时间段内火场边界的蔓延速度为:
(3)
图1 连续三个观测时刻的火场示意图
得到蔓延速度后,即可预测当前时刻往后一定时间段内的火蔓延情况。对于方向,+1时刻的火场边界位置为
,+1=,*(+1-)
(4)
地理元胞自动机(Geograph Cellular Automaton,GeoCA)具有许多特性,使其在地理时空建模中具有吸引力。利用简单的局部规则,它们有可以模拟复杂的现象。此外,它们本质上具有空间性,其结构与各种数字源提供的地理空间数据集兼容。标准元胞自动机是一个四元组即
=(,,,)
(5)
式中:为正整数,为元胞的基本单元;为元胞自动机中元胞的状态,每个元胞状态是一个有限或无限的元素;为元胞自动机的邻域,邻域的状态是决定该元胞是否进行转换,邻域的值通过土地覆盖数据赋予;为元胞自动机空间中的转换规则。当元胞自动机在二维空间上时,元胞空间结构与栅格空间结构极为相似,故以500m分辨率元胞自动机来模拟草原火蔓延。
根据草原火蔓延过程中各个阶段,每个元胞的状态用连续的整数值来表示():=0:“不可燃烧”状态是指此类元胞的属性本身不具有可燃性,不具备向相邻元胞扩散的能力。如沙地、水体、公路等;=1“未燃烧”状态是指该元胞具有可燃性,由于草原火未将该元胞所表示单元点燃;=2:“燃烧”状态是指元胞具备可燃性并且可以向其相邻的扩散;=3:“熄灭”状态是指该元胞具有可燃性并且该元胞已经燃烧完毕。
草原火蔓延空间可以解释为二维元胞空间,针对草原火蔓延提出以下元胞转换规则:1:如果元胞(,)在时刻的状态为2,该元胞相邻元胞状态不存在0或3 则其状态根据式(6)计算。
2:如果元胞(,)在时刻的状态为2,在+1时刻元胞(,)状态转变为3,则表示该元胞烧毁,元胞状态值保持在3不变,并且不会向相邻元胞扩散。
3:如果元胞(,)在时刻的状态为2,但是相邻元胞状态值为0或3,则元胞(,)不会传递。
在+1时刻,元胞(,)状态利用规则对相邻元胞进行扩散,即
(6)
草原火蔓延过程中不同路径上可燃物、气象、地形要素的影响对草原火蔓延速度产生的差异,毛贤敏对王正非模型中风和地形因子的细化,结合文献[10]可燃物修正系数的研究成果对模型进行修正。王正非火蔓延模型涉及四个参数:初始蔓延速度、可燃物配置格局修正系数、风力修正系数、地形坡度修正系数。
=
(7)
式中:为火蔓延速度m/min;为初始蔓延速度m/min;是可燃物配置格局修正系数(无因次);为风力修正系数(无因次);(1cos)为地形坡度修正系数(无因次)。
251 蔓延范围面积相对偏差Bias
从监测到初始火点开始的一定时间段内,火蔓延范围的面积可由蔓延边界包含在内的像元面积相加而得。若模型模拟得到的火蔓延面积为S,参考数据得到的该面积为S′,则蔓延范围面积的相对误差Bias由式(8)计算得到。Bias越小,蔓延预测结果越精确,反之亦然。
Bias=(-′)′
(8)
2.5.2 蔓延边界均方根误差BRMSE(Boundary Root Mean Square Error)
(9)
式中BRMSE单位为像元(Pixel),Res为像元分辨率。
2.5.3 像元平均绝对误差CMAPE(Cell Mean Absolute Percentage Error)
模型模拟的像元总数为N,模型模拟区域与实际火情区域相交区域像元数为N,CMAPE为N与N的差的比值。CMAPE0%表示完美模型,CMAPE大于 100 %则表示劣质模型。
(10)
草原火蔓延过程较为复杂,在算法实现过程中输入数据的空间分辨率不同,草原火蔓延预测计算统一采用500m像元分辨率和WGS_1984_Albers投影系统。
首先,使用ArcGIS提取DEM数据中包含的坡度数据、利用IDL平台将GRAPES_MESO中国及周边区域3小时数值预报产品中表示经度方向上的风向量u和纬度方向上的风向量v计算得出风速和风向数据,将提取的坡度数据和风速数据分级赋值(见表1和表2)。将Himawari-8数据经过几何校正、重投影、图像配准得到更为精准的分析数据,发生火灾时中红外第7(3.74-3.96μm)波段在时相上的亮温变化明显,亮温会迅速升高,而第14(11.0-11.30μm)通道亮温变化较小,产生背景亮温差。因此根据时相变化设置阈值检测火点,通过运算获得影像过火区以及火点,作为草原火蔓延的现势火情数据。
保证草原火蔓延模型与辅助数据之间的兼容性,算法统一采用IDL语言编写,包括草原火蔓延主程序和地形数据、气象数据标准化函数。程序首先读取Himawari-8数据,之后进行数据标准化和分级赋值处理。地形数据根据在相关研究中将坡度分划为6个方向的修正公式,根据研究区范围内DEM数据,通过修正公式计算获得分级数值(见表1),风力等级标准根据国家风力等级划分标准赋值(见表2),可燃物数据与草地类型数据一一对应赋值分级。
表1 地形坡度修正系数(值)
表2 风力修正系数(值)
我国北方草原当年10月至来年4月容易发生火灾。在此基础上,通过内蒙古草原区可燃物样本采集进行燃烧实验,获得适用于草原区的王正非模型可燃物配置系数。使得王正非模型在模拟草原火蔓延的研究中有更高的精度和广泛应用性,由此彰显模型的优势。玉山对内蒙古草原70种主要草本可燃物燃烧特性进行测定,并分析草本间在各指标上的差异,采用主成分分析法对可燃物燃烧性的各项指标进行分析后对草本的燃烧性进行排序。利用聚类分析法对可燃物燃烧强度、速度、碳排放量、燃烧速率、产烟能力和点燃难易度等进行分析后,以内蒙古草原高火险区草地类型为基础将可燃物修正系数进行细分(见表3)。
表3 可燃物更正系数对应表(Ks值)
草原火蔓延模拟预测的实质是草原火蔓延边界,其计算最大的困难在于计算不同时刻草原火蔓延的速度。现实情况下由于气象、可燃物类型、地形等影响因子不同。在火蔓延的不同阶段,即使在同一方向上蔓延速度也不尽相同。
3.2.1 草原火蔓延预测速度求解
草原火蔓延过程中火头部分延伸和发展最快,所以利用火头速度计算草原火蔓延过程中的最大速度。具体的,以Himawari-8数据为基础,使用时序线性外推法求得的当前时刻的蔓延速度作为下一时段的蔓延初速度();而下一时段的蔓延实际速度()则由、 可燃物配置格局修正系数、通过GRAPES_MESO计算得到的预报风力修正系数、以及地形坡度修正系数相乘后得到。
欲更好的反映实际情况,将栅格空间近似为元胞空间,每个元胞即为代替的栅格。以元胞自动机为为基础,将可燃物层、地形层、气象层等因子赋值到元胞中,利用转换规则和元胞状态判断火行为是否满足持续发生条件。利用元胞自动机得到火行为发生结果与王正非模型和时序线性外推法结合起来,建立最优速度求解模型。同时利用Hamawari-8静止卫星多时相的特点,采用动态预测模型,对火蔓延预测结果进行实时的更新和修正,以避免误差积累,得到更为可靠的预测结果。具体来说,当新的火情监测数据产生后,利用最近两次(最新和上一监测时刻)的火情数据更新各个方向上的蔓延速度,再次利用最优速度模型对往后一定时间段内的火蔓延速度进行预测,从而更新预测数据。
草原火从发生到发现的时间内草原火已形成一定的规模,可燃物燃烧已达到一定的程度。利用Hamawari-8卫星10min分辨率为观测周期的特点,假设建立4个以10min为时间间隔的时间段,在第一个周期内即0-10min未发现火情,在第二个周期10-20min内发现火情则利用第三个周期20-30min利用时序线性外推法计算在本周期内火头蔓延的速度并且利用大兴安岭森林防火中心提供的火头、火翼火尾蔓延速度关系经验资料(见表4),预测各个方向的火蔓延速度。在Himawari-8卫星未观测到火情的时间段内,火场已经拥有一定的规模和强度,利用王正非模型计算该时间段内的速度。所以算法获取该时间段前一定时间的可燃物数据、气象数据、地形数据计算卫星未观测到火情前一定时间内火蔓延的速度。在第四周期30-40min 将线性时序外推法求得的当前时刻的蔓延速度作为下一时段的蔓延初速度(V)而下一时段的蔓延实际速度(V)则以V为基础,进行各项影响因子的修正后得到最优速度模型,式中参数与式(7)相同。
=***
(11)
表4 火头、火翼、火尾蔓延速度关系表
3.2.2 基于最优火蔓延速度的边界预测
如图所示,在火蔓延过程中,火头通过元胞空间蔓延速度是最快的,而蔓延边界上任意一点的蔓延速度与的比例关系()可近似为二者到初始火点的距离之比,即
=*()
(12)
图2 蔓延边界上任意一点F的蔓延速度与火头H速度关系示意图
据此,可以逐像元估算草原火沿任一方向蔓延到任一位置所需的时间。在给定的时间内,则能估算任一方向蔓延所能到达的位置。
时序线性外推法要求计算相邻两次监测时间段内,各个方向上的火蔓延速度。因此,需要确定在每个监测时刻各方向火场边界所处的位置。首先,以初始火点为原点建立直角坐标系,为正东方向,为正北方向;然后利用初始火点和当前时刻其它所有火点的经纬度,结合目标分辨率,将所有火点位置坐标化(,)。坐标化方法由式(13)给出。
(13)
式中,为火点的坐标位置;和分别为点的实际经纬度;为目标分辨率0005°。
′=*+*
(14)
′=-*+*
(15)
={′|′=0′≥0}
(16)
图3 通过坐标旋转计算火场边界位置
算法使用IDL语言编写并实现可视化,针对实际火情,算法将产生以监测到的初始火点时刻向后一定时间段内(150min)的火情蔓延结果。蔓延结果每10分钟更新一次,以500米空间分辨率的二值化图像形式呈现。另外,算法将按时间顺序叠合预测时间段内的多个火情图。算法的精度验证是使用阈值法提取Himawari-8火点作为参考数据,与蔓延预测结果进行比较,从而计算出蔓延面积相对偏差、蔓延边界均方差和像元平均绝对误差3个精度指标对蔓延模型进行评价。
3.3.1 应用案例
研究案例选取2016年3月29日位于内蒙古锡林郭勒盟东乌珠穆沁旗萨麦苏木发生的草原大火。萨麦苏木位于内蒙古高原东部,地势总体上呈北高南低,由东向西倾斜,海拔在800-1956m之间,东乌珠穆沁旗全旗地形坡度在20°以下。地理位置介于北纬45°70′到46°58′,东经116°25′到118°09′。萨麦苏木春秋季多西北风,月平均风速3.6m/s,太阳辐射强烈,湿度低等因素致使乌珠穆沁草原火灾频繁并且萨麦草原北接蒙古国,蒙古国邻接草原火灾频繁,过境大火经常性侵袭,乌珠穆沁草原成为内蒙古乃至全国草原火最为严重的区域之一,如图4所示。
图4 案例区示意图
3.3.2 精度验证
算法使用IDL编写并实现可视化,对萨麦苏木草原火做案例分析。算法从3月29日3时20分开始到5时50分共历时150min。每10min更新一次模拟结果,蔓延结果见图5。
图5 模型模拟结果图
根据模型模拟结果与实际过火数据做验证分析,叠加Himawari-8第7、14波段解译数据与模型模拟结果进行比较,在验证分析中利用统计学方法对模型模拟结果与实际过火结果统计得到模型精度。得到蔓延范围面积相对偏差Bias各时间点变化在1.49%-42.87%(见图6);蔓延边界均方根误差BRMSE在1.43-22.14个像元(见图7)。
图6 蔓延范围面积相对偏差变化图
图7 蔓延边界均方根误差变化图
如图8(a)所示,模拟的整体形状接近实际火场,在空间位置上与实际位置产生微偏差。图8(b)中浅绿色部分为漏报区域、蓝色部分为与实际重叠部分、黄色为误报区域。像元平均绝对误差CMAPE为18.09%;模型正确模拟区域为84.68%(见表6)。模型模拟结果整体与实际过火区域之间的偏移,是由于气象数据不全且卫星数据分辨率较低所引起。
图8 模型精度对比图
表5 蔓延精度像元数统计表
研究方法应用于2016年3月29日3点20分位于内蒙古锡林郭勒盟东乌珠穆沁旗萨麦苏木发生的草原大火。分析了本文提出的方法在实际案例中的可用性,得出以下结果:
1)以Himawari-8静止气象卫星结合时序线性外推法计算代替实验获取草原火蔓延初速度,并将其运用到王正非草原火模型中存在一定可行性。通过验证得到蔓延范围面积相对偏差Bias在1.49%-42.87%;蔓延边界均方根误差BRMSE在1.43-22.14个像元;像元平均绝对误差CMAPE为18.09%;模型正确模拟区域为84.68%。
2)GRAPES_MESO区域气象数值预报产品可用于草原火蔓延模拟预测。弥补草原火蔓延模型在气象数据缺失地区的所存在的不足。
3)本文在研究参阅大量相关国内外文献后,以萨麦苏木草原大火为研究对象,以计算机技术为基础,IDL系统平台,进行草原火蔓延模拟可视化研究,得到了较好的可视化结果。
模拟草原火蔓延时着火区可燃物量的变化,风力的变化等都是面临的问题。因而,建立火蔓延模型变的尤为重要。在建立草原火蔓延模型的过程中,模型的选择,模型参数的细化都需要结合区域因素综合考虑。同时结合多源、高时空分辨率的影像数据、增加更为详细的气象数据,能够更好的为草原火蔓延预测研究提供数据基础。