0　引 言

火山喷发会给全球的气候系统和生态环境带来灾难性和连锁性后果。2022年1月14—15日汤加火山大爆发，造成汤加国整体失联，并在环太平洋数十个国家地区引发海啸，其释放的能量远高于二战时期广岛原子弹的威力，危害波及半个地球^［1］。因此，对火山活动进行有效监测并提供及时预警，是一个亟待解决的重要问题^［2］。SO₂气体是火山活动的重要产物，会在岩浆的孕育及迁移过程中沿火山口逸出地表，并且其含量变化规律与岩浆上升动力学紧密相关，对于预测火山喷发的性质有重要的指示意义。因此，SO₂的排放速率是一个评估火山喷发风险的重要参数^［3］。

鉴于此，多种监测火山SO₂气体排放量的光学遥感技术相继问世，包括主动探测和被动探测两种类型。主动探测方式主要指激光雷达技术，如拉曼散射激光雷达或差分吸收激光雷达等，其主要优点是探测精度高，但存在时间分辨率低、空间分辨率低等缺点^［4］。被动探测的典型代表是近十几年兴起的SO₂紫外成像技术，其在时间分辨率和空间分辨率方面有显著优势，可以实时、快速并且精确地获取火山烟羽的二维图像，并直观地展现出视场中任意位置的SO₂浓度图像，能够对火山早期活动过程中烟羽动力学的空间演化进行实时分析。例如喷发前岩浆的浅层侵入，以及爆炸性喷发期间气体的快速上升和释放等，进而实现火山喷发的早期预警^［5］。

获得SO₂浓度图像以后，如何反演SO₂的排放速率，是SO₂紫外成像技术面向实际应用的一大难题，光流算法的提出提供了很好的解决思路。光流是指时变图像中对应点的亮度模式的运动速度，通过光流算法可以把空间物体的运动转化为传感器成像平面上的像素运动^［6］。光流算法获得的光流不仅包含图像目标的运动信息，还包含了整个图像丰富的空间三维结构的信息，因此被广泛应用于**航天^［7］、信息科学^［8］、交通监管^［9］、气象^［10］、环境监测^［11］、医学^［12］等多个领域。近年来，光流算法逐渐应用于火山SO₂排放速率的反演^［13-15］。

光流算法应用于气体运动信息的提取，尚存在诸多亟待解决的技术难题。光流算法基于图像连续帧之间的亮度和局部相位守恒的假设计算流速，因此光流在光照改变和遮挡条件下的鲁棒性以及光流的实时性运算严重阻碍了该技术的发展。并且在光流计算过程中需要使用光流约束方程对每个图像像素进行数值求解，然而光流约束方程面临重大挑战，因为其试图从一个单一的约束条件中找到垂直速度分量和水平速度分量，因此需要对流场施加空间上的一致性^［6］。此外，在实际的火山监测过程中，相机拍摄的整个烟羽的羽流可能存在湍流，从而导致光流算法的计算精度迅速降低。而且，在火山喷口源头附近，羽流扩散不明显导致的图像对比度问题同样会使得光流算法产生巨大误差^［16］。

鉴于此，从SO₂相机、光流算法以及神经网络的原理出发，提出了一种基于神经网络的光流算法，以识别和纠正潜在的非物理光流引起的误差，从而实现SO₂排放速率的精确反演。基于神经网络的光流算法的提出，是机器学习应用于传统光流算法的重要体现。其优点主要包括：1) 不受某种特定的光流算法限制；2) 具有强大的抗干扰能力；3) 可以实时反映羽流的空间变化；4) 可以解决光流算法边缘羽流截面排放速率反演的工程难题。

1　理 论

1.1　紫外相机

火山喷发时，排放的烟羽中不仅存在SO₂气体，还包含火山灰颗粒物，这些颗粒物的米氏散射和消光作用会对SO₂的探测造成显著的干扰。利用双通道SO₂紫外成像技术可以极大地消除颗粒物造成的影响，从而提高对SO₂浓度反演的准确性。SO₂相机由紫外灵敏相机和两个探测带宽为10 nm的紫外窄带滤光片组成，中心波长分别为310 nm和330 nm。其中，中心波长为310 nm的探测通道为SO₂吸收敏感通道，包含SO₂和气溶胶颗粒物，定义为A通道；另一个探测通道几乎不存在SO₂气体，只含有颗粒物，可用于颗粒物散射的校正，定义为B通道。

SO₂紫外相机的工作原理是朗伯比尔定律^［17］，其表达式为

式中I'和I0分别为穿过烟羽前后的光强，λ为波长，τ为光学厚度。

A通道受SO₂的吸收作用及颗粒物的消光作用的双重影响，其光学厚度τA的计算公式为

式中σ为SO₂分子的吸收截面，S为气体柱浓度，τm为颗粒物的光学厚度。

在B通道中，几乎不存在SO₂的吸收，只受颗粒物消光作用的影响，因此B通道的光学厚度τB只是颗粒物的光学厚度，即

联立式 (2) 和式 (3)，即可求得SO₂气体的光学厚度τSO2，即

式中IA'、IA0分别表示A通道太阳散射光穿过烟羽前后的光强，而IB'、IB0分别为B通道太阳散射光穿过烟羽前后的光强。

太阳光在不同太阳天顶角时穿过臭氧层的平均光程不同，导致到达近地面的紫外波段的太阳散射光的光谱强度差异很大。因此，SO₂紫外相机测得的光学厚度，尤其是τA，受太阳天顶角影响很大，需要对反演曲线实时定标，才能实现SO₂气体的准确探测。

差分吸收光谱 (DOAS) 定标法^［18］是借助光谱仪测量烟羽中某点位置的SO₂浓度时间序列，并将其与SO₂相机在相同空间位置测得的SO₂光学厚度进行匹配，采取函数拟合的方式获得定标曲线。由于该方法可减小颗粒物的影响，且能够校正光稀释效应的影响，因此逐渐成为SO₂紫外相机定标的主流方法。

利用DOAS法获得的定标曲线，可以将SO₂光学厚度图像转换为SO₂浓度图像，再利用光流算法，从浓度图像中提取烟羽的速度信息，浓度信息与速度信息相结合，即可得到SO₂气体的排放速率^［5］。横截面的SO₂排放速率Φ(𝓁) 的计算公式为

式中f为相机的焦距，i为羽流横截面k个样本数据中的一个，SSO2为SO₂浓度，veff为光流法获得的羽流横截面法线方向的速度分量，dpl为相机到烟羽的距离，Δs为羽流横截面的积分步长。

图1为SO₂相机拍摄烟羽的示意图。其中，相机拍摄的烟羽二维图像为图像A，取二维图像A中的直线NQ作为羽流横截面计算SO₂排放速率。在M点处，烟羽的实际速度为Vtrue，而在图像A上的速度为Veff。由于观测角度的原因，M点的烟羽速度相当于比实际的速度低估了一个余弦分子的乘积。相机与火山喷口距离超过10 km，M点可大致认为处于火山喷口的上方，M点的视线长度增加了一个余弦因子的乘积，大致抵消了速度被低估的那个余弦因子。因此，在选择适当观测角度的前提下，排放速率反演时可以忽略速度和距离从三维到二维转换的影响^［5］。

图 1. SO₂相机拍摄烟羽示意图

Fig. 1. Schematic diagram of smoke plume taken by SO₂ camera

下载图片 查看所有图片

1.2　神经网络和光流算法

SO₂气体排放速率的反演需要将紫外相机测量的气体浓度信息与羽流中气体的流速信息相乘，并对两者乘积沿羽流横截面积分［见式 (5)］。光流算法可以通过跟踪连续帧中不同的图像特征，在像素级监测气体的运动状态，从而计算得到整个烟羽的流速图像信息。因此，光流算法是紫外相机获取SO₂气体排放速率的关键技术。

1.2.1　光流算法

光流是空间运动的物体在成像平面上的像素级运动的瞬时速度，由场景中目标的运动、相机的移动或两者共同的运动产生，通常利用相邻帧之间的亮度恒定以及对应的取帧时间连续这两个假设来计算^［19-21］。在二维平面中，假设前一帧的时间为t，后一帧的时间为 t+dt。则前一帧位置I处的像素点的光强可以表示为I(x,y,t)，后一帧对应位置的光强为I(x+dx,y+dy,t+dt)。因此

将式 (6) 用泰勒级数展开，可以得到

式中o(x,y,t) 是泰勒级数展开式的高阶项，在“微小”运动的情况下可以忽略。因此

基于空间一致假设，可求得水平方向以及垂直方向的速度u、v，其计算公式为

式中∂I∂x、∂I∂y和It可由图像数据获得，而 (u,v) 即为所求光流。其中，u、v分别为

1.2.2　神经网络

神经网络是一种类比生物大脑中的神经网络信息传递并可以对传递的信息进行加工处理的数学运算模型，由数量众多的节点相互连接而成，由于其运转方式与人的神经网络信息传递过程类似，因此也被称为人工神经网络。本研究所使用的神经网络是反向传播 (BP) 神经网络，其模型包括输入层，隐含层和输出层，其中隐含层至少含有一个^［22］，其简单示意图如图2所示。

图 2. 神经网络简单示意图

Fig. 2. Simple diagram of neural network

下载图片 查看所有图片

假设f为激活函数，w_ij为i和j节点之间的权重，b_j为节点j的阈值，x_j为每个节点的输出值，则每个节点输入与输出的关系^{［23，24］}为

选择的激活函数f(x) 为

对于w_ij，有

同样对于b_j，有

在计算调整完隐含层以及输出层之间的权重值和输出层的阈值后，需要对输入层和隐含层的阈值进行调整计算。假设w_ki是输入层第k个节点和隐含层第i个节点之间的权重值，那么

再根据梯度下降法，对隐含层与输出层之间的权重值和阈值进行调整，可以得到

同样，对输入层和隐含层之间的权重值和阈值进行调整，可以得到

1.2.3　神经网络光流算法

神经网络光流算法是将BP神经网络融入光流算法的新型算法，可以优化光流算法求得的SO₂气体的流动速度，从而优化SO₂的排放速率的反演结果。该方法需要运用光流法求得成像平面上临近像素的光流速度，分别作为训练的输入输出数据，以获得理想的权重更新模型，并将该模型应用于整个图像来获得所需像素的光流流速。该方法旨在识别所有成功约束的运动矢量，从这些矢量中得出对羽流截面的速度矢量估计，并将之替换烟羽中的非物理运动矢量，从而实现对整个烟羽排放速率的精确反演。

2　排放速率反演

2.1　实验及仪器参数

火山SO₂气体排放速率监测实验的数据是Gliß等^［15］利用其自主研发的型号为Hamamatsu C8484-16C 的SO₂紫外相机，于2015年9月16日07：06至07：22，在距离Etna火山烟羽源约10.3 km的Milo镇采集，实验过程中主要仪器信息如表1所示。由于整个实验过程中用于排放速率反演所采集图片信息的时间较短 (约15 min)，太阳高度角变化并不明显，因此认为实验期间定标曲线无变化。

2.2　排放速率反演流程

图3为获取SO₂排放速率的流程图。利用SO₂相机的310 nm和330 nm两个通道分别采集天空背景以及火山烟羽图像，根据朗伯比尔定律得到两通道对应的光学厚度图像，而后对两光学厚度图像做差值计算获取SO₂光学厚度图像。采用光流法对连续的SO₂光学厚度图像进行处理可得气体的速度场分布信息，然后利用神经网络算法对速度场进行学习，即可获得神经网络光流算法的速度图像。此外，利用DOAS法获取的定标曲线将SO₂光学厚度图像转换为SO₂浓度图像，然后借助式 (5)，结合浓度图像和速度图像，即可反演得到SO₂的排放速率。

图 3. SO₂排放速率流程图

Fig. 3. Flow chart of SO₂ emission rate

下载图片 查看所有图片

2.3　浓度反演

SO₂紫外相机A、B两通道采集的天空背景及火山烟羽图像经电子偏移以及暗电流校正处理，并通过朗伯比尔定律计算后，得到的光学厚度图像分别如图4 (a)、(b) 所示。结合图4 (c) 所示的DOAS方法获得的定标曲线，做两通道光学厚度图像的差值图像，即SO₂气体的光学厚度图像，进行线性转换，即可得到SO₂气体的浓度图像，如图4 (d) 所示。

图 4. SO₂浓度反演。(a) 310 nm光学厚度图像; (b) 330 nm光学厚度图像; (c) 定标曲线; (d) SO₂浓度图像

Fig. 4. SO₂ concentration inversion. (a) 310 nm optical density image; (b) 330 nm optical density image;(c) calibration curve; (d) SO₂ concentration image

下载图片 查看所有图片

2.4　速度获取

SO₂相机具有较高的空间分辨率和时间分辨率，能够实现火山烟羽SO₂浓度图像的实时获取。利用SO₂浓度图像反演得到SO₂排放速率，是火山监测领域极为重要的工程问题。传统的速度估计方法主要有三种：风速测量、交叉相关法以及光流算法。其中风速测量获得的速度是烟羽周围大气运动的速度，无法体现烟羽运动的真实情况；交叉相关法无法解析气体烟羽内部复杂的运动状态，会使得反演结果不准确；光流算法在运算速度以及计算精度等方面都具有显著的优势，是估算火山SO₂气体扩散运动速度场的主流算法，本研究所使用的是在火山排放速率反演方面具体良好效果的Farneback光流算法^［14］。

Farneback光流算法的工作原理是将图像中每个像素的局部邻域表示为多项式进行展开，然后计算每个连续帧的多项式之间的位移，并通过强制执行某个运动模型 (例如仿射运动模型)，估计二维位移场^［15］。此外，Farneback算法还采用了迭代方法，从粗特征开始估计运动，然后使用逐步细化的特征进行细化，因而具有较高的计算精度。

Farneback光流算法在理论中被证明有良好的准确性，但在实际的烟羽气体运动评估中，光流计算结果容易受图像对比度和湍流影响^［16］。为抑制这两个因素的影响，本研究提出一种基于神经网络的光流算法，对Farneback算法计算产生的密集位移矢量场进行分析训练，推导出一个经过神经网络迭代的速度矢量，用来取代位移矢量场中的速度矢量。

为消除异常值的影响，在设计神经网络训练模型之前需添加一个3 pixel× 3 pixel的中值滤波器对SO₂光学厚度图像进行滤波。在训练模型时，每一次训练的输入输出数据为当前帧图像中的相隔20列的羽流截面的气体光流数据，从而保证训练数据的实时性。图5显示了07：12第55 s的SO₂浓度图像，包括浓度图像上的光流及图像坐标，烟羽中竖直的线条表示接下来用于排放速率反演的横截面。每个位移矢量场包含10 pixel× 10 pixel。为突出羽流运动细节，每个显示的矢量长度被扩大了3倍。

图 5. 原始光流图像 (包括用于排放速率反演的羽流横截面)

Fig. 5. Original optical flow image (including plume cross-section for emission rate inversion)

下载图片 查看所有图片

图6为采用交叉相关法、Farneback光流算法以及神经网络光流算法获得的速度曲线图，其流速方向为羽流截面对应的法线方向。可以看出，与交叉相关法相比，Farneback光流算法以及神经网络光流算法更能反映出SO₂气体运动速度的瞬时变化。神经网络光流算法与Farneback光流算法得到的SO₂气体速度变化在时间序列上吻合性很好，证明神经网络光流算法继承了Farneback光流算法能够解析气体烟羽内部复杂运动状态的优点，同时，神经网络光流算法因其优化解的优点极大地削弱了湍流和图像低对比度对反演结果的异常影响。另外，神经网络光流法 (标准差为1.24) 的数据稳定性优于Farneback光流算法 (标准差为2.23) 的数据稳定性，显示了神经网络光流算法对湍流以及图像低对比度的强大的抗干扰能力。

图 6. 交叉相关法、Farneback光流算法以及神经网络光流法的速度

Fig. 6. Speed of cross correlation method, Farneback optical flow and neural network optical flow method

下载图片 查看所有图片

2.5　排放速率获取

为进一步验证模型的准确性并且对比交叉相关法、Farneback光流算法以及神经网络光流法的优劣，对三种方法反演的SO₂排放速率进行了对比，如图7所示。交叉相关法反演结果虽然符合火山烟羽排放的周期性，但对SO₂排放速率的瞬时变化不敏感。在无规律湍流以及图像低对比度的双重影响下，Farneback光流算法反演得到的SO₂排放速率存在异常的峰值，而神经网络光流算法明显削弱了这两个因素的影响，使得反演结果更加符合火山烟羽排放的周期性。此外，神经网络光流算法反演的SO₂排放速率与Gliß等^［15］的反演结果的最大值、最小值以及均值的比值都接近1，进一步证明了神经网络光流算法反演SO₂排放速率的准确性。

图 7. 交叉相关法、Farneback光流算法以及神经网络光流法排放速率反演

Fig. 7. Emission rate retrieval of cross correlation method, Farneback optical flow and neural network optical flow method

下载图片 查看所有图片

在靠近图像边缘的像素由于不符合“微小”运动假设，因此采用Farneback光流算法反演图像边缘处的SO₂排放速率时，会产生较大误差，而神经网络光流法可以很好地解决此问题。图8为图像边缘处 (第1343列) 羽流截面的排放速率反演结果。由图8与图7对比可以看出，在边缘效应影响下，Farneback光流算法计算得到的排放速率已经完全偏离实际值，而神经网络光流算法依然表现出较高的计算精度。

图 8. 边缘羽流横截面SO₂排放速率

Fig. 8. SO₂ emission rate of edge plume cross-section

下载图片 查看所有图片

为更详细客观地展现出交叉相关法、Farneback光流算法与神经网络光流算法三种方法在反演边缘排放速率的优劣，对1243～1343列的SO₂平均排放速率进行对比，如图9所示。由图可以明显看出，交叉相关算法得到的反演结果比另外两种方法在1243～1303列的正常区域的反演结果高1.5～1.6倍，其误差稳定在50%～60%之间。而Farneback光流算法与神经网络光流算法在远离图像边缘 (1243～1315列) 时，两者具有较好的一致性，进一步证明神经网络光流算法的计算精度较高。然而，越靠近图像边缘 (即在第1315列之后)，利用Farneback光流算法反演的SO₂平均排放速率越小，误差越大，在最边缘处，误差高达94%，这是由于在图像边缘，用于反演的上一帧图像像素的运动信息与当前帧图像像素对应的运动信息难以匹配，从而造成光流反演的巨大误差。但在图像低对比度及烟羽湍流效应的双重影响下，神经网络光流算法依然保持较高的准确度，其反演误差不超过5%，并且其边缘5列的平均误差不超过10%。由此可见，神经网络光流算法能够实现边缘羽流截面排放速率的精确反演。

图 9. 第1243～1343列羽流横截面SO₂平均排放速率 (a) 及其误差 (b)

Fig. 9. Mean SO₂ emission rates (a) and errors (b) for plume cross-sections in columns 1243-1343

下载图片 查看所有图片

3　结 论

介绍了SO₂相机的物理机理以及光流法与神经网络算法的数学原理，详细描述了SO₂浓度图像的获取方法，细致地阐述了利用Farneback算法及神经网络光流算法计算SO₂气体排放速率的实现过程。以Gliß等^［15］采集的Etna火山烟羽图像为处理对象，分别用交叉相关法、Farneback光流算法以及神经网络光流算法反演得到了火山SO₂的排放速率。实验结果对比显示了神经网络光流算法的科学性，以及可清晰反映羽流实时空间变化的优势，表明神经网络光流算法可以约束不规则物理运动，能够有效减小湍流和图像低对比度影响的优点。最后，对边缘羽流截面的排放速率进行计算，结果显示，在靠近图像边缘位置处，交叉相关法获得的排放速率的误差处于50%～60%之间，Farneback光流算法获得的排放速率误差超过90%，而神经网络光流算法的误差只有5%，并且神经网络光流算法边缘5列的平均误差不超过10%，因而证明神经网络光流算法具有可实现边缘羽流截面SO₂排放速率精确反演的显著优势。

与传统光流算法相比，神经网络光流算法能够有效抑制烟羽湍流、图像低对比度，以及图像边缘效应的影响，且具有较高的计算精度和工程可靠性。SO₂紫外相机本身具有高时间分辨与高空间分辨率等优势，而神经网络光流算法具有像素级精确检测物体运动信息以及高速准确优化求解的技术优点，将神经网络光流算法应用于SO₂紫外相机技术可实现对羽流排放速率的实时、快速、精确的反演。本工作将推动SO₂紫外相机在工业烟囱及船舶尾气遥感监测中的应用。