灾害规模、受灾程度,以及灾害的空间分布等信息是灾后应急指挥和决策的依据,特别是在发生严重地震灾害时,通讯与电力中断往往导致灾区与外界失去联系,如何快速获取灾情信息进行紧急救援成为指挥决策的迫切需求。长期以来,地震工作者一直探索着强震巨灾发生后如何快速有效地获取灾情信息。当前,国内外利用遥感技术获取地震灾害信息主要集中在对震后建筑物、桥梁、断层等地表形变破坏以及崩塌、滑坡等地震诱发次生灾害的研究(张景发等,2002; 单新建等,2005; 郭华东等,2011; 王晓青等,2015; 张小咏等,2016),但地震所带来的破坏还包括水、电、气、通讯、交通等生命线工程的损失(陆鸣等,2008; 刘爱文等,2012; 李东平等,2017)。一旦发生强震,灾区的电力系统可能受损,这种破坏不仅能反映地震造成的直接破坏,如房屋倒塌、道路中断等明显破坏,还能反映表面没有破坏而实际造成的经济损失,如电力中断导致的生产停止,这种软破坏对经济和民生也是一种直接的损失。而灯光遥感图像可以直接反映这类信息,在地震破坏范围确定、受灾人口、经济损失评估等有着重要的应用价值。

1976年,美国军事气象卫星(Defense Meteorological Satellite Program,简称DMSP)开始搭载的Operational Lines can System(简称OLS)传感器可以在夜间获取灯光影像,星下点空间分辨率1 km,具有每天覆盖全球的观测能力,幅宽达2 960 km,能够探测到城市灯光甚至小规模居民地、车流等发出的低强度灯光,并使之明显区别于黑暗背景(陈晋等,2003; 何春阳等,2006)。由于灯光直接与人类活动密切相关,DMSP/OLS已经广泛应用于多个方面,包括灯光与经济活动之间的相关性(Elvidge et al,1997),使用灯光数据定量估算美国、日本、中国,甚至全球人口(Sutton,1997; Kimberly,2010); 灯光与地震灾害之间的相关性,如日本和美国合作研制的地震灾害评估系统(Kohiyama et al,2004),且该评估系统已被应用到日本阪神地震评估中(Hayashi et al,2000)。

由于灯光遥感图像在地震灾害评估方面具有重要的应用价值,本文利用DMSP/OLS 数据对2001年印度7.9级地震和2011年“3·11”日本9.0级大地震后灯光遥感图像破坏信息进行提取,并讨论了灯光遥感数据用于强震巨灾后灾情信息获取的应用潜力和局限性。

北京时间2001年1月26日11时16分,印度古吉拉特邦发生7.9级地震,此次地震造成大量建筑物倒塌,2.5万人死亡,16.7万人受伤。2011年3月11日13时46分,日本东北部近海发生9.0级特大地震,并引发高达40.5 m海啸,造成15 985人死亡,13 346人失踪。这2次大地震发生后震区内电力、通信设备受到严重破坏。本文收集这2次地震前后灯光遥感图像用于分析,具体数据见图1、图2。从图中可以看出,夜间灯光遥感图像相对于白天的对地观测卫星图像非常简单,仅有灯光和背景2种信息,还可以清楚识别出城市和较大居民点,包括大城市和中小城市。对比印度和日本地震前后灯光图像(图1、图2),不同国家、不同经济发展水平、不同能源管理方式的灯光亮度差别很大,日本经济发达、灯光成片且亮度大,而印度灯光相对稀疏。同样,日本地震造成的灯光破坏也更加明显。

灯光代表了人类活动,正常情况下同一地区的灯光相对稳定(图2a,b,c)。地震发生后,震区灯光强度明显减弱,灯光范围减小或者消失,发生了显著变化(图2d),地震破坏引起灯光变化是震害信息的一个明显表征。基于上述现象分析,可以通过检测地震前后灯光遥感图像的变化提取震害信息。目前,基于像素的遥感图像

变化检测常用的方法有差值法和比值法,考虑到灯光遥感图像空间分辨率较低,且由地震引起的变化明显,本文选定差值法,该方法的关键在于如何确定差值的阈值,如果图像差值满足一定概率分布,基于数理统计的显著性检验方法是获得差值分割阈值的经典方法(李春干,代华兵,2017; 高波等,2012),也避免了人为给定一个阈值的随意性。

为了分析基于显著性检验方法进行变化检测的合理性,分别对2个研究区的不同时相图像进行了辐射归一化,然后计算了不同时相图像的差值,并根据差值直方图拟合得到高斯分布曲线(图3)。从图3可发现:(1)地震前不同时相图像差值完全满足正态分布;(2)图3a和3b中对称轴左侧的差值分布基本一样,而图3b中对称轴右侧明显不对称;(3)日本和印度地震前后图像差值规律一样(图3b,c),整体上符合正态分布,但差值在大于0方向的分布明显多于小于0方向,且在大于0方向有很多数值超过高斯分布曲线。

进一步分析可知,正常情况下灯光遥感图像没有明显变化,震前不同时相灯光遥感图像的差值主要是随机噪声。从统计意义上来讲,遥感图像每个像素都可以看作是一个统计样本,因此震前不同时相图像差值的直方图分布应该满足高斯分布,日本地震前图像差值的统计结果也证明了这一点(图3a)。当地震引起灯光变化后,地震

图1 2001年印度古吉拉特邦7.9级地震前后灯光遥感图像
Fig.1 Nighttime light images before and after the 2001 Gujarat M7.9 earthquake,India

图2 2011年日本9.0级地震前后灯光遥感图像
Fig.2 Nighttime light images before and after the 2011 Japan M9.0 earthquake

前后图像的差值直方图虽然整体上仍然符合正态分布规律,但是在大于0的方向很多差值分布在高斯分布曲线之外,这部分像素是由于地震引起灯光熄灭或者亮度减弱导致地震前后差值很大,从而出现图3b和3c中的非对称现象。日本地震的破坏范围大,对应的差值直方图在高斯拟合曲线之外的像素多; 而印度地震破坏相对小,在高斯分布曲线之外的像素相对较少,但表现的规律和特征都一样。因此,可以通过地震前后差值直方图拟合得到高斯概率密度分布函数,进而利用该函数的显著性检验的方法来提取地震导致的灯光破坏区域。

图3 不同时相灯光遥感图像差值直方图分布和高斯拟合曲线
Fig.3 The histogram and Gauss fitting curves of differences of nighttime light images

根据显著性检验,当概率分布函数和置信度确定后,可计算出概率密度函数的自变量值。本研究中,自变量是地震前后图像的差值,置信度对应的自变量数值是差值图像分割的阈值。置信度的设置通常采用0.9、0.95和0.99,0.99的置信度一般用于医学等要求非常高的领域,在遥感影像的变化检测中0.9的置信度结果存在较大的虚警率(Kohiyama et al,2004; 高波等,2012)。本文的灯光遥感图像分辨率较低,且由于多种原因使得图像存在较大的噪声,因此置信度设置为0.95,这样检测的结果不会导致虚警率或漏检率过高。由于变化区域仅仅体现在差值直方图大于0的方向,在利用显著性检验确定阈值时,选择右侧检验。

利用地震前后灯光遥感图像进行震害信息提取流程如图4所示。由于获取的DMSP/OLS灯光数据是DN值图像,没有进行地表辐射量的定量化转换和去云和去噪声等处理,所以先进行辐射归一化消除成像过程中产生的误差。具体采用统计量法,该方法假设同一地区的遥感图像具有相同的概率分布,即假设图像具有相同的均值和方差,所以辐射归一化可以通过源图像和参考图像概率分布的统计量进行匹配。根据这一假设则有:

((b₂-μ₂))/((b₁-μ₁))=(σ²)/(σ¹)(1)

对式(1)进行变换,得到统计量法的归一化模型为

b₁=(σ₁)/(σ₂)b₂-(σ₁)/(σ₂)μ₂+μ₁(2)

式中:b₁为原图像归一化后的像素灰度值; μ₁,σ₁分别为参考图像的像素均值和标准差; b₂,μ₂,σ₂分别为原图像的像素灰度值、均值和标准差。

辐射归一化后进行地震前后图像的差值计算,同时统计差值图像的直方图,并拟合得到高斯分布函数,高斯函数采用:

y=1/((2π)^1/2×a) e^{(-(x-b)2)/(a2)}(3)

利用高斯拟合的概率分布函数进行显著性检验,计算置信度为0.95条件下的差值图像的分割阈值,最后根据阈值对差值图像进行分割提取震害区域。

图4 基于灯光遥感图像的震害信息提取流程
Fig.4 Extraction flow chart of earthquake damage information based on nighttime light images

根据本文的计算方法和计算流程,首先对地震前后灯光遥感图像进行辐射归一化。在统计归一化转换参数时,为了降低云及其他噪声的影响,在图像中选择海洋上空若干不受云和噪声影响的区域作为统计样本,统计样本的均值与方差见表1,再根据式(2)进行辐射归一化,并计算归一化图像的差值,统计差值图像的归一化直方图,同时拟合得到高斯概率分布函数(图3),表1给出了拟合的高斯概率函数的参数。根据高斯拟合曲线计算置信度为0.95条件下的差值图像的分割阈值(表1),印度地震和日本地震分别对应差值为10和22。最后,根据阈值对差值图像进行阈值分割提取震害,差值大于阈值的像素作为提取的震害区域(图5)。

表1 印度和日本地震前后辐射归一化统计参数及高斯拟合参数和分割阈值
Tab.1 The parameters of radiometric normalization,Gauss fitting parameters and segmentation thresholds of the difference images before and after the India and Japan earthquakes

图5 基于灯光遥感影像的日本地震(a)和印度地震(b)灯光震害信息提取结果
Fig.5 Results of the lighting damages of Japan earthquake(a)and India earthquake(b)based on the nighttime light images

日本地震提取结果表明,破坏灯光主要分布在青森县、秋田县、岩手县、宫城县、茨城县(图5),其中宫城县和茨城县破坏灯光面积最大,岩手县和宫城县破坏灯光比例最大,达90%以上。图6a是基于灯光变化的震害破坏范围和地震动分布图。从图中可以看出,宫城县和茨城县中部地震动达到55～63 gal,对应位置的灯光也出现大范围破坏; 岩手县破坏灯光呈现带状分布,对应地震动的高值区,强烈的地震动导致大规模的电力设施破坏,因此,灯光破坏与地震动强度存在较好的一致性。本文提取的岩手县、宫城县的破坏灯光面积最大,与震前灯光遥感图像相比,破坏比例达到90%以上,这与日本气象厅发布的岩手县和宫城县受灾范围大,全县停电相一致; 茨城县破坏灯光较多,青森县、秋田县提取的破坏灯光相对较少,福岛县、山形县、新潟县只是出现少量灯光破坏,与日本气象厅发布发布的局部停电也大致一致。由于本文使用的3月12日DMSP/OLS灯光图像的空间采样为1 km,灯光遥感分辨率低使得这些县内的一些人口密度小的地区灯光无法识别,无法提取灯光破坏。对比图2a,b,c和图6b可以看出,正常情况下震前灯光遥感图像与人口密度空间分布具有非常高的相关性。对日本全年平均灯光强度与人口密度空间分布进行统计分析(以县为统计单元),见图7。从图7也可以看出,灯光强度与人口密度具有很好的相关性,因此,通过灯光破坏区域提取也能间接的获取地震受灾人口的规模和空间分布,这也是利用灯光遥感图像提取震灾害信息的优势。

图6 “3·11”日本地震灯光破坏、地震动和人口分布图
Fig.6 The lighting damage,ground motion and population distribution related to the “3·11” Japan earthquake

图7 日本年平均灯光强度与人口密度空间分布的散点图
Fig.7 Scatter diagram of the distribution of the average light intensity and population density in Japan

从“3·11”日本地震和印度地震应用效果来看,利用灯光遥感图像提取的破坏灯光能够从宏观上反映地震灾害的总体情况,特别是灯光破坏强弱及空间分布与受灾人口规模与空间分布也密切相关,对于震后第一时间的应急救援具有重要意义。但目前使用的DMSP/OLS灯光数据还存在一些缺陷,无论是空间分辨率还是辐射分辨率都比较低,不能满足对灯光分布和强度的精细化观测:(1)空间分辨率在星下点为1 km,边缘图像分辨率2.7 km,数据量化等级为6 bit,DN值0～64,这极大降低了数据对灯光的描述能力;(2)由于数据不能进行绝对辐射定标和大气校正等定量化处理,使得数据不能定量反演灾害程度;(3)灯光数据受云的影响大,遥感器波段设置应该考虑去云波段的设置。我国有关部门正在进行灯光遥感卫星的论证,应提高空间分辨率、辐射分辨率,同时提高数据定量化处理和去云处理的能力。

灯光是人类活动的标志,强烈地震发生后的灯光变化反映了电力供应与使用方面的变化,本质上是灾害造成的发电、输电等设施和线路的破坏。因此,灯光遥感影像可以作为地震破坏的一种显性载体。本文基于对灯光遥感影像和震后灯光变化原因分析,提供了一种利用灯光遥感影像进行地震灾害信息提取的方法,并在2001年印度地震、2011年“3·11”日本地震进行应用试验,发现灯光遥感图像目标单一,除了受云等噪声影响,主要包括灯光和背景信息,可以作为一种地震灾情快速获取的信息来源,利用灯光遥感影像不仅能够准确提取地震破坏范围和强度,而且能够直接反映了震后受灾人口规模与空间分布,对于震后应急救援具有重要意义。

2011年“3·11”日本大地震的地震动数据由日本消防厅Hasokawa博士提供,审稿专家提出了中肯的修改意见和建议,在此表示衷心感谢!