2. 广东省粤电集团有限公司 天生桥一级水电开发有限责任公司水力发电厂, 贵州 兴义 562400
2. Tianshengqiao-I Hydro-Power Development Limited Company Hydropower Plant, Guangdong Yudean Group Company Limited, Xingyi Guizhou 562400, China
通过放大视频中的微小运动,可以展现一些肉眼不可见的运动信息[1],在可视化、临床诊断、监控等领域[2-4]都有着广阔的应用。2005年Liu等[5]提出了一种针对视频的运动放大技术 (Motion Magnification,MM),该方法对目标的特征点进行聚类、跟踪,然后再放大运动幅度,但该方法计算复杂度高,并且相关参数难以控制,在实际应用中采用得不多。2012年Wu等[6]提出了一种称为欧拉视频放大 (Eulerian Video Magnification, EVM) 方法,该方法不是显式地跟踪和估计粒子的运动,而是通过分析整个场景图像的像素点的值随时间的变化,从中分离出感兴趣的频段并进行增强,从而达到放大动作信息的目的。2013年Wadhwa等[7]提出了一种基于相位的运动放大技术 (Phase-Based Video Motion processing,PBVM),该方法通过放大相位来改变运动目标的位移变化,能够有效抑制噪声,支持更大的放大倍数。但是,EVM和PBVM方法均需要用户为每个待放大的视频提供一组输入参数才能进行放大,对于未知的视频,这需要耗费时间去找到合适的参数,对于非专业的用户将无法使用该方法操作新视频。
针对运动放大参数难以确定的问题,Sushma等[8]提出了一种半自动的视频中微小运动放大 (Semi-Automated Magnification of subtle motion in video,SAM) 的方法,通过使用从前两帧视频获得的信息自动地确定这些参数。此外,在文献[9]中,通过功率谱估计得出运动的频率等信息,再进行自动运动放大处理,取得了不错的效果。但是,这些方法只能提取出一组固定的参数,这样提取的运动信息是不准确的。
借鉴现有的方法,本文提出一种基于时频分析的微小运动自动检测及放大方法,可以实现对非平稳微小运动的全自动检测及放大。首先, 从给定的视频序列中各像素点的时间波形入手,通过S变换 (Stockwell Transform)[10]求出各个时间波形在各个时间点的瞬时频率; 然后, 通过聚类方法[11]求出每一时间点的瞬时频率; 最后, 以此瞬时频率设计出一个滤波器来实现运动放大。
1 微小运动全自动放大框架本文对视频在时间上存在多频率分布的情况进行分析,以光点视频sims为例进行分析。视频sims是将文献[6]中提供的sim4视频的4种频率在时间域上合并而成,得出的光点视频总长为560帧,其中每140帧分别以2 Hz、3 Hz、5 Hz、7 Hz的频率左右摆动。图 1给出了本文多频率下运动信息自动检测及放大的整体框架,主要包含S变换和动态滤波。
具体来说,对于包含非平稳运动信息的视频,进行运动的自动检测及放大的算法流程如下:
步骤1 针对视频中的所有像素点,对于视频f(x, y, t) 中像素点 (i, j) 在时域信号fij(t) 进行S变换,求出关于信号fij(t) 的时频表示图;
步骤2 根据信号的时频表示计算关于点 (i, j) 在时间上瞬时频率值Fij(t);
步骤3 通过聚类方法对所有像素点的Fij(t) 进行运动层次分析,求出关于时间的瞬时中心频率值F(t);
步骤4 以瞬时中心频率值F(t) 和S变换的信号带宽B(t) 设计出对应的动态带通滤波器,以此实现非平稳运动信号的自动检测及放大。
2 信号时频分析及滤波器设计对于视频中运动信息的多频率分布的情况,属于非平稳信号的表现形式,直接对这类信号进行频谱分析会忽略信号频率在时间上的分布情况,因此需要分析信号在任意时刻的频率特征。对此,本文采用S变换这种时频表示方法对微小运动点的时间波形进行分析,并根据信号在时间轴上瞬时频率来设计动态的带通滤波器,从而实现非平稳运动信号的放大。
2.1 非平稳运动信号S变换对于信号频率分析的问题,傅里叶变换只给出了信号中频谱分量的信息,不能定位频率成分在信号中出现的位置。因此,应用时频表示 (Time Frequency Representation, TFR) 进行分析是比较合适的。现在有很多不同的时频表示方法,比如短时傅里叶变换 (Short-Time Fourier Transform, STFT)、连续小波变换 (Continuous Wavelet Transform, CWT) 等,而Stockwell[10]提出的S变换集中了短时傅里叶变换和小波变换的优点,免去了窗函数的选择和克服了窗宽固定的缺陷。因此,我们选择使用S变换这种时频表示方法。通过S变换,得出信号的频率及其对应的时间位置,然后提取出信号的频率,从而确定带通滤波器通带范围。
对于给定的时域信号x(t),其S变换的定义公式如下:
$S(t,f)=\int_{-\infty }^{\infty }{x(\tau )}\omega (t-\tau ){{\text{e}}^{-i2\pi f\tau }}\text{d}\tau $ | (1) |
其中:f表示频率,dτ表示时间位移量,τ为窗函数的中心点,S变换采用的窗函数为高斯窗函数,如式 (2) 所示:
$\omega (t)=\frac{\left| f \right|}{\sqrt{2\pi }}{{\text{e}}^{-{{t}^{2}}{{f}^{2}}/2}}$ | (2) |
微小运动频率的自动检测主要体现在时间序列上频率变化时,通过S变换确定信号在不同时刻的频率分布情况,从而检测出该运动频率及其对应的时间点。图 2以一个时间上不同频率的中心向四周灰度渐变的模糊化光点视频sims为例来分析非平稳运动的S变换 (其中fs为信号的采样频率。在本例中,视频帧率为30帧每秒即为信号采样频率fs,下同)。为了直观展示视频运动的效果,将视频切取一条线,然后沿着时间轴依次展开得到其时空切片图,图 2(b)给出了光点视频sims的时空切片图。
由图 2(d)可知,所有频率的信号在经过S变换后,图中都可以清晰地看出随时间变化而变化的情况,即通过S变换,我们可以知道某个频率及其所对应的时间位置。
通过对视频中N个像素点的时间波形进行S变换,可以得出N个T时间长度的频率矩阵,即N×T的矩阵; 然后在每个时间点上采用K-means聚类算法对所有频率值进行聚类分析,求出微小运动在每个时间点上的中心频率值,得出1×T的矩阵。在此基础上,通过确定视频在每一时刻运动信息的瞬时频率,来设计一个关于时间序列的动态带通滤波器,即为本文所需的时频滤波器。
2.2 动态滤波器的设计首先对S变换的逆变换进行分析,将S变换沿时间轴方向进行积分,将得到信号x(t) 的频谱X(f)。利用高斯窗函数所包含的面积等于1的特性,有:
$\int_{-\infty }^{\infty }{\omega (t\text{-}\tau ,f)}\text{d}\tau =\int_{-\infty }^{\infty }{\frac{\left| f \right|}{\sqrt{2\pi }}{{\text{e}}^{-{{(t-\tau )}^{2}}{{f}^{2}}/2}}}\text{d}\tau =1$ | (3) |
因此,将s(t, f) 沿着时间轴t积分,可得
$\begin{align} & \int_{-\infty }^{\infty }{s(t,f)\text{d}t}\text{=} \\ & \int_{-\infty }^{\infty }{x(\tau )}\left[ \int_{-\infty }^{\infty }{\frac{\left| f \right|}{\sqrt{2\pi }}{{\text{e}}^{-{{\left( t-\tau \right)}^{2}}{{f}^{2}}/2}}\text{d}\tau } \right]{{\text{e}}^{-i2\pi f\tau }}\text{d}t=X(f) \\ \end{align}$ | (4) |
式 (4) 表示S频谱是可逆的,同时也为S变换提供一个简单的逆变换:
$x\left( \tau \right)=\int_{-\infty }^{\infty }{\left[ \int_{-\infty }^{\infty }{S\left( t,f \right)}\text{d}t \right]{{\text{e}}^{i2\pi f\tau }}}\text{d}f=\int_{-\infty }^{\infty }{X\left( f \right){{\text{e}}^{i2\pi f\tau }}}\text{d}f$ | (5) |
由此可以设计一个基于S域的滤波器H(t, f) 来对信号x(t) 进行处理。
${{x}_{\text{filter}}}\left( \tau \right)=\int_{-\infty }^{\infty }{\left[ \int_{-\infty }^{\infty }{S\left( t,f \right)\cdot H\left( t,f \right)}\text{d}t \right]{{\text{e}}^{i2\pi f\tau }}}\text{d}f$ | (6) |
其中H(t, f) 表示时频滤波器,时频滤波器在此可以根据信号的瞬时频率设计成理想的时频二维滤波器或动态Butterworth滤波器。根据本文视频运动放大的需要,按照文献[12]设计出的所需的动态带通滤波器。
令fr(t) 表示信号用S变换分析所得的瞬时频率值,在视频序列中,时间变量t表示成离散变量,即为视频帧序号,则所需设计的滤波器表示如式 (7) :
$H\left( t,f \right)=\left\{ \begin{align} & 1,f\in \left[ fr\left( t \right)-B\left( t \right)/2,fr\left( t \right)+B\left( t \right)/2 \right] \\ & 0,其他 \\ \end{align} \right.$ | (7) |
其中B(t) 为时频滤波器的带宽,带宽B(t) 一般根据信号的能量在视频平面上沿瞬时频率的聚集程度来确定。本文中带宽B通过经验公式 (8) 求取:
$B=\left( {{F}_{\text{max}}}+F \right)/2-\left( {{F}_{\text{min}}}+F \right)/2$ | (8) |
式中:F为聚类得到的瞬时中心频率,Fmax、Fmin分别为瞬时中心频率所在类中的最大频率值和最小频率值,下同。对于图 2中给定点进行S变换后的波形,根据式 (7) 设计的动态理想带通滤波器结果如图 3(a)所示,滤波后的时频表示如图 3(b)所示,可以看出每个时刻点的有效频率得到保留,而其他干扰频率成分均被滤波掉了,然后对滤波后的信号进行放大处理,即可实现微小运动自动检测及放大的效果。
对将上述非平稳的微小运动检测方法、动态带通滤波器的设计方法与文献[6]的欧拉视频放大方法结合起来,得出基于非平稳微小运动自动检测及放大方法的总体框架如图 4所示。
该框架中,ωl和ωh代表带通滤波器的通带截频,通过S变换后的信号确定;α代表视频放大的倍数,λc表示α为0时的截断波长 (其确定方法沿用文献[5]中的方法)。
本文对视频进行非平稳微小运动自动放大处理包括以下几个步骤:
1) 参数估计:① 从视频中提取时间序列;② 对每个时间序列应用S变换这种时频表示进行分析;③ 估计带通滤波参数以及运动放大参数。
2) 运动放大:① 将视频序列进行拉普拉斯金字塔空域分解,得到不同空域分辨率的视频序列;② 对每个尺度的空域图像进行时域动态带通滤波处理,实时得到感兴趣的若干频带;③ 对感兴趣的频带信号进行线性放大;④ 通过将放大信号与原始信号进行重建处理,合成运动放大后的图像。
通过上述分析可知,本文首先从视频本身自动地检测出运动信息,然后进行动态带通滤波处理,最后进行运动放大,实现了非平稳微小运动的全自动检测及其放大。
4 实验结果及分析文中的实验素材主要来源于文献[6]提供的视频,部分实验素材来自于本文作者拍摄及处理的视频。首先验证本文方法同样适用于单一频率的运动目标; 其次对含有多种频率的视频进行实验,并分别对其抗噪性能进行分析。
4.1 估计的频率参数通过本文方法对所有实验视频进行实验,其中间过程的部分频率参数如表 1所示。
通过表 1的对比可以看出,与原始的EVM方法和SAM方法相比,本文方法得出的大部分视频的带通滤波器通频带范围与之基本保持一致或包含在其通频带频率范围内。根据带通滤波器频带的讨论,包含在通频带范围内信号得到放大,而相应的不在通频带范围的信号将会被衰减。
因此,通过本文方法使得检测出的运动信号的频率更为准确,从而设计出的带通滤波器包含中心频率且有相对较窄的通频带,在实现微小运动放大的同时能够更好地抑制噪声信号的放大。
表 1中数据只是某一时刻的,而所有的滤波参数是根据运动目标的中心频率动态调整的,其滤波效果更为准确。准确的带通滤波参数在实现微小运动放大时能够抑制噪声信号,得出更好的运动放大效果。
4.2 单一频率视频实验结果首先对单一运动频率的视频flower进行实验分析,flower视频是一个在室内相对较暗的光照下进行轻微运动的太阳能摇摆花,该视频是我们自己拍摄的未知运动放大参数的视频。由于不知道运动放大参数,很难得出带通滤波器的通带范围,用原始的欧拉视频放大方法难以将该视频进行运动放大处理。因此,本文应用提出的动态滤波方法进行运动信息检测及放大,实验结果如图 5所示。
实验结果表明,本文基于非平稳信号动态滤波的方法能够实现对未知运动信息的视频进行自动检测及放大,并且取得了不错的运动放大效果,由运动放大前后的时空切片图可以清晰地看出运动细节的变化。接下来对时间域上存在多频率分布的视频进行实验分析。
4.3 多频率视频实验结果视频sims是一个包含了4种频率运动信息的视频,不同时刻对应的频率可能不相同,用欧拉视频放大方法可以使用包括所有运动运动频率的带通滤波器来实现,但由于通频带较宽,会将更多含有噪声信息的信号放大。而本文基于S变换的动态滤波方法能够实时地根据检测出的频率调整带通滤波器的通频带范围,使得检测出的运动信息更准确,避免了过宽的通频带下更多的噪声被放大。图 6给出了使用欧拉视频放大方法和本文方法对视频sims进行放大的实验结果,其中欧拉方法的放大倍数分别取20和50,EVM方法中滤波器带宽选择[1 Hz, 8 Hz],因为该通频带包含了视频中的所有频率。
由图 6可知,两种方法均可实现微小运动放大。比较图 6(b)、(d)可以看出本文的自动频率检测方法可以同欧拉放大一样对视频进行放大,并且无需进行频率参数选择,可以自动地进行动态滤波。比较图 6(c)、(e),当放大倍数增加到50倍时,由于欧拉方法所取的通频带较宽,使得放大后的结果在正常信号部分中出现黑影现象,与原文献[6]中放大倍数限制的问题相对应,而本文动态滤波方法在每一时间点均自动设计合适的带通滤波器,可以有效保证检测出的运动信息的准确性,因而呈现出更好的视觉放大效果。此外,视频中不同频率连接处出现黑影现象,这是由于视频在连接处频率变换不平缓造成的。
接下来对三个单一频率的视频 (baby、shadow和wrist) 拼接在一起的时间上多频率视频BSW进行实验分析,BSW视频的1到302帧是baby,303到482帧是shadow,483到838帧是wrist。这样该视频在时间上将呈现多频率分布,利用本文方法对其进行运动放大,其结果如图 7所示。
由图 7可知,微小运动细节得到了明显的增强。对比图 7(b)、(d)的XT切片图可知,在视频BSW中的三段不同频率的视频均被进行了运动放大,并且可以很直观地看出运动放大后的效果。因此,本文方法对合成的视频BSW取得了不错的运动放大效果,同时也证明了该方法对放大在时间上存在多频率的视频的有效性。该视频应用欧拉视频放大处理,同样只能让通频带包含这三种频率 (同图 6实验),再进行放大处理。这样的处理方法通频带较宽,会将更多含有噪声的信号放大,而该处理方法都是建立在已知运动频率的情况下,对于未知视频,欧拉视频放大方法将无法获取滤波器参数。
4.4 抗噪性能分析为了定量地比较视频放大的抗噪性能,本文基于现有的图像信噪比 (Signal-to-Noise Ratio, SNR) 评价算法[14-15],自定义一种视频区域信噪比SNR的评价方法, 该SNR定义为选择视频中的一个平滑区域,其方差作为噪声功率,平均值作为有用信号,从而计算出该选定区域的信噪比SNR。接下来详细介绍具体计算过程。
视频的每一帧图像M×N可以看作图像信息和噪声的叠加,如式 (9) 所示:
$D\left( x,y \right)=S\left( x,y \right)+N\left( x,y \right)$ | (9) |
在第一帧视频的光滑区域选择M1×N1大小的块,记作Dpatch(i, j),计算该区域的平均值:
$Patc{{h}_{\rm{mean}}}=\frac{1}{{{M}_{1}}{{N}_{1}}}\sum\limits_{i=1}^{{{M}_{1}}}{\sum\limits_{j=1}^{{{N}_{1}}}{{{\mathit{\boldsymbol{D}}}_{\rm{patch}}}\left( i,j \right)}}$ | (10) |
其中该区域噪声均方误差为:
$Patc{{h}_{\rm{MSE}}}=\frac{1}{{{M}_{1}}{{N}_{1}}}\sum\limits_{i=1}^{{{M}_{1}}}{\sum\limits_{j=1}^{{{N}_{1}}}{\{{{\mathit{\boldsymbol{D}}}_{\rm{patch}}}(i,j)}}-Patc{{h}_{\rm{mean}}}{{\}}^{2}}$ | (11) |
因此,第一帧视频选定区域的信噪比SNR表示为:
$SNR = 10{\kern 1pt} \;{\rm{lg}}\frac{{\frac{1}{{MN}}\sum\limits_{i = 1}^{{M_1}} {\sum\limits_{j = 1}^{{N_1}} {\mathit{\boldsymbol{D}}_{{\rm{patch}}}^2(i,j)} } }}{{Patc{h_{{\rm{MSE}}}}}}$ | (12) |
为了评价放大处理对视频的噪声放大结果,根据式 (12) 计算视频sims放大后相对原视频的信噪比SNR的差值。得出的SNR差值越大,表明该方法的抗噪性能越好,视频放大对噪声敏感性影响越弱。
为了证明该评价方法的有效性,选用baby视频进行论证,讨论EVM方法在不同放大倍数下的结果,由于baby视频帧数过多,选取其中的第1帧~第250帧的结果如图 8所示。由图 8中视频baby的实验表明,视频放大20倍时的信噪比SNR比放大10倍时的值要小得多,即放大20倍时的视频噪声相对较大。对于EVM方法,随着放大倍数的增加,视频的噪声信号也被线性地放大,视频baby的SNR实验与此相符,证明了该方法的有效性。
对于sims视频的实验,信噪比SNR的结果如图 9所示,同样选取其中的第1帧~第250帧的结果作为示例。由图 9可知,在相同的放大倍数 (α=20) 情况下,本文基于S变换的动态滤波方法比原始的欧拉方法具有更好的抗噪性能。因为基于S变换的动态带通滤波方法具有更精确的滤波范围,能够精确地放大每个频率的信号,而使用欧拉方法放大该视频,只能选择使用包含该视频中的所有频率的通频带1 Hz~8 Hz,通频带较宽,放大了额外的噪声信息。
通过视频flower、视频sims以及视频BSW的实验结果可知,本文基于S变换的动态滤波放大方法不仅适用于单一频率或者频率变化不大的视频,而且对多频率的视频也能动态地进行放大处理,并且能够实现对未知的视频进行运动放大,因此具有更好的实用性。
4.5 时间复杂度分析对本文方法的计算复杂度进行定量的分析,并将其与己有方法的计算复杂度进行对比分析,如表 2所示。
由表 2可知,相比原始的EVM方法与SAM方法,本文方法的计算复杂度有所增加,其原因是计算运动信号频率以及进行动态滤波时增加了计算复杂度。然而,本文增加一点计算复杂度换来的优势是能够对包含非平稳运动的视频进行自动检测及放大,不需要人为选取运动放大参数,这与需要用户耗费时间去手动找到合适的参数相比,更具有实际使用价值。此外,对于微小运动放大计算复杂度的优化,文献[16]中已经提出一种加速算法,可以更快地实现微小运动的放大。因此,对于寻找参数耗时的问题,可以通过其他优化方法解决。
5 结语本文针对视频在时间域上包含多频率分布的情况,提出了基于S变换动态滤波的自动检测及放大视频中非平稳微小运动的方法。通过S变换,对不同时刻呈现不同频率的运动信号进行分析处理,得出随时间变化的动态频率值,以此设计出动态带通滤波器,最后实现运动放大处理,动态地展现视频中微小运动的运动情况,并在一定程度上抑制噪声干扰。此外,基于现有的一些图像信噪比评价算法,自定义一种视频的SNR来分析本文方法的抗噪性能。实验结果表明,本文方法在实际视频放大处理中得出很好的运动放大效果,不仅适用于单一频率或者频率变化不大的视频,而且对多频率的视频放大效果也比较理想,具有更好的实用性。后续将继续研究非平稳信号运动放大的有关问题,并从空域多分辨率角度进行讨论。
[1] | RUBINSTEIN M, WADHWA N, DURAND F, et al. Revealing invisible changes in the world[J]. Science, 2013, 339(6119): 519-519. |
[2] | PARK S, KIM D. Subtle facial expression recognition using motion magnification[J]. Pattern Recognition Letters, 2009, 30(7): 708-716. doi: 10.1016/j.patrec.2009.02.005 |
[3] | BALAKRISHNAN G, DURAND F, GUTTAG J. Detecting pulse from head motions in video[J]. Computer Vision and Pattern Recognition, 2013, 9(4): 3430-3437. |
[4] | POH M Z, MCDUFF D J, PICARD R W. Non-contact, automated cardiac pulse measurements using video imaging and blind source separation[J]. Optics Express, 2010, 18(10): 10762-10774. doi: 10.1364/OE.18.010762 |
[5] | LIU C, TORRALBA A, FREEMAN W T, et al. Motion magnification[J]. ACM Transactions on Graphics, 2005, 24(7): 519-526. |
[6] | WU H, RUBINSTEIN M, SHIH E, et al. Eulerian video magnification for revealing subtle changes in the world[J]. ACM Transactions on Graphics, 2012, 31(4): Article No. 65. |
[7] | WADHWA N, RUBINSTEIN M, DURAND F, et al. Phase-based video motion processing[J]. ACM Transactions on Graphics, 2013, 32(4): Article No. 80. |
[8] | SUSHMA M, GUPTA A, SIVASWAMY J. Semi-automated magnification of small motions in videos[C]//PReMI 2013:Pattern Recognition and Machine Intelligence, LNCS 8251. Berlin:Springer, 2013:417-422. |
[9] | 雷林, 李乐鹏, 李准, 等. 自动检测及放大视频中的微小运动[J]. 小型微型计算机系统, 2016, 37(9): 2120-2124. ( LEI L, LI L P, LI Z, et al. Automated detection and magnification of small motion in videos[J]. Journal of Chinese Computer Systems, 2016, 37(9): 2120-2124. ) |
[10] | STOCKWELL R G, MANSINHA L, LOWE R P. Localization of the complex spectrum:the S transform[J]. IEEE Transactions on Signal Processing, 1996, 44(4): 998-1001. doi: 10.1109/78.492555 |
[11] | 张建辉. K-means聚类算法研究及应用[D]. 武汉: 武汉理工大学, 2007. ( ZHANG J H. Research an application of K-means clustering algorithm[D]. Wuhan:Wuhan University of Technology, 2007. ) |
[12] | 赵淑红, 朱光明. S变换时频滤波去噪方法[J]. 石油地球物理勘探, 2007, 42(4): 402-406. ( ZHAO S H, ZHU G M. Time-frequency filtering to denoise by S transform[J]. Oil Geophysical Prospecting, 2007, 42(4): 402-406. ) |
[13] | RUBINSTEIN M. Analysis and visualization of temporal variations in video[D]. Cambridge:Massachusetts Institute of Technology, 2014. |
[14] | TURAGA D S, CHEN Y, CAVIEDES J. No reference PSNR estimation for compressed pictures[C]//Proceedings of the 2002 International Conference on Image Processing. Piscataway, NJ:IEEE, 2004:173-184. |
[15] | 杨柱中, 周激流, 郎方年. 用噪声检测算法改进理想低通滤波器[J]. 计算机应用, 2014, 34(10): 2971-2975. ( YANG Z Z, ZHOU J L, LANG F N. Improving ideal low-pass filter with noise detection algorithm[J]. Journal of Computer Applications, 2014, 34(10): 2971-2975. doi: 10.11772/j.issn.1001-9081.2014.10.2971 ) |
[16] | 李乐鹏, 雷林, 孙水发, 等. 视频微小运动放大的加速方法[J]. 计算机工程与应用, 2015, 51(24): 195-200. ( LI L P, LEI L, SUN S F, et al. Improved video small motion magnification processing[J]. Computer Engineering and Applications, 2015, 51(24): 195-200. doi: 10.3778/j.issn.1002-8331.1503-0166 ) |