CN102707340A - 一种基于视频图像的降雨量测量方法 - Google Patents

Abstract

一种基于视频图像的降雨量测量方法，对视频摄像机采集的降雨视频，采用对数图像处理LIP方法，用灰色调替代亮度值，建立灰色调约束条件提取候选雨滴；利用主成分分析法PCA构建候选雨滴的倾斜方向的概率密度分布函数，并推导求取该分布函数峰值的Mean-shift方法，通过约束候选雨滴的倾斜方向去除干扰；通过场景背景和雨滴前景的色彩向量夹角约束去除散焦雨滴；利用摄像机参数标定计算每个聚焦雨滴的尺寸，构建雨滴谱，利用Gamma雨滴谱分布模型拟合推导降雨量。本发明打破传统的降雨量测量方法，创新的提出基于视频的测量方法，获得的降雨量数据具有很高的时间分辨率、空间分辨率，经济快捷。

Description

一种基于视频图像的降雨量测量方法

技术领域

本发明属于机器视觉、视频图像处理领域，特别涉及基于视频图像的降雨量测量方法，主要用于自动化气象测量中，为一种基于视频图像的降雨量测量方法。

背景技术

降雨量数据一般通过雨量计或气象雷达获得。光学雨量计^[1]感知经过雨的激光光束的衰减程度来测定降雨量；微波雨量计^[2]测量由雨滴导致的多普勒偏移，利用雨滴终极速度和雨滴尺寸大小的关系估计降雨量；还有一种常用的翻斗式雨量计^[3]，当翻斗内积水量达到一定量时，翻斗倾倒并输送一个脉冲信号，从而记录降雨量。气象雷达^[4]测量降雨量，首先建立雷达反射率和降雨量的关系，然后通过反演计算推导降雨量。

不同于以上方法，本发明提出利用视频图像测量降雨量。目前，监控摄像机日益普及，分布密度较大，性能也逐步增强，采用基于视频的方法获得的降雨量数据具有较大的空间分辨率；另一方面，视频图像通常数十帧/秒，因而获得的降雨量数据也具有很高的时间分辨率。采用视频图像进行降雨量测量主要涉及两方面内容：一是从视频中检测出雨滴，二是从图像中雨滴占据的像素个数推导出真实的降雨量数值。

视频雨滴检测有多种方法。文献^[5]研究了雨的动力学和光学特性并提出利用时空相关性和光度模型检测雨滴，通过拟合同一雨滴中的像素在受雨滴影响前后的亮度差值与背景亮度的线性比值辨别是否雨滴，但雨滴通常面积较小，拟合易受噪声干扰；文献^[6][7]对整个视频利用聚类方法辨别雨滴和背景，不能实时处理。已有文献中，进行雨滴检测并用于测量视频降水量的研究甚少。为进行降雨量计算，更为关键的是在检测出的雨滴中区分聚焦雨滴和散焦雨滴，而现有的视频雨滴检测方法很少考虑这方面内容。文献^[8]利用雨滴长宽比约束来区分聚焦和散焦雨滴，但是雨滴宽度通常较小，测量误差较大。

参考文献

[1]F.V.Brock and S.J.Richardson,Meteorological Measurement Systems.New York:Oxford Univ.Press,2001.

[2]T.J.Mansheim,A.Kruger,J.Niemeier,etc.“A Robust Microwave Rain Gauge,”IEEETrans.Instrument and Measurement,vol.59,no.8,pp:2204-2210,2010.

[3]E.Habib,W.F.Krajewski,and A.Kruger,“Sampling errors of tipping bucket rain gaugemeasurements,”ASCE J.Hydrol.Eng.,vol.6,no.2,pp.159–166,Mar./Apr.2001.

[4]R.Teschl,W.L.Randeu,and F.Teschl.”Weather Radar Estimates of Rainfall Adjustedto Rain Gauge Measurements Using Neural Networks,”In Proc.of Int.Joint Conf.onNeural Networks,Vancouver,Canada,2006,pp:5126-5131.

[5]Garg K,Nayar S K.Detection and removal of rain from videos.In Proc.2004 IEEEComputer Society Conference on Computer Vision and Pattern Recognition,Washington,2004:I-528-I-535.

[6]Zhang X P,Li H,Qi Y Y,et al.Rain removal in video by combining temporal andchromatic properties.In Proceedings og International Conference on Multimedia andExpo,Toronto,Canada,2006:461-464.

[7]刘鹏,徐晶,刘家锋,唐降龙.一种受雨滴污染视频的快速分析方法.自动化学报.2010,36(10):1371-1378.

[8]Garg K,Nayar S K.When does a camera see rain?In Proc of the 10th Int.conf.oncomputer vision.Beijing,China,2005:1067-1074.

[9]郑娇恒，陈宝君.雨滴谱分布函数的选择：M-P分布和Gamma分布的对比研究.气象科学.2007,27(1):17-19.

[10]Safaee Rad R,Smith K C,Benhabib B,Tchoukanov I.Application of Moment andFourier Descriptors to the Accurate Estimation of Elliptical Shape Parameters.PatternRecognition Letters,1992,13(7):2465-2468.

发明内容

本发明要解决的问题是：现有的降雨量测量方法时间和空间分辨率较低，现有的雨滴检测方法易受噪声干扰、实时性不足，现有的聚焦散焦雨滴区分方法误差较大。

本发明的技术方案为：一种基于视频图像的降雨量测量方法，对视频摄像机采集的降雨视频，采用对数图像处理LIP方法，用灰色调替代亮度值，建立灰色调约束条件提取候选雨滴；利用主成分分析法PCA构建候选雨滴的倾斜方向的概率密度分布函数，并推导求取该分布函数峰值的Mean-shift方法，通过约束候选雨滴的倾斜方向去除干扰；通过场景背景和雨滴前景的色彩向量夹角约束去除散焦雨滴；利用摄像机参数标定计算每个聚焦雨滴的尺寸，构建雨滴密度分布函数，利用Gamma雨滴谱分布模型拟合推导降雨量，具体如下：

1）检测候选雨滴：利用视频中连续三帧图像：第n－1帧，第n帧和第n＋1帧，在LIP框架下建立灰色调约束，提取得到第n帧图像的候选雨滴；

2）构建雨滴倾斜方向概率分布密度函数：利用PCA方法，提取第n帧图像每个候选雨滴的倾斜方向θ_i以及长a_i、宽b_i，i＝{1,2,…,P_n}，P_n为第n帧图像中检测出的候选雨滴总个数，然后构建雨滴倾斜方向概率密度分布函数f(θ)；

3）确定降雨方向，去除噪声干扰：针对步骤2)中获得的f(θ)，用Mean-shift算法提取它的峰值，所述峰值即为降雨方向，对候选雨滴实施方向约束，不符合降雨方向的判定为噪声；

4）去除散焦雨滴：计算候选雨滴中每个像素的场景背景和雨滴前景的色彩向量夹角γ(x,y)，对每个候选雨滴，检查它所包含像素的平均γ值是否小于设定的阈值，若是则判定为散焦雨滴，否则为聚焦雨滴；

5）构建雨滴密度分布函数：利用视频摄像机参数计算每个聚焦雨滴的尺寸，构建雨滴密度分布函数N(d)；

6）利用步骤5)中的雨滴密度分布函数N(d)，以及Gamma雨滴谱分布模型拟合推导降雨量R。

步骤1)中的灰色调约束是：

f_n(x,y)<f_n-1(x,y)&f_n(x,y)<f_n+1(x,y)    (1)

|f_n(x，y)Θf_n-1(x，y)|_E＞th & |f_n(x，y)Θf_n+1(x，y)|_E＞th    (2)

其中f_n-1(x，y)、f_n(x，y)、f_n+1(x，y)是n-1、n、n+1连续三帧图像的灰色调函数，th为灰色调阈值，同时满足式(1)、(2)的像素则判定为候选雨滴；

灰色调函数f(x，y)与入射光强度即输入图像F(x，y)的关系为：

$f(x,y)=M(1-\frac{F(x,y)}{F_{\max }})---\left(3\right)$

其中F_max是人眼视觉系统的饱和亮度值，取为255，M是自定义的灰色调函数值域的最大值，取255，由式（3）和连续三帧图像，得到三帧图像各自的灰色调函数；

式(2)中Θ代表LIP框架下灰色调空间的减法操作，两个灰色调函数f(x，y)和g(x，y)的减法定义为：

$f(x,y)\mathrm{\&Theta;g}(x,y)=M\frac{f(x,y)-g(x,y)}{M-g(x,y)}---\left(4\right)$

式(2)中“| |_E”表示LIP框架下灰色调空间的取模操作，定义为：

${\left|f(x,y)\right|}_E=\left\{\begin{array}{ccc}f(x,y);& \mathrm{if}& f(x,y)\&GreaterEqual;0\\ \mathrm{\&Theta;f}(x,y);& \mathrm{if}& f(x,y)<0\end{array}\right.---\left(5\right).$

步骤2)中提取候选雨滴的倾斜方向θ_i以及长a_i、宽b_i的方法是：

令N_i表示第i个雨滴占有的像素个数，X_i为N_i*2的矩阵：

$X_i=\left[\begin{array}{cc}{x}_{i1}& y_{i1}\\ x_{i2}& y_{i2}\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ x_{{\mathrm{iN}}_i}& y_{{\mathrm{iN}}_i}\end{array}\right]---\left(6\right)$

其中（x_ik，y_ik），是该雨滴中第k个像素相对于该雨滴质心的坐标，k=1,2,…N_i，设向量u₁ ⁱ=[u_1x ⁱ,u_1y ⁱ]^T,u₂ ⁱ＝[u_2x ⁱ,u_2y ⁱ]^T是矩阵X_i ^TX_i的两个特征向量，分别对应于特征值λ₁ ⁱ,λ₂ ⁱ，λ₁ ⁱ＞λ₂ ⁱ，u₁ ⁱ即表示第一主成分方向，则雨滴的倾斜方向θ_i为：

${\mathrm{\&theta;}}_i={\tan }^{-1}\left(\frac{{u_{1y}}^i}{{u_{1x}}^i}\right)---\left(7\right)$

长a_i、宽b_i为：

$a_i=4\sqrt{{{\mathrm{\&lambda;}}_1}^i},$ $b_i=4\sqrt{{{\mathrm{\&lambda;}}_2}^i}---\left(8\right)$

雨滴倾斜方向概率密度分布函数f(θ)的构建为：

$f\left(\mathrm{\&theta;}\right)=\underset{i=1}{\overset{P_n}{\mathrm{\&Sigma;}}}\frac{N_i}{\sqrt{2\mathrm{\&pi;}}{h}_i}k\left({\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2\right)---\left(9\right)$

其中k(X)为核函数， $k\left(X\right)=e^{-\frac{1}{2}X},$ 对应式（9）则 $X={\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2,$ 核函数带宽h_i用于描述θ_i的不确定性，计算方法为：

$h_i=C\frac{b_i}{a_i}+H---\left(10\right)$

C是将视频图像的长宽比值转换至核函数带宽时的尺度缩放，H表示核函数基本带宽，C、H均为常数。

步骤3)中求取f(θ)峰值的Mean-shift方法为：

步骤1:设定迭代初始值θ₀=0;

步骤2:令θ₁＝m(θ₀)+θ₀;

步骤3:若|θ₁-θ₀|<ε，令降雨方向θ_main=θ₁;否则令θ₀=θ₁并回到步骤2进行循环判断；

其中，

$m\left(\mathrm{\&theta;}\right)=\frac{\underset{i=1}{\overset{P_n}{\mathrm{\&Sigma;}}}\frac{2N_i}{\sqrt{2\mathrm{\&pi;}}{h_i}^3}{\mathrm{\&theta;}}_{i}g\left({\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2\right)}{\underset{i=1}{\overset{P_n}{\mathrm{\&Sigma;}}}\frac{2N_i}{\sqrt{2\mathrm{\&pi;}}{h_i}^3}g\left({\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2\right)}-\mathrm{\&theta;}---\left(12\right)$

g(X)=-k′(X)

-k'(X)是核函数k(X)的导数再取负，核函数 $k\left(X\right)=e^{-\frac{1}{2}X},$ 对应式（12）有 $X=$ ${\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2;$

根据降雨方向θ_main，对候选雨滴实施方向约束：

如果第i个雨滴的方向θ_i相对于降雨方向θ_main的差值大于设定的阈值θ_T，则判定为干扰噪声：

|θ_i-θ_main|>θ_T    (13)

阈值θ_T设为10°。

步骤4)中色彩向量夹角γ(x,y)计算方法为：

$\mathrm{\&gamma;}(x,y)={\cos }^{-1}\left(\frac{{\stackrel{\&RightArrow;}{I}}_r(x,y)\&CenterDot;{\stackrel{\&RightArrow;}{I}}_b(x,y)}{\left|{\stackrel{\&RightArrow;}{I}}_r(x,y)\right|\left|{\stackrel{\&RightArrow;}{I}}_b(x,y)\right|}\right)---\left(16\right)$

其中，为当前帧被雨滴遮挡的像素的色彩向量，是像素未被雨滴遮挡时的色彩向量，

直接从当前帧图像获得，

通过对前一帧和后一帧的像素色彩向量

取平均获得：

${\stackrel{\&RightArrow;}{I}}_b(x,y)=\frac{{\stackrel{\&RightArrow;}{F}}_{n-1}(x,y)+{\stackrel{\&RightArrow;}{F}}_{n+1}(x,y)}{2}---\left(17\right).$

摄像机参数满足z₀≤fd_min，f为摄像机焦距，z₀为焦平面距离，d_min为期望捕获的最小雨滴的直径，取为0.3mm；

根据雨滴的终极速度v和其直径d的关系：

步骤5)中雨滴密度分布函数N(d)构建为：

对于落在焦平面上的雨滴，即聚焦雨滴，其直径d_i为：

$d_i=2{\left(\frac{z_0{a}_i}{200\mathrm{fT}}\right)}^2---\left(19\right)$

T为摄像机曝光时间，

得到每个检测到的聚焦雨滴的直径d_i后，计算直径d_i落在[d-1/2Δd，d+1/2Δd]区间内的雨滴数目，从而构建雨滴密度分布函数N(d)，Δd是单位尺度间隔，取1mm，利用连续多帧图像进行统计，以获得稳定鲁棒的N(d)。

步骤6)中降雨量计算方法为根据雨滴密度分布函数N(d)，利用Gamma雨滴谱分布模型拟合推导降雨量R：

$R=0.0071N_0\frac{\mathrm{\&Gamma;}(\mathrm{\&mu;}+4.67)}{{\mathrm{\&omega;}}^{\mathrm{\&mu;}+4.67}}---\left(26\right)$

其中，N₀,μ和ω是Gamma雨滴谱分布模型的三个参数，根据雨滴密度分布函数N(d)求得：

$\mathrm{\&mu;}=\frac{3M_4{M}_2-4{M_3}^2}{({M_3}^2-M_4{M}_2)}---\left(23\right)$

$\mathrm{\&omega;}=\frac{(\mathrm{\&mu;}+4)M_3}{M_4}---\left(24\right)$

$N_0=\frac{{\mathrm{\&omega;}}^{\mathrm{\&mu;}+4}{M}_3}{\mathrm{\&Gamma;}(\mathrm{\&mu;}+4)}---\left(25\right)$

$M_i=\underset{0}{\overset{8}{\&Integral;}}{d}^{i}N\left(d\right)\mathrm{dd},i=\mathrm{2,3,4}---\left(22\right).$

本发明补充现有的基于雨量计或气象雷达的降雨量测量方法，提出一种基于视频的测量方法，具有较高的时间和空间分辨率；针对视频雨滴检测算法存在的易受噪声干扰、实时性不足等问题，提出采用符合人眼视觉的对数图像处理LIP(Logarithmic ImageProcessing)方法，用灰色调替代传统的亮度值，分析受雨滴影响的像素光度特征，建立灰色调约束条件提取候选雨滴，同时，利用主成分分析法PCA构建候选雨滴的倾斜方向的概率密度分布函数，并推导求取该概率密度分布函数峰值的Mean-shift方法，通过约束候选雨滴的倾斜方向去除干扰，本发明基于灰色调和倾斜方向双重约束的雨滴检测方法，可实时进行并有效去除噪声干扰；另外，针对聚焦和散焦雨滴辨别问题，提出色彩向量夹角约束，可有效区分两者，从而为视频降雨量计算提供准确的雨滴数据。

本发明创新的提出了一种基于视频图像的降雨量测量方法，利用雨滴检测算法，检测出视频图像中的候选雨滴，并去除噪声干扰和散焦雨滴，利用摄像机参数标定计算每个聚焦雨滴的尺寸，构建雨滴密度分布函数，即雨滴谱(Raindrop Size Distribution,RSD)，再利用Gamma雨滴谱分布模型拟合推导降雨量。本发明打破传统的降雨量测量方法，创新的提出基于视频的测量方法，获得的降雨量数据具有很高的时间分辨率、空间分辨率，经济快捷。

附图说明

图1为本发明降雨量测量算法流程图。

图2为本发明降雨场景中某像素的时序灰色调变化。

图3为本发明图像中聚焦雨滴和散焦雨滴的色彩向量差别。

图4为本发明摄像机参数设置需满足的条件。

图5为小雨情况下获得的雨滴谱和拟合的雨滴谱。

图6为中雨情况下获得的雨滴谱和拟合的雨滴谱。

图7为大雨情况下获得的雨滴谱和拟合的雨滴谱。

具体实施方式

本发明通过视频处理的方法构建雨滴谱RSD，再利用气象学上常用的Gamma分布模型拟合推导降雨量。另外，为了精确估计降雨量，需要区分聚焦与散焦雨滴，本发明从光学原理出发，推导了区分聚焦与散焦雨滴的准则。

本发明降雨量测量算法流程如图1所示，其中F_n-1、F_n、F_n+1为连续三帧亮度图像，详细实施方式如下：

1.利用连续三帧图像检测出当前帧的雨滴。被雨滴遮挡后的像素亮度往往高于该处原先的背景亮度，并且雨滴降落速度较快，因而同一像素位置在相邻两帧很少同时被雨滴遮挡，其时序灰色调变化呈抖动状态，如图2所示，因此同时满足式(1)、(2)的像素可被判定为候选雨滴：

f_n(x,y)<f_n-1(x,y)&f_n(x,y)<f_n+1(x,y)    (1)

|f_n(x，y)Θf_n-1(x，y)|_E＞th & |f_n(x，y)Θf_n+1(x，y)|_E＞th    (2)

其中f_n-1(x，y)、f_n(x，y)、f_n+1(x，y)是n-1、n、n+1连续三帧图像的灰色调函数，th为设定的灰色调阈值.LIP框架中灰色调函数与入射光强度即输入图像F(x，y)的关系如式(3)所示：

$f(x,y)=M(1-\frac{F(x,y)}{F_{\max }})---\left(3\right)$

其中F_max是人眼视觉系统的饱和亮度值，一般情况取为255。M是自定义的灰色调函数值域的最大值，本发明中取255，由式（3）和连续三帧图像，得到三帧图像各自的灰色调函数。

式(2)中Θ代表LIP框架下灰色调空间的减法操作，灰色调函数f(x，y)和g(x，y)的减法定义为：

$f(x,y)\mathrm{\&Theta;g}(x,y)=M\frac{f(x,y)-g(x,y)}{M-g(x,y)}---\left(4\right)$

式(2)中| |_E表示LIP框架下灰色调空间的取模操作：

${\left|f(x,y)\right|}_E=\left\{\begin{array}{ccc}f(x,y);& \mathrm{if}& f(x,y)\&GreaterEqual;0\\ \mathrm{\&Theta;f}(x,y);& \mathrm{if}& f(x,y)<0\end{array}\right.---\left(5\right)$

|fΘg|_E实际代表了灰色调f和g的欧式距离。采用灰色调函数运算而不直接采用亮度图像运算的优势在于，亮度图像运算譬如亮度减法在整个亮度范围内均匀线性操作，而灰色调函数运算如灰色调距离随图像灰色调值自适应调节，具有非线性特征，与人眼视觉规律相符合。

2.利用PCA方法构建雨滴方向分布概率密度函数。令N_i表示第i个雨滴占有的像素个数，X_i为N_i*2的矩阵：

$X_i=\left[\begin{array}{cc}{x}_{i1}& y_{i1}\\ x_{i2}& y_{i2}\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ x_{{\mathrm{iN}}_i}& y_{{\mathrm{iN}}_i}\end{array}\right]---\left(6\right)$

其中（x_ik，y_ik），k＝1,2,…Ni是该雨滴中第k个像素相对于该雨滴质心的坐标.设向量u₁ ⁱ=[u_1x ⁱ,u_1y ⁱ]^T,u₂ ⁱ＝[u_2x ⁱ,u_2y ⁱ]^T是矩阵X_i ^TX_i的两个特征向量，分别对应于特征值λ₁ ⁱ,λ₂ ⁱ(λ₁ ⁱ＞λ₂ ⁱ)，u₁ ⁱ即表示第一主成分方向，则雨滴的倾斜方向θ_i可由式(7)计算.

${\mathrm{\&theta;}}_i={\tan }^{-1}\left(\frac{{u_{1y}}^i}{{u_{1x}}^i}\right)---\left(7\right)$

另外，雨滴的长a_i和宽b_i可通过椭圆拟合得到（参见参考文献[10]）：

$a_i=4\sqrt{{{\mathrm{\&lambda;}}_1}^i},$ $b_i=4\sqrt{{{\mathrm{\&lambda;}}_2}^i}---\left(8\right)$

构建降雨方向的概率密度分布函数f(θ)，由高斯核函数叠加的形式给出：

$f\left(\mathrm{\&theta;}\right)=\underset{i=1}{\overset{P_n}{\mathrm{\&Sigma;}}}\frac{N_i}{\sqrt{2\mathrm{\&pi;}}{h}_i}k\left({\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2\right)---\left(9\right)$

其中P_n为第n帧图像中检测出的候选雨滴个数，k(X)为核函数，

对应式(9)则

核函数带宽h_i用于描述θ_i的不确定性，考虑到越接近圆形的雨滴，其方向估计越难准确，因此采用式(10)来定义该不确定性：

$h_i=C\frac{b_i}{a_i}+H---\left(10\right)$

其中C、H均为常数，C是将长宽比值转换至核函数带宽时的尺度缩放，H表示核函数的基本带宽。

3.用Mean-shift算法提取降雨方向，对候选雨滴区域实施方向约束，去除噪声干扰。式(9)可视为一维核密度函数，Mean-shift方法是求取其极值的有效方法.但是不同于一般的核密度函数，式(9)中核函数带宽h_i随θ_i而变，推导Mean-shift向量m(θ)：

$\&dtri;f\left(\mathrm{\&theta;}\right)=\underset{i=1}{\overset{P_n}{\mathrm{\&Sigma;}}}\frac{2N_i}{\sqrt{2\mathrm{\&pi;}}{h_i}^3}({\mathrm{\&theta;}}_i-\mathrm{\&theta;})g\left({\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2\right)$

                                   (11)

$=\left[\underset{i=1}{\overset{P_n}{\mathrm{\&Sigma;}}}\frac{2N_i}{\sqrt{2\mathrm{\&pi;}}{h_i}^3}g\left({\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2\right)\right]\&times;m\left(\mathrm{\&theta;}\right)$

$m\left(\mathrm{\&theta;}\right)=\frac{\underset{i=1}{\overset{P_n}{\mathrm{\&Sigma;}}}\frac{2N_i}{\sqrt{2\mathrm{\&pi;}}{h_i}^3}{\mathrm{\&theta;}}_{i}g\left({\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2\right)}{\underset{i=1}{\overset{P_n}{\mathrm{\&Sigma;}}}\frac{2N_i}{\sqrt{2\mathrm{\&pi;}}{h_i}^3}g\left({\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2\right)}-\mathrm{\&theta;}---\left(12\right)$

g(X)＝-k′(X)

-k′(X)是核函数k(X)的导数再取负，核函数 $k\left(X\right)=e^{-\frac{1}{2}X},$ 对应式(12)有 $X=$ ${\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2;$

降雨降水方向θ_main通过如下方式获得：

Step1：设定迭代初始值θ₀＝0；

Step2：令θ₁＝m(θ₀)+θ₀；

Step3：若|θ₁-θ₀|＜ε，停止循环并令θ_main＝θ₁；否则令θ₀＝θ₁并回到step2.

如果第i个雨滴的方向θ_i相对于降雨方向θ_main的差值大于设定的阈值θ_T，则判定为干扰噪声：

|θ_i-θ_main|>θ_T    (13)

阈值θ_T设为10°。

4.通过场景背景和雨滴前景的色彩向量夹角约束去除散焦雨滴。若雨滴落在摄像机焦平面上，则被雨滴遮挡的像素的亮度I_r(x，y)如式(14)所示，其中E_r(x，y)是由雨滴本身带来的照度，E_b(x，y)是由背景带来的照度，τ为该像素被雨滴遮挡的时间，T为摄像机曝光时间；若雨滴不在摄像机焦平面上，则其被雨滴遮挡的像素I_r′(x,y)如式(15)所示，此时，雨滴由于散焦的原因，其带来的照度变成E’_r(x，y)，且根据光学原理有，E’_r(x，y)<E_r(x，y)。

I_r(x,y)=τE_r(x,y)+(T-τ)E_b(x,y)    (14)

I′_r(x,y)=τE_r′(x,y)+(T-τ)E_b(x,y)    (15)

实际上I_r(x，y)是不需要通过式(14)来计算的，因为I_r(x，y)就是当前帧被雨滴遮挡的这个像素的色彩向量，可以直接从当前帧图像获得。将E_r(x，y)、E’_r(x，y)、E_b(x，y)、I_r(x，y)和I’_r(x，y)从亮度扩展到RGB三通道色彩向量空间，分别用

表示。根据式(14)、(15)，相对于聚焦雨滴，散焦雨滴在色彩空间的和

向量夹角更小，如图3所示，同时，由于因此散焦雨滴的和

的向量夹角相对于聚焦雨滴的也更小。因此可采用式(16)来区分聚焦雨滴和散焦雨滴：

$\mathrm{\&gamma;}(x,y)={\cos }^{-1}\left(\frac{{\stackrel{\&RightArrow;}{I}}_r(x,y)\&CenterDot;{\stackrel{\&RightArrow;}{I}}_b(x,y)}{\left|{\stackrel{\&RightArrow;}{I}}_r(x,y)\right|\left|{\stackrel{\&RightArrow;}{I}}_b(x,y)\right|}\right)---\left(16\right)$

其中，为当前帧被雨滴遮挡的像素的色彩向量，

是像素未被雨滴遮挡时的色彩向量。

可以直接从当前帧图像获得，

可由通过对前一帧和后一帧的像素色彩向量

取平均获得：

${\stackrel{\&RightArrow;}{I}}_b(x,y)=\frac{{\stackrel{\&RightArrow;}{F}}_{n-1}(x,y)+{\stackrel{\&RightArrow;}{F}}_{n+1}(x,y)}{2}---\left(17\right)$

对于检测到的每个候选雨滴，检查它所包含像素的平均γ值是否小于设定阈值，若是，则判定其为散焦雨滴并移除。

5.利用已知的摄像机参数计算雨滴的尺寸，构建雨滴密度分布函数。为了获得准确的雨滴密度分布函数，摄像机参数应设置成能够观察到最小的雨滴，如图4所示，由图可知，摄像机参数需满足：

z₀≤fd_min    (18)

其中，z₀为焦平面距离，f为摄像机焦距，d_min为期望捕获的最小雨滴的直径，本发明中取为0.3mm。对于落在焦平面上的雨滴，即通过步骤1）~4）最终确定的聚焦雨滴，其直径由式(19)给出：

$d_i=2{\left(\frac{z_0{a}_i}{200\mathrm{fT}}\right)}^2---\left(19\right)$

式(19)利用到了雨滴终极速度v和直径d的关系：

$v=200\sqrt{d/2}---\left(20\right)$

得到每个检测到的聚焦雨滴的直径d_i后，可计算直径d_i落在[d-1/2Δd，d+1/2Δd]区间内的雨滴数目，d由式（20）计算得到，从而构建雨滴密度分布函数N(d)，Δd是单位尺度间隔，通常取1mm。为获得稳定鲁棒的N(d)，可以利用连续多帧图像进行统计。

6.根据雨滴密度分布函数N(d)，利用Gamma雨滴谱分布模型拟合推导降雨量R：

气象学上常用的雨滴密度分布函数——Gamma分布为：

N(d)=N₀d^μe^-ωd    (21)

其中，N₀,μ和ω是Gamma分布的三个参数。根据前面步骤获得的N(d)，可通过参考文献[9]记载的矩方法求得三个参数，从而实现基于Gamma分布模型的雨滴谱拟合：

$M_i=\underset{0}{\overset{8}{\&Integral;}}{d}^{i}N\left(d\right)\mathrm{dd},i=\mathrm{2,3,4}---\left(22\right)$

$\mathrm{\&mu;}=\frac{3M_4{M}_2-4{M_3}^2}{({M_3}^2-M_4{M}_2)}---\left(23\right)$

$\mathrm{\&omega;}=\frac{(\mathrm{\&mu;}+4)M_3}{M_4}---\left(24\right)$

$N_0=\frac{{\mathrm{\&omega;}}^{\mathrm{\&mu;}+4}{M}_3}{\mathrm{\&Gamma;}(\mathrm{\&mu;}+4)}---\left(25\right)$

降雨量R与N(d)的关系如式(26)所示，可利用拟合出来的参数N₀,μ和ω，通过式(26)计算R：

$R=3.6\&times;{10}^{-3}\&times;\frac{\mathrm{\&pi;}}{6}\underset{0}{\overset{\&infin;}{\&Integral;}}N\left(d\right)v\left(d\right)d^3\mathrm{dd}$

$=0.0071N_0\frac{\mathrm{\&Gamma;}(\mathrm{\&mu;}+4.67)}{{\mathrm{\&omega;}}^{\mathrm{\&mu;}+4.67}}---\left(26\right)$

本发明算法的实施效果如图5、6、7所示。利用本发明方法对2012年2月12日11:30-13:10、2012年2月22日13:00-13:37以及2012年3月22日9:42-9:59三个时段的降雨情况进行测量。图5、6、7分别显示了小雨、中雨、大雨三种情况下由本发明算法获得的雨滴谱和Gamma拟合后的雨滴谱，三种情况下利用本发明算法测得的瞬时降雨量分别为1.0mm/h、4.9mm/h、11.0mm/h，可见两个分布较为吻合，说明通过视频图像可以提取出可靠的雨滴谱分布。另外，由本发明方法获得的三个时段的平均降雨量结果分别为0.8mm/h、8.7mm/h、4.7mm/h，与同时段利用量雨筒获得的数据1.1mm/h、6.5mm/h、3.5mm/h较为吻合，验证了本发明方法的有效性。

Claims (7)

1.一种基于视频图像的降雨量测量方法，其特征是对视频摄像机采集的降雨视频，采用对数图像处理LIP方法，用灰色调替代亮度值，建立灰色调约束条件提取候选雨滴；利用主成分分析法PCA构建候选雨滴的倾斜方向的概率密度分布函数，并推导求取该分布函数峰值的Mean-shift方法，通过约束候选雨滴的倾斜方向去除干扰；通过场景背景和雨滴前景的色彩向量夹角约束去除散焦雨滴；利用摄像机参数标定计算每个聚焦雨滴的尺寸，构建雨滴密度分布函数，利用Gamma雨滴谱分布模型拟合推导降雨量，具体如下：

1）检测候选雨滴：利用视频中连续三帧图像：第n－1帧，第n帧和第n＋1帧，在LIP框架下建立灰色调约束，提取得到第n帧图像的候选雨滴；

2）构建雨滴倾斜方向概率分布密度函数：利用PCA方法，提取第n帧图像每个候选雨滴的倾斜方向θ_i以及长a_i、宽b_i，i＝{1,2,…,P_n}，P_n为第n帧图像中检测出的候选雨滴总个数，然后构建雨滴倾斜方向概率密度分布函数f(θ)；

3）确定降雨方向，去除噪声干扰：针对步骤2)中获得的f(θ)，用Mean-shift算法提取它的峰值，所述峰值即为降雨方向，对候选雨滴实施方向约束，不符合降雨方向的判定为噪声；

4）去除散焦雨滴：计算候选雨滴中每个像素的场景背景和雨滴前景的色彩向量夹角γ(x,y)，对每个候选雨滴，检查它所包含像素的平均γ值是否小于设定的阈值，若是则判定为散焦雨滴，否则为聚焦雨滴；

5）构建雨滴密度分布函数：利用视频摄像机参数计算每个聚焦雨滴的尺寸，构建雨滴密度分布函数N(d)；

6）利用步骤5)中的雨滴密度分布函数N(d)，以及Gamma雨滴谱分布模型拟合推导降雨量R。

2.根据权利要求1所述的一种基于视频图像的降雨量测量方法，其特征是步骤1)中的灰色调约束是：

f_n(x,y)<f_n-1(x,y)&f_n(x,y)<f_n+1(x,y)    (1)

|f_n(x，y)Θf_n-1(x，y)|_E＞th & |f_n(x，y)Θf_n+1(x，y)|_E＞th    (2)

其中f_n-1(x，y)、f_n(x，y)、f_n+1(x，y)是n-1、n、n+1连续三帧图像的灰色调函数，th为灰色调阈值，同时满足式(1)、(2)的像素则判定为候选雨滴；

灰色调函数f(x，y)与入射光强度即输入图像F(x，y)的关系为：

$f(x,y)=M(1-\frac{F(x,y)}{F_{\max }})---\left(3\right)$

其中F_max是人眼视觉系统的饱和亮度值，取为255，M是自定义的灰色调函数值域的最大值，取255，由式（3）和连续三帧图像，得到三帧图像各自的灰色调函数；

式(2)中Θ代表LIP框架下灰色调空间的减法操作，两个灰色调函数f(x，y)和g(x，y)的减法定义为：

$f(x,y)\mathrm{\&Theta;g}(x,y)=M\frac{f(x,y)-g(x,y)}{M-g(x,y)}---\left(4\right)$

式(2)中“||_E”表示LIP框架下灰色调空间的取模操作，定义为：

${\left|f(x,y)\right|}_E=\left\{\begin{array}{ccc}f(x,y);& \mathrm{if}& f(x,y)\&GreaterEqual;0\\ \mathrm{\&Theta;f}(x,y);& \mathrm{if}& f(x,y)<0\end{array}\right.---\left(5\right).$

3.根据权利要求1所述的一种基于视频图像的降雨量测量方法，其特征是步骤2)中提取候选雨滴的倾斜方向θ_i以及长a_i、宽b_i的方法是：

令N_i表示第i个雨滴占有的像素个数，X_i为N_i*2的矩阵：

$X_i=\left[\begin{array}{cc}{x}_{i1}& y_{i1}\\ x_{i2}& y_{i2}\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ x_{{\mathrm{iN}}_i}& y_{{\mathrm{iN}}_i}\end{array}\right]---\left(6\right)$

其中（x_ik，y_ik），是该雨滴中第k个像素相对于该雨滴质心的坐标，k＝1,2,…N_i，设向量u₁ ⁱ=[u_1x ⁱ,u_1y ⁱ]^T,u₂ ⁱ＝[u_2x ⁱ,u_2y ⁱ]^T是矩阵X_i ^TX_i的两个特征向量，分别对应于特征值λ₁ ⁱ,λ₂ ⁱ，λ₁ ⁱ＞λ₂ ⁱ，u₁ ⁱ即表示第一主成分方向，则雨滴的倾斜方向θ_i为：

${\mathrm{\&theta;}}_i={\tan }^{-1}\left(\frac{{u_{1y}}^i}{{u_{1x}}^i}\right)---\left(7\right)$

长a_i、宽b_i为：

$a_i=4\sqrt{{{\mathrm{\&lambda;}}_1}^i},$ $b_i=4\sqrt{{{\mathrm{\&lambda;}}_2}^i}---\left(8\right)$

雨滴倾斜方向概率密度分布函数f(θ)的构建为：

$f\left(\mathrm{\&theta;}\right)=\underset{i=1}{\overset{P_n}{\mathrm{\&Sigma;}}}\frac{N_i}{\sqrt{2\mathrm{\&pi;}}{h}_i}k\left({\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2\right)---\left(9\right)$

其中k(X)为核函数， $k\left(X\right)=e^{-\frac{1}{2}X},$ 对应式（9）则 $X={\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2,$ 核函数带宽h_i用于描述θ_i的不确定性，计算方法为：

$h_i=C\frac{b_i}{a_i}+H---\left(10\right)$

C是将视频图像的长宽比值转换至核函数带宽时的尺度缩放，H表示核函数基本带宽，C、H均为常数。

4.根据权利要求1所述的一种基于视频图像的降雨量测量方法，其特征是步骤3)中求取f(θ)峰值的Mean-shift方法为：

步骤1:设定迭代初始值θ₀=0;

步骤2:令θ₁＝m(θ₀)+θ₀;

步骤3:若|θ₁-θ₀|<ε，令降雨方向θ_main=θ₁;否则令θ₀=θ₁并回到步骤2进行循环判断；

其中，

$m\left(\mathrm{\&theta;}\right)=\frac{\underset{i=1}{\overset{P_n}{\mathrm{\&Sigma;}}}\frac{2N_i}{\sqrt{2\mathrm{\&pi;}}{h_i}^3}{\mathrm{\&theta;}}_{i}g\left({\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2\right)}{\underset{i=1}{\overset{P_n}{\mathrm{\&Sigma;}}}\frac{2N_i}{\sqrt{2\mathrm{\&pi;}}{h_i}^3}g\left({\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2\right)}-\mathrm{\&theta;}---\left(12\right)$

g(X)=-k′(X)

-k'(X)是核函数k(X)的导数再取负，核函数 $k\left(X\right)=e^{-\frac{1}{2}X},$ 对应式（12）有 $X=$ ${\left|\right|\frac{\mathrm{\&theta;}-{\mathrm{\&theta;}}_i}{h_i}\left|\right|}^2;$

根据降雨方向θ_main，对候选雨滴实施方向约束：

如果第i个雨滴的方向θ_i相对于降雨方向θ_main的差值大于设定的阈值θ_T，则判定为干扰噪声：

|θ_i-θ_main|>θ_T    (13)

阈值θ_T设为10°。

5.根据权利要求1所述的一种基于视频图像的降雨量测量方法，其特征是步骤4)中色彩向量夹角γ(x,y)计算方法为：

$\mathrm{\&gamma;}(x,y)={\cos }^{-1}\left(\frac{{\stackrel{\&RightArrow;}{I}}_r(x,y)\&CenterDot;{\stackrel{\&RightArrow;}{I}}_b(x,y)}{\left|{\stackrel{\&RightArrow;}{I}}_r(x,y)\right|\left|{\stackrel{\&RightArrow;}{I}}_b(x,y)\right|}\right)---\left(16\right)$

其中，

为当前帧被雨滴遮挡的像素的色彩向量，

是像素未被雨滴遮挡时的色彩向量，直接从当前帧图像获得，

通过对前一帧和后一帧的像素色彩向量

取平均获得：

${\stackrel{\&RightArrow;}{I}}_b(x,y)=\frac{{\stackrel{\&RightArrow;}{F}}_{n-1}(x,y)+{\stackrel{\&RightArrow;}{F}}_{n+1}(x,y)}{2}---\left(17\right).$

6.根据权利要求1所述的一种基于视频图像的降雨量测量方法，其特征是摄像机参数满足z₀≤fd_min，f为摄像机焦距，z₀为焦平面距离，d_min为期望捕获的最小雨滴的直径，取为0.3mm；

根据雨滴的终极速度v和其直径d的关系：

步骤5)中雨滴密度分布函数N(d)构建为：

对于落在焦平面上的雨滴，即聚焦雨滴，其直径d_i为：

$d_i=2{\left(\frac{z_0{a}_i}{200\mathrm{fT}}\right)}^2---\left(19\right)$

T为摄像机曝光时间，

得到每个检测到的聚焦雨滴的直径d_i后，计算直径d_i落在[d-1/2Δd，d+1/2Δd]区间内的雨滴数目，从而构建雨滴密度分布函数N(d)，Δd是单位尺度间隔，取1mm，利用连续多帧图像进行统计，以获得稳定鲁棒的N(d)。

7.根据权利要求1所述的一种基于视频图像的降雨量测量方法，其特征是步骤6)中降雨量计算方法为根据雨滴密度分布函数N(d)，利用Gamma雨滴谱分布模型拟合推导降雨量R：

$R=0.0071N_0\frac{\mathrm{\&Gamma;}(\mathrm{\&mu;}+4.67)}{{\mathrm{\&omega;}}^{\mathrm{\&mu;}+4.67}}---\left(26\right)$

其中，N₀,μ和ω是Gamma雨滴谱分布模型的三个参数，根据雨滴密度分布函数N(d)求得：

$\mathrm{\&mu;}=\frac{3M_4{M}_2-4{M_3}^2}{({M_3}^2-M_4{M}_2)}---\left(23\right)$

$\mathrm{\&omega;}=\frac{(\mathrm{\&mu;}+4)M_3}{M_4}---\left(24\right)$

$N_0=\frac{{\mathrm{\&omega;}}^{\mathrm{\&mu;}+4}{M}_3}{\mathrm{\&Gamma;}(\mathrm{\&mu;}+4)}---\left(25\right)$

$M_i=\underset{0}{\overset{8}{\&Integral;}}{d}^{i}N\left(d\right)\mathrm{dd},i=\mathrm{2,3,4}---\left(22\right).$

Images