激光雷达中基于组合纹理的低空风切变识别算法

蒋立辉，高 浩，庄子波，熊兴隆

(1.民航气象研究所中国民航大学，天津300300;2.智能信号与图像处理天津市重点实验室中国民航大学，天津300300)

1 引言

风切变表现为空间两点之间风的矢量差，即在统一高度或不同高度短距离内风向和风速的变化，飞机在风向风速不断发生变化的环境中飞行，相应地就要发生突然性的空速变化，空速变化引起了升力变化，升力变化又引起了飞行高度的变化。如果遇到的是空速突然减小，而飞行员又未能立即采取措施，飞机就要掉高度，以至发生事故。通过对严重事故的分析确认后，低空风切变是引起这些飞机失事的主要原因。

目前对于低空风切变的识别这方面的课题研究还不是很成熟，如何有效地识别风切变日益成为保障民航安全运输的一个重要课题，国内外在此方面也取得了一定的进展。Potts认为下击暴流最可靠的特征是地面辐散，在多普勒速度图上表现为在雷达径向上的“牛眼”回波［1］。1987年Merritt根据下击暴流到达地面后向四周辐散的特性，提出了基于切变段提取的下击暴流探测算法［2］。Hermes在研究自动识别、追踪阵风锋的方法中指出，当阵风锋与雷达径向不垂直时，辐合区在多普勒雷达上显示为速度沿径向减小的区域［3］。

近年来，由于激光雷达在探测低空风切变中具有精度高、分辨率高等特点，基于二维图像的方法对风切变进行识别成为研究的新思路，高志光采用了Hu矩和灰度共生矩阵的串联组合特征对风切变进行了识别［4］。然而Hu矩和灰度共生矩阵都只能刻画图像的全局特征，在同类风切变风向变化较大，尤其是风向与雷达径向垂直时并不能得到理想的识别效果。

文中对激光雷达探测的二维风切变图像进行了纹理分析，将风切变的纹理特征分为局部特征和全局特征，局部的纹理特征可以提取到风场局部风速变化较为微弱的信息，全局的纹理特征代表风场全局风速的变化，而对局部风速变化小的区域较为不敏感，二者都包含有代表风场类别的信息，文中采用典型相关分析对二者进行融合，最后通过最近邻分类器进行匹配识别。识别的流程如图1所示。

图1 风切变识别的流程图Fig.1 Flow chart of the proposed method

2 样本的仿真

由于真实的风切变数据难以获取，本节首先基于大气流场理论，运用Fluent流场软件对微下击暴流、低空急流、海陆风进行了仿真模拟。进而通过多普勒激光雷达速度方位(VAD)扫描方式得到三种风切变的仿真激光雷达图像。

2.1 风切变风场的建模

低空风切变是指600 m以下空气层中风向和风速突然改变的现象。三种低空风切变的风场特性可概括如下:

下击暴流是雷暴云中的强烈下沉气流，在冲击地面后向各个方向辐散。

低空急流是对流层下层的一个准管状的狭长强风带，在强风带中心风速最大，上下和两侧逐渐减小。

海陆风切变是指海风从海面吹向大陆，与大陆上相反方向的背景风汇聚而形成的风切变，在海风锋两侧风向相反。

根据以上三种低空风切变的特性，分别选用“壁面射流”模型和“平面自由射流”模型对风场进行数学描述，选用相应的湍流模型和边界条件，具体参见文献［5］，通过Fluent数值模拟仿真后便可得到三种低空风切变的三维风场，三种风场的尺寸均为8000 m×8000 m×1000m，如图2所示。

图2 三种风切变的三维流场仿真图Fig.2 The three － dimensional field simulation figure

2.2 仿真激光雷达图像的生成

多普勒激光雷达速度方位扫描(VAD)方式［6］如图3所示。

图3 雷达圆锥扫描示意图Fig.3 The scanning way of lidar

雷达波束在一个固定仰角α上，沿方向角θ进行扫描。设去向的径向速度为正，来向的径向速度为负，雷达位于坐标原点。其中，Vk为水平风，Vf为垂直风，θ0为水平风和X轴的夹角。

由图4可知，径向速度v:

式(1)描述了雷达径向风速与三维实际风场各个方向上风速分量的关系。这样就可以由仿真得到的三维风场数据得到基于激光雷达的变化风场探测数据。待测风场的激光雷达PPI图像如图4所示，周围的蓝色区域代表无有效地风场数据，图像的高度代表雷达径向距离库，水平宽度代表方位角，单位库长为100 m，方位角分辨率为1°。

图4 三种风切变的激光雷达PPI图像Fig.4 The lidar PPI image of three kinds of wind shear

3 组合纹理特征提取

3.1 局部纹理特征

局部二进制模式(Local Binary Pattern，LBP)［7］是一种基于灰度描述图像纹理特征的不相关算子，它通过对图像任意一点与其周围点的灰度值大小关系来表征图像的局部纹理特征，具有对灰度值的线性变化不敏感的特性。

LBP算子通常由参数(P，R)来表示，其中P表示邻域内包含的像素个数，R表示邻域半径，对于不同的(P，R)值，LBP算子也不相同，图5为两种不同LBP算子。

图5 两种不同的LBP算子Fig.5 Two different LBP operator

对于任意的LBP算子，它的编码公式为:

式(2)中，gc为中心像素;gi(i=0，…，P －1)为以中心像素为圆心;R为半径的圆周上的P个像素值。LBP算子的计算过程如图6所示。

图6 LBP算子的计算过程Fig.6 The calculation process of LBP operator

为了得到旋转不变的LBP算子，文献［8］对原始的LBP算子做了改进，即不断旋转圆形邻域得到一系列的LBP值，取其最小值作为该邻域的LBP值，计算公式如下:

式(3)中，LBPriP，R表示旋转不变得 LBP 算子，ROR(x，i)为旋转函数，表示将x循环右移i位。风切变图像经过LBP算子运算后，可得到对应的LBP编码图像，如图7所示。对得到的LBP编码做直方图统计，将归一化后统计直方图做为最终的局部特征向量。

图7 三种风切变的LBP编码图像Fig.7 The LBP coded image of three kinds of wind shear

3.2 全局纹理特征

灰度－梯度共生矩阵(Gray－Gradinent Co－occurrence Matrix，GGCM)是利用整幅图像的灰度和梯度综合信息来提取纹理特征。灰度－梯度共生矩阵的元素H(i，j)定义为在归一化的灰度图像F(m，n)和归一化的梯度图像G(m，n)中共同具有灰度为i和梯度为j的总像素点数。以图像的总像素点数归一化后得到概率为p(i，j)。从灰度－梯度共生矩阵中，可以提取15个二次特征参量［9］。

在雷达风切变图像中，灰度代表径向风速值，梯度代表径向风速的切变值，为了克服传统的梯度算子对噪声敏感的不足，文中采用最小二乘拟合法分别计算雷达的径向速度的径向切变和径向速度的切向切变(除以r是为了要归一化)，并将径向方位组合切变G做为最终的梯度值，计算公式如下:

式(4)中，vi为第i个点的径向速度;ri为相应点距雷达的距离;n为参与计算的同一径向上的像素点数。

式(5)中，vj为第j个点的径向速度;θj为相应地方位角;r为距离雷达的距离;m为参与计算的相同径向距离上像素点数。

由最小二乘求得的三种风切变图像的切变值图像如图8所示。

图8 三种风切变的切变值图像Fig.8 The shear value images obtained by LS

3.3 特征融合

典型相关分析(Canonical Correlation Analysis，CCA)是一种提取两个特征集间线性相关的强有力的工具。通过局部特征和全局特征的提取，我们得到两组特征集:X=［x1，…，xn］，Y=［y1，…，yn］，其中 xi∈Rp，yi∈Rq，i=1，2，…，n。n 为样本的个数。CCA的目标可看为分别对特征集X和特征集Y寻找两组基向量wx∈Rp和wy∈Rq。使得随机变量u=wTxx和v=wTyy之间的相关性达到最大，即相关系数ρ(u，v)最大。CCA可以通过下面的式子描述:

式(7)是一种最优化问题，可以通过Lagrange乘子法求得最优解(wxi，wyi)，i=1，…，d(d≤min(p，q))，具体的求解方法参见文献［10］，这样便求出了两组投影矩阵 Wx=［wx1，…，wxd］，Wy=［wy1，…，wyd］，融合后的特征可以定义为:

其中，ti即为全局纹理特征和局部纹理特征融合后的组合纹理特征(Global－Local Texture Features，GLTF)。

4 应用最近邻分类器进行分类识别

从激光雷达风切变图像中提取到LBP特征，灰度－梯度共生矩阵特征并构成融合特征后，可以通过最近邻分类器进行分类识别。

最近邻分类器是按距离来分类的，这种分类器较简单且花费时间很小，其分类的思想为:如果未知样本f的特征和已知类别样本fk的特征之间的欧氏距离最小，则判断f与fk同类，假设有c个类别w1，w2，…，wc，每类有表明类别的样本 Ni个，则可以规定wi类的判别函数为:

决策规则为:若gj(f)=(f)，i=1，2，…，c，则 f∈wj。

5 结果分析

文中采用自建立的雷达扫描图像样本库对常见的三种低空风切变进行分类识别，每种激光雷达风切变图像80幅，分别为待测风场位于激光雷达不同方位下的扫描图像，总计240张图像，每幅图像的尺寸分别代表径向距离和方位角。识别过程采用交叉验证的策略，以避免人为选择样本所导致的识别结果的偶然性，即随机选取20幅作为训练样本，剩下的作为测试样本，重复进行50次。

为了突出文中算法的优势，将其与文献［4］中采用的Hu+GLCM特征、LBP特征、GGCM特征所得的识别率进行对比。Hu+GLCM包括Hu矩的7个形状特征和GLCM的4个纹理特征;LBP算法采用p=8，R=1的旋转不变LBP算子，包含36个特征;GGCM包括小梯度优势、大梯度优势、灰度分布的不均匀性、梯度平均、灰度熵、混合熵、惯性、逆差矩等15个全局统计特征;GLTF由 LBP特征和GGCM特征经CCA融合而成，特征维数为15。实验结果如图9和表1所示。

图9 四种算法的识别率曲线Fig.9 The recognition rate curve of four algorithms

表1 四种算法识别率的对比Tab.1 Recognition rate comparison of the four algorithms %

图9给出了50次分别应用Hu+GLCM特征，LBP特征，GGCM特征和GLTF特征进行风切变识别时的识别曲线，表1给给出了50次识别的平均识别率，S1、S2、S3分别代表微下击暴流、低空急流、海陆风三种风切变，“±”后面的数字是识别率的标准差。

从表1中可以看出，对三种低空风切变的识别中，Hu+GLCM算法、LBP算法和GGCM算法的平均识别率只有83.31%，93.52%和 92.53%，运用文中所述的组合纹理特征进行识别，平均识别率达到99.02%，识别效果得到了很大改善。

6 结论

本文提出了基于组合纹理的低空风切变算法，从激光雷达探测到的低空风切变图像中提取LBP特征和GGCM特征，通过CCA融合得到组合纹理特征，LBP特征反应图像局部的纹理，可以提取到风场局部变化较为微弱的信息，GGCM特征反应图像的全局纹理，代表风场全局风速的变化。通过对三种低空风切变的识别，该算法取得了较高的识别率，而且特征维数较少，达到了较好的分类识别效果。

［1］ Potts R.Microburst precursors observed with Doppler radar［C］//Conference on Radar Meteorology，24th，Tallahassee，FL，1989:158 －162.

［2］ Merritt M W.Automated detection of microburst windshear for terminal doppler weather radar［C］//1987 Cambridge Symposium.International Society for Optics and Photonics，1988:61 －69.

［3］ Hermes L G，Witt A，Smith SD，et al.The gust－ front detection and wind－shift algorithms for the terminal doppler weather radar system［J］.Journal of Atmospheric and O-ceanic Technology，1993，10(5):693 －709.

［4］ Jiang Lihui，Gao Zhiguang，Xiong Xinglong，et al.The study on type recognition of low altitude wind shear based on lidar image processing.［J］.Infrared and Laser Engineering，2012，41(12):3410 －3415.(in Chinese)蒋立辉，高志光，熊兴隆.基于激光雷达图像处理的低空风切变类型识别研究［J］.红外与激光工程，2012，41(12):3410－3415.

［5］ Liu Shuo.The research on numerical simulation of high rise building outdoor flow field［D］.Harbin:Harbin Institute of Technology，2007.(in Chinese)刘朔.高层建筑室外气流场的数值模拟研究［D］.哈尔滨:哈尔滨工业大学，2007.

［6］ Wang Bangxin，Shen Fahua，Sun Dongsong，et al.Beam scanning and wind field measurement of direct－detection doppler lidar［J］.Infrared and Laser Engineering，2007，36(1):69 －72.(in Chinese)王邦新，沈法华，孙东松，等.直接探测多普勒激光雷达的光束扫描和风场测量［J］.红外与激光工程，2007，36(1):69 －72.

［7］ Ojala T Pietikainen M，Harwood D.A comparative study of texture measures with classification based on feature distribution［J］.Pattern Recognition，1996，29(1):51－59.

［8］ Ojala T Pietikainen M，Maenpaa T.Multiresolution gray scale and rotation invariant texture analysis with local binary patterns［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2002，24(7):971 －987.

［9］ Xia Deshen，Jin Sheng，Wang Jian.Fractal dimension and GGCM meteorology cloud pictures recognition［J］.Journal of Nanjing University of Science and Technology，1999，23(4):289 －292.(in Chinese)夏德深，金盛，王建.基于分数维与灰度梯度共生矩阵的气象云图识别［J］.南京理工大学学报，1999，23(4):289－292.

［10］ DR Hardoon，SSzedmak，JShawe － Taylor.Canonical correlation analysis:an overview with application to learning methods［J］.Neural Computation，2004，16(12):2639－2664.