笔者自2013年以来开展了基于社会视频源的地震影响场勾画方法研究,基于多次震例探索了震后在海量的社会视频监控信息源支持下,利用定性分析烈度经验矩阵,建立监控视频与地震烈度的相对关系,在视频监控中通过模型提取震害和烈度信息,建立了基于视频运动估计的地面强震动定量分析方法,目前利用该方法可计算出监控点的近似地震动峰值加速度值。基于该原理开发了地震视频灾情判读系统(图1),并应用于2014年浙江文成—泰顺4.2级地震、2017年浙江磐安3.5级地震等应急处置工作(李东平等,2017)。

ASP是微软公司开发的代替CGI脚本程序的一种应用,它可以与数据库和其它程序进行交互,是一种简单、方便的编程工具。ASP是一种服务器端脚本编写环境,可以用来创建和运行动态网页或web应用程序,它包含HTML标记、普通文本、脚本命令以及COM组件等(胡秀源,2011),ASP的以上特点比较适合应急指挥中心的工作环境。地震视频灾情判读系统就是基于ASP开发的。

图1 地震视频灾情判读系统界面
Fig.1 Interface of earthquake video disaster assessment system

在近年来的几次地震应急中,笔者及团队不断改进视频分析方法,总结了在不同烈度下的探头晃动和人员反应表现形式,并推导出地震烈度经验矩阵,进一步进行以视频运动估计的精确定量计算。但该方法只对高烈度区、图像抖动明显的监控点识别较好,而低烈度区监控点的计算值与实际差距较大。在对算法进行了部分优化,引入场景深度估计方法后,图像帧之间像素的运动距离和速度值计算准确性均有明显提高。

地震视频分析的步骤如下:①对地震视频逐帧进行分析,利用宏块分析法计算出相邻帧之间发生位移的宏块。②对帧间图像发生位移的宏块进行统计,近而可以计算出地震视频中每一帧图像发生震动的像素个数。例如:导入的视频分辨率为320×384像素,每一帧图像共有20×24个16×16宏块。通过对地震发生时振动视频某帧图像进行运动计算分析,得出该帧图像的振动宏块为M块,该帧图像发生震动的像素为M×16×16个。③计算地震视频片段所有图像帧中上述运动像素的运动矢量均值,并将其作为该视频地震烈度的参考值(常安琪,2015)

由于不同的地震视频具有不同的地震加速度和速度值,为此需要构建地震视频数据库。在该数据库中,每个地震视频样本均以预先通过地震设备得到的Ground Truth地震数据(烈度、速度和加速度等)作为人工标注值。将具有相似烈度(比如某一范围内的烈度)的地震视频归为同一类,并通过支持向量机(Support Vector Machine,简称SVM)对不同烈度的地震视频样本进行学习,而输入的特征值即为该视频的运动矢量均值。利用训练好的SVM模型,对不同特征值输入的地震视频样本输出其烈度的预测。笔者及团队自2013年开始至今,已完成了1 000余个地震视频的建库工作。

通过地震视频计算地震的速度和加速度,需要知道地震视频中不同像素点之间代表的真实距离值。由于不同监控摄像机架设的高度和俯仰角都不同,因此像素之间的真实距离难以准确获取,主要是因为:不同监控场景中任意2个像素点的距离所对应的真实世界中的距离都是不同的,并且随着距离监控相机光心的距离变化而变化,无法通过参考图像中的某些物体尺寸做等价距离推算。一种近似的解决方案是通过机器学习技术估计场景深度信息,而这可通过单幅图像深度估计算法获得。

图2 场景深度估计算法
Fig.2 Algorithm for the depth of evaluation scenarios

场景深度估计算法步骤如图2所示。首先,对地震视频中的每幅图像进行无监督分割操作以生成超像素(Superpixel)区域。再以每个超像素为中心,尺寸为224×224像素的矩形图像块作为输入,训练一个具有5个卷积层和4个全连接层(5conv+4fc)的卷积神经网络(Convolutional Neural Network,简称CNN),如图3所示。在5个卷积层中分别使用了11×11,5×5和3×3的卷积模板,并采用修正线性单元(Rectified Linear Units,简称ReLU)作为激活函数以及2×2的池化操作; 在前2个全连接层中使用ReLU作为激活函数,第3个全连接层使用f(x)-(1+e^-z)-1作为激活函数,而最后一个全连接层则没有使用。最终的网络输出为对应每个超像素的一维深度值(杨宇翔等,2015)。

对该CNN模型进行训练时,样本采用的是NYU v2 Kinect 数据集和Make3D数据集,上述数据集不仅具有不同场景类别的图像样本,而且每个图像样本均有对应的采集自Kinect深度相机和激光雷达的真实距离数据。通过上述2个数据集训练得到的CNN模型,可输出每个超像素所对应的深度预测值Z_P。

图3 CNN模型结构
Fig.3 CNN model framework

然而,仅通过CNN模型推理得到的深度图通常存在相邻超像素的深度跳跃的情况。为了获取更为平滑而准确的深度图,笔者采用了条件随机场(Conditional Random Field,简称CRF)模型,通过最大后验概率(Maximum A Posteriori,简称MAP)估计最小化条件概率分布Pr(y|x):

Pr(y|x)=1/(z(x))exp[-E(y,x)](1)

式中:z(x)是z以观察序列x为条件的概率归一化因子。将CNN模型输出的每个超像素的深度预测值Z_P作为数据项U(y_p,x)(x为超像素特征,即超像素的深度值)的核心部分。在构建平滑项U(y_p,y_q,x)时,则通过一个具有全连接层的单层神经网络模型训练多种不同类型的超像素相似度R_pq,比如:相邻超像素p和q颜色均值之差,颜色直方图之差以及纹理特征值之差等。最后的能量函数E(y,x)定义为:

E(y,x)=∑_p∈NU(y_p,x)+∑_(p,q)∈δV(y_p,y_q,x)

=∑_p∈N(y_p-z_p)²+∑_(p,q)∈δR_pq(y_p-y_q)²(2)

最终经过CRF优化得到的场景深度估计结果如图4所示。在图4a和4c所对应的场景深度图4b和4d中,像素对比度越弱的代表距离相机越近,像素对比度越强则代表距离越远。

图4 基于CNN和CRF优化的场景深度估计
Fig.4 The depth of evaluation scenarios based on CNN and CRF optimized

在得到场景深度图之后,便可以估计所得到图像中像素之间所代表的大致真实距离。地震视频片段中每幅图像的地震速度和加速度均可由其运动矢量对应的真实运动距离分别计算得到。在地震视频片段中,第N帧图像中每个运动像素的速度V^N可由其运动矢量MV所对应的真实距离D^N(即运动矢量起点像素和终点像素所对应的深度值之差)除以时间T(即1/F,其中F为当前图像帧率)计算得到:

V^N=D^N_i/T(3)

地震视频的加速度a^N则是由相邻图像帧N和N+1之间的速度差除以时间T计算求得:

a^N=(V^N-V^N+1)/T(4)

由于该模块需要逐帧计算场景深度,因而运算速度较为缓慢。针对基于运动矢量统计截取的地震视频中可能存在一些误检,比如夜间车辆的灯光以及白天的车辆移动等场景时常会因为存在较大的像素运动而被误检为地震。为此,我们加入一些处理方法将上述可能的误检场景剔除。

2017年4月12日2时25分,浙江省杭州市临安区发生4.2级地震,微观震中位于临安区河桥镇与潜川镇交界处的九龙坑(30.08°N,119.34°E),震源深度15 km。震中附近区域震感强烈,部分房屋出现墙体裂缝、掉瓦等震害现象。地震发生后,浙江省地震局立即启动Ⅲ级应急响应。现场震害评估与科学考察组通过对25个调查点资料的综合整理分析,参考震源机制解和震中附近地震构造等资料,确定了此次地震的最高烈度为Ⅴ度,等震线呈椭圆形,长轴方向呈NE向,为7.4 km,短轴为3.8 km,V度区总面积约22.6 km²,主要涉及到临安市潜川镇和於潜镇。

地震发生后,浙江省地震应急指挥中心通过视频监控共享系统,在震后2 h内获取了震中附近发震时刻77个视频点位图像数据,通过对比摄像数据中人员反应、镜头晃动等情况,对烈度做了初步判定,并利用地震视频灾情判读系统对视频点的地面强震动进行了定量分析。震后4 h内得到了基于社会视频监控的烈度分布,计算得到的最高烈度与现场实际调查相同,等震线优势方向与现场工作结果一致,V度区面积和现场工作结果接近。

工作流程如下:首先,浙江省地震应急指挥中心启动了快速评估系统,绘制了基于经验模型的地震影响场分布图(图5),在模型计算得出的IV度区及其以上区域进一步选取了监控探头。浙江区域的烈度衰减模型采用经过修正的经典模型(李东平,2011):

I_a=3.701+1.019·M-0.93·ln(5+D)

I_b=3.469+0.968·M-0.846·ln(3+D)(5)

式中:D为震中距,M为震级。

图5 基于经验模型绘制的地震影响场分布图
Fig.5 Seismic influence field map based on empirical model

通过初步筛选,剔除了部分非高清视频和受干扰视频,选取了49个点位进行数据分析。个别摄像头处理结果和现场科考调查队得出的数据有一定差距,这是因为对比2013年以来的视频资源发现,每年的高清监控比例在不断增加,大量的高清探头替代了非高清探头(图6),获取的数据质量也大幅度提高。从第一时间获取的视频监控中,列出20个典型的摄像头,获取的资料包括标准的经纬度信息和视频分析得出的烈度加速度、现场调查烈度值(表1)。由于在Ⅴ度区范围内没有强震台,因此笔者采用的对比数据是现场调查数据。从表1可以看出视频分析加速度与现场调查烈度值有一定差距,但已经比较接近(图7)。

表1 视频分析加速度值与现场调查对比
Tab.1 The comparison of video analysis acceleration value and actual survey

图6 部分摄像头监控截图
Fig.6 Screenshots of part of social videos surveillance

图7 摄像头监控视频处理波形
Fig.7 Processing waveform of social video surveillance

第一时间获取的视频监控包含大量的非高清探头,丢帧现象普遍,在软件分析过程中,由于场景深度算法无法做到仪器分析精度,而且受场地条件的影响和安装附着物的影响,造成了分析结果与实际结果有一定的差距,分析处理结果仍不完善。在应急过程中,笔者及团队还获取了相当数量的互联网摄像头,但是这些摄像头的地理位置相对模糊,成像质量也不是很好,分析较为困难,所以在第一时间分析时并没有采用这些摄像头,只是在应急期结束后的总结过程中采用了部分这样的摄像头作为补充。

空间分析计算并未打包在软件中,由软件导出的数据在ArcGIS平台中完成,采用经验贝叶斯克里金法模拟震害分布。经验贝叶斯克里金法是一种地学统计插值方法,可自动执行构建有效克里金模型过程中的较复杂的步骤,主要通过估计基础半变异函数来解释所引入的误差。考虑到半变异函数估算的不确定性,该方法相对其他克里金方法所预测的标准误差更小,可准确预测一般程度上不稳定的数据,对于样本点不丰富的情况也更准确(李晓军等,2014)。在插值过程中,为避免由于样本点分布不均匀造成的插值算法失真,更好地反映震害分布,在视频点位相对少的地区人为添加了少量经验点来完善插值图形。插值后V度区结果较现场科考结果偏大,但较好地反映出了影响场的优势方向,也与现场科考结果接近。震中以北由于交通条件所限,现场调查选点较少,但视频监控不受交通限制,因此较好地弥补了现场工作缺陷(图8)。

图8 震后4小时视频烈度判读结果
Fig.8 The result of video intensity interpretation after 4 hours of earthquake