一种基于小波变换和支持向量机的涡结构识别方法

摘要

本发明涉及一种基于支持向量机(SVM，SupportVectorMachine)的涡结构识别方法，在特征提取方面，提取折射率场数据经小波变换后系数矩阵的统计量作为涡结构特征量，与传统的直接利用小波变换后系数矩阵相比，剔出了大量的冗余信息，降低了计算量。在分类器设计方面，提出了基于支持向量机的分类方法，利用结构风险最小原则使分类面不仅能将涡结构正确分开，而且使分类间隔最大，比传统的经验风险最小原则分类方法从原理上降低了误识率。本发明能更加准确地表征和区分湍流涡结构，为导弹、飞机等机载光学设备的光学效应精确建模奠定了基础。

说明书

技术领域

本发明涉及一种基于小波变换和支持向量机的涡结构识别方法，可用于分析大气中高速飞行的导弹和飞机的光学头罩与来流之间形成的气动光学效应，为图像校正、控制湍流提供基础。

背景技术

高速成像制导导弹和飞机等在大气层内高速飞行时，光学头罩周围会形成高速复杂流场，致使目标光波通过高速流场时光波相位发生改变，目标图像产生偏移、模糊、抖动和能量衰减，这种现象称为气动光学效应。气动光学效应的研究重点之一是畸变波前同湍流涡动力学特性之间的关系。按照湍流的涡旋学说，湍流的脉动与混合主要由大尺度的涡结构造成，具有高度的各向异性，且其流动、演变具有重复性和可预测性；小尺度涡几乎是各向同性的，而且不同流动的小尺度涡有许多共性。大尺度涡造成了大部分的密度脉动，在光学传输中造成了大部分的相位脉动，是影响气动光学传输效应的主要因素。因此要按照不同的涡结构尺度建立精确的光学传输模型，就必须进行涡结构的识别。

涡结构识别的关键在于特征提取和分类器设计。传统的涡结构识别方法，在特征提取方面通常是依据小波变换具有多分辨分析特性，直接进行小波分解，利用小波系数矩阵作为特征量使得特征冗余多，计算量大；在分类器设计方面传统方法通常直接根据经验确定小波分解水平来区分大、小尺度涡结构，识别率低。

发明内容

本发明的目的是：克服传统的涡结构识别方法的不足，提出一种基于小波变换和支持向量机的涡结构识别方法。

本发明的技术方案是：把机载光学窗口处湍流折射率场经小波变换后得到各频带的系数矩阵，以其统计量作为涡结构的特征量进行识别。提取的特征有小波系数矩阵的均值、方差、能量、信息熵、对角线惯性矩、绝对均匀性和聚类重要性。再利用支持向量机方法建立分类器，根据结构风险最小原则使分类面不仅能将涡结构正确分开，而且分类间隔最大。

具体步骤如下：

(1)利用小波变换对描述涡结构的折射率进行小波分解，得到不同尺度上的小波系数矩阵。

(2)对小波系数矩阵进行特征提取，小波系数矩阵的统计量作为涡结构的特征值。

(3)根据所提取特征值，利用支持向量机的方法建立最佳分类判别函数，其建立准则采用利用风险最小化原则，精确识别大小涡结构。

(4)利用步骤(3)建立的最佳分类判别函数，重复步骤(1)、(2)对涡结构样本进行涡结构识别。

(5)利用小波反变换进行大、小尺度涡重建。

本发明的基本原理是：导弹和飞机周围的绕流流场中存在不同尺度的涡结构，而具有“数字显微镜”之称的小波分析具有良好时、频特性以及多尺度多分辨分析特性，是研究精细结构的理想工具，可以很好地反映涡结构。传统的方法直接利用小波变换的系数作为涡结构的特征量，这样造成特征冗余大，存储和计算资源的浪费。表示涡结构的矩阵是一个随机变量，用统计量来描述更为科学，因此本发明将折射率场经小波变换后得到各频带的系数矩阵的统计量作为涡结构特征量，并以此进行识别，提取的特征有小波系数矩阵的均值、方差、能量、信息熵、对角线惯性矩、绝对均匀性和聚类重要性，这些涡结构的特征就构成了涡结构模式的样本，简称涡结构样本，记为x_i；涡结构的类别记为y_i∈{-1，+1}，i＝1，2...l，l为样本个数。已知涡结构类型的样本包括大尺度涡结构样本和小尺度涡结构样本，大尺度涡结构样本由流场的折射率求得，流场的折射率是通过利用大涡模拟的方法求解Navier-Stokes方程(简写N-S方程)，得到密度场数据再经过Gladstone-Dale变换(简称GD变换)求得，小尺度涡结构的样本通过N-S方程的精确数值解求得，N-S方程的精确数值解中既包含大尺度涡模式又包含小尺度涡模式，从N-S方程的精确数值解中减去大尺度涡模式便得到小尺度涡模式的密度场描述，在经过GD变换便可得到小尺度涡结构。但由于建立N-S方程和求解精确数值运算量巨大，需要多节点计算机同时并行计算多天，因此本发明提出基于小波变换和支持向量机方法识别方法，大大提高计算效率。利用N-S方程精确解出的涡结构作为已知类别的样本建立最佳分类判别函数，其原理是因为涡结构样本属于二维线性不可分问题，通过核函数进行非线性变换，将低维不可分的情况转化为高维可分空间，即通过函数映射x→z＝(x)，把对利用特征x的二维线性不可分问题转化为高维线性可分问题，因此定义高维线性判别函数f＝sgn[w^T·z+b]，w、b为权值向量及偏置，对特征x识别问题转换为对特征z的识别问题。如图3所示，表示了涡结构分类问题的基本原理，支持向量机的主要思想是利用核函数将输入向量x映射到一个高维特征空间z，并在该空间内构造一个最优分类超平面来逼近分类函数F＝w^T·z+b＝0。因为原x特征空间为非线性空间，映射到z线性空间，线性空间的线性函数表达式为w^T·z+b＝0，最优分类超平面的构造最终归结为在原空间上求解一个凸二次规划问题。图3中○和□分别代表两类涡结构样本，F与F₁之间的间隔为，F₁与F₂之间的距离为分类间隔，若使分类间隔最大，即要最小。因为样本不完全线性可分，可求解$\frac{1}{2}{w}^{T}w+h\underset{i=1}{\overset{l}{\Sigma }}{\xi }_i$的最小值，其中h为惩罚参数，h越大表示对错误分类的惩罚越大；ξ为分类误差。采用拉格朗日乘子法求解这个具有线性约束的二次规划问题，目标函数表示为$\underset{\lambda }{\max }\underset{w,b}{\min }\{L(w,b,\lambda )=\frac{1}{2}{\left|\right|w\left|\right|}^2+h\underset{i=1}{\overset{l}{\Sigma }}{\xi }_i-\underset{i=1}{\overset{l}{\Sigma }}{\lambda }_i[y_i(w^T·z_i+b)-1]\},$进一步转化得$w=\underset{i=1}{\overset{l}{\Sigma }}{\lambda }_i{y}_i{z}_i,$由此构成的最佳分类判别函数为$f=\mathrm{sgn}[\underset{i=1}{\overset{l}{\Sigma }}{\lambda }_i{y}_i(z_i^T·z)+b]=\mathrm{sgn}[\underset{i=1}{\overset{N_{\mathrm{NSV}}}{\Sigma }}{\lambda }_i{y}_i(z_i^T·z)+b].$虽然将特征x映射高维特征空间，维数有所提高，但也只考虑在高维特征空间的点积运算，故不必明确知道映射函数，构造内积核函数为$K(x_i·x)=z_i^T·z.$就可得到非线性支持向量机的判决函数$f=\mathrm{sgn}[\underset{i=1}{\overset{N_{\mathrm{NSV}}}{\Sigma }}{\lambda }_i{y}_{i}K(x_i,x)+b],$N_NSV表示支持向量的个数，从最佳分类判别函数可以看出，i为样本序号，l表示所有样本的个数，即所有特征值λ_i所对应的样本。而经过推导只有不为零的λ_i才对分类有贡献，它所对应的样本就是支持向量x_NSV.在计算分类函数时，不需计算l个样本，只需计算N_NSV样本即可，也就是支持向量x_NSV决定了分类超平面，而其他的样本(非支持向量)并不影响分类。因此支持向量机分类判决函数是一种用支持向量集合对所有样本数据的压缩表示，降低了维数，这是其他机器学习技术(例如神经网络技术)所不具备的优点，最优分类面就是使分类面不但能将两类正确分开(训练错误率为0)，而且使分类间隔最大，即分类器具有最好的泛化能力。

本发明与现有技术相比的优点在于：

(1)在特征提取方面，把折射率场数据经小波变换后系数矩阵的统计量作为涡结构特征量，包括信息熵、绝对均匀性和聚类重要性等七个特征，比传统的直接利用小波系数矩阵作为特征量，不仅剔出了大量的冗余信息，降低了计算量。

(2)在分类器设计方面，提出了基于支持向量机的分类方法，利用结构风险最小原则使分类面不仅能将涡结构正确分开，训练错误率为0，而且使分类间隔最大，比传统的经验风险最小原则分类方法更具有泛化能力。本发明能更加准确地表征和区分湍流涡结构，为导弹、飞机等机载光学设备的光学效应精确建模奠定了基础。

附图说明

图1为本发明的涡结构识别方法流程图。

图2为塔式小波变换图，D为低频部分，H为高频，下标1，2表示图像的第一级和第二级分解。塔式小波变换是不断地对每一级的低频子带进行滤波和采样处理，保留每级中的高频子带部分。

图3为最优分类超平面。○和□分别代表两类涡结构样本，F为分类超平面，F₁和F₂分别为各类中距离分类超平面最近的涡结构样本且平行于分类超平面的平面，它们之间的距离即记为分类间隔。

具体实施方式

如图1～2所示，本发明的具体方法如下：

1利用Daubechines小波变换对折射率进行小波分解，得到不同尺度上的小波系数矩阵。

原始数据为光学窗口附近的折射率场数据。小波变换方式采用塔式信号分解方式，如图2所示，D为低频部分，H为高频，下标1，2表示图像的第一级和第二级分解。塔式小波变换是不断地对每一级的低频子带进行滤波和下采样处理，保留每级中的高频子带部分。小波函数采用Daubechies小波，小波分解水平N根据具体问题进行选择，与所要求的细节有关，N越大分解水平越高，得到的细节越多。假定对于某个分解水平M(0≤M≤N)，N为最高的分解水平，分解从水平N到N-M，由综合考虑本实施例分解的折射率数据时N，M都取为3，即分解水平分别为3，2，1，0。

小波分解过程为根据Daubechies小波函数和尺度函数确定{a_n，m}，{b_n，m}-L₁≤n≤0，-L₂≤m≤0，{a_n，m}，{b_n，m}是无限序列的截断，这两个序列可通过现有的小波分析方法得到。L₁，L₂相当于滤波窗口的大小，L₁＝L₂＝2N-1。对已知的折射率场数据[c_k；n，m]，k＝N-1，...N-M；-L₃^(k)≤n≤L₄^(k)；-L₅^(k)m≤L₆^(k)，其中k为分解的水平，L₃，L₄分别表示数据的起始和终止的行数；L₅，L₆分别表示数据的起始和终止的列数。本实施例L₃^(N)＝L₄^(N)＝L₆^(N)＝128。令$L_3^{\left(k\right)}=L_3^{(k+1)}/2,$$L_4^{\left(k\right)}=(L_1+L_4^{(k+1)})/2;$$L_5^{\left(k\right)}=L_5^{(k+1)}/2,$$L_6^{\left(k\right)}=(L_2+L_6^{(k+1)})/2.$小波分解的系数表示为

$c_{k;n,m}=\underset{l=-L_1+n}{\overset{n}{\Sigma }}\underset{j=-L_2+m}{\overset{m}{\Sigma }}{a}_{l-2n,j-2m}{c}_{k+1;l,j}$

${d^p}_{k;n,m}=\underset{l=-L_1+n}{\overset{n}{\Sigma }}\underset{j=-L_2+m}{\overset{m}{\Sigma }}{b^p}_{l-2n,j-2m}{c}_{k+1;l,j};p=\mathrm{1,2,3}---\left(1\right)$

2对小波系数矩阵进行特征提取。

因为大尺度涡是会产生较大的折射率变化，是气动光学效应的主要原因，故只计算低频系数c_k；n，m的特征值，包括均值、方差、能量、信息熵、对角线惯性矩、绝对均匀性和聚类重要性作为涡结构的特征。

1)均值

${x_1}^{\left(k\right)}\text{=}\underset{n=L_3^k}{\overset{L_4^{\left(k\right)}}{\Sigma }}\underset{m=L_5^k}{\overset{L_6^{\left(k\right)}}{\Sigma }}n·m·c_{k;n,m}---\left(2\right)$

2)方差

${x_2}^{\left(k\right)}=\frac{1}{n·m}\underset{n=L_3^k}{\overset{L_4^{\left(k\right)}}{\Sigma }}\underset{m=L_5^k}{\overset{L_6^{\left(k\right)}}{\Sigma }}{(c_{k;n,m}-{x_1}^{\left(k\right)})}^2---\left(3\right)$

3)能量

${x_3}^{\left(k\right)}=\underset{n=L_3^k}{\overset{L_4^{\left(k\right)}}{\Sigma }}\underset{m=L_5^k}{\overset{L_6^{\left(k\right)}}{\Sigma }}{c^2}_{k;n,m}---\left(4\right)$

4)信息熵

${x_4}^{\left(k\right)}=\underset{m=L_3^k}{\overset{L_4^{\left(k\right)}}{\Sigma }}\underset{m=L_5^k}{\overset{L_6^{\left(k\right)}}{\Sigma }}{c^2}_{k;n,m}{\log }_{10}{c}_{k;n,m}---\left(5\right)$

5)对角线的惯性矩

${x_5}^{\left(k\right)}=\underset{n=L_3^k}{\overset{L_4^{\left(k\right)}}{\Sigma }}\underset{m=L_5^k}{\overset{L_6^{\left(k\right)}}{\Sigma }}{(m-n)}^2{c}_{k;n,m}---\left(6\right)$

6)绝对均匀性

${x_6}^{\left(k\right)}=\underset{n=L_3^k}{\overset{L_4^{\left(k\right)}}{\Sigma }}\underset{m=L_5^k}{\overset{L_6^{\left(k\right)}}{\Sigma }}\frac{1}{1+{(n-m)}^2}{c}_{k;n,m}---\left(7\right)$

7)聚类重要性

${x_7}^{\left(k\right)}=\underset{n=L_3^k}{\overset{L_4^{\left(k\right)}}{\Sigma }}\underset{m=L_5^k}{\overset{L_6^{\left(k\right)}}{\Sigma }}{(n+m-{{2x}_1}^{\left(k\right)})}^4{c}_{k;n,m}---\left(8\right)$

3根据所提取的特征和已知涡结构样本类别，利用支持向量机的方法建立最佳分类判别函数。

因为折射率表示的涡结构x是非线性的，要是使其线性可分，就必须对x作非线性变换z＝φ(x)，使其成为线性可分，则在映射后z空间中求取权值向量w及偏置b，使其满足线性可分。涡结构识别中输入x_i映射成z_i，输出模式表示为y_i，分别表示大尺度涡结构和小尺度涡结构{z_i，y_i}，i＝1，2...l，z_i∈R^d，yⁱ∈{-1，+1}，z_i，y_i分别表示第i个样本所对应的特征和样本类别，l为样本总数，z_i∈R^d表示z_i属于d维实数空间，y_i样本类别有两类分别表示大涡结构和小涡结构。。而非线性变换的映射函数并不需要知道显式形式，只需计算核函数。定义大涡和小涡分别表示为：

若y_i＝1，   (w^T·z_i)+b≥1

若y_i＝-1，  (w^T·z_i)+b≤-1

以上二式合写为

y_i[(w^T·z_i)+b]-1≥0    (9)

最优分界面方程为：

F＝w^T·z+b＝0          (10)

由此可见求解最优分界面F应在式(9)的约束下，即使两类涡结构样本满足y_i(w^T·z_i+b)≥1，i＝1，2...l，这样分类间隔为，要使分类间隔最大，即要最小，因为样本不完全线性可分，可求解$\frac{1}{2}{w}^{T}w+h\underset{i=1}{\overset{l}{\Sigma }}{\xi }_i$的最小值，采用拉格朗日乘子法求解这个具有线性约束的二次规划问题，目标函数表示为：

$\left(\begin{array}{c}\underset{\lambda }{\max }\underset{w,b}{\min }\{L(w,b,\lambda )=\frac{1}{2}{\left|\right|w\left|\right|}^2+h\underset{i=1}{\overset{l}{\Sigma }}{\xi }_i-\underset{i=1}{\overset{l}{\Sigma }}{\lambda }_i[y_i(w^T·z_i+b)-1]\}\\ s·\mathrm{t\mu }\ge \lambda \ge 0\end{array}\right)---\left(11\right)$

其中h为惩罚参数，h∈[0，1]，越大表示对错误分类的惩罚越大，通过$\frac{\partial L}{\partial \xi }=0$求得；ξ为分类误差；λ_i为拉格朗日乘子，μ为λ_i的最大值。由此得到

$\frac{\partial L}{\partial w}=0\to w-\underset{i=1}{\overset{l}{\Sigma }}{\lambda }_i{y}_i{z}_i=0$

$\frac{\partial L}{\partial b}=0\to \underset{i=1}{\overset{l}{\Sigma }}{\lambda }_i{y}_i=0$

可得

$w=\underset{i=1}{\overset{l}{\Sigma }}{\lambda }_i{y}_i{z}_i---\left(12\right)$

至此可求得w，b，λ，式中λ_i可能是：①λ_i＝0；②0＜λ_i＜μ；③λ_i＝μ。只有后两者所对应的λ_i对w有贡献，对应的样本成为支持向量，计为x_NSV，也就是对最优分界面、判别函数有贡献，所对应的学习方法称之为支持向量机。在支持向量中满足③所对应的z_i称为边界支持向量，实际上是错分的训练样本点；满足②所对应的z_i称为标准支持向量。

将式(12)带到式(11)，进一步处理得

$\left(\begin{array}{c}\underset{\lambda }{\max }[\omega \left(\lambda \right)=\underset{i=1}{\overset{l}{\Sigma }}{\lambda }_i-\underset{i=1}{\overset{l}{\Sigma }}\underset{j=1}{\overset{l}{\Sigma }}{\lambda }_i{\lambda }_j{y}_i{y}_j(z_i·z_j)]\\ s·t\underset{i=1}{\overset{l}{\Sigma }}{\lambda }_i{y}_i=0,\mu \ge {\lambda }_i\ge 0\end{array}\right)---\left(13\right)$

由式(13)可以求得各λ_i，最优分界面方程为：

$F=w^T·z+b=\underset{i=1}{\overset{l}{\Sigma }}{\lambda }_i{y}_i{z_i}^T·z+b=0---\left(14\right)$

由此构成的最佳分类判别函数为：

$f=\mathrm{sgn}[\underset{i=1}{\overset{l}{\Sigma }}{\lambda }_i{y}_i(z_i^T·z)+b]\mathrm{sgn}[\underset{i=1}{\overset{N_{\mathrm{NSV}}}{\Sigma }}{\lambda }_i{y}_i(z_i^T·z)+b]---\left(15\right)$

涡结构是一种线性不可分的情况，无法知道映射函数z＝φ(x)的显式形式，可以利用内积核函数的非线性映射将输入向量映射到高维空间，支持向量机能在此高维空间中给出最佳分类超平面。构造内积核函数为

$K(x_i·x)=z_i^T·z---\left(16\right)$

它表示特征空间中的一个内积。最佳分类判别函数进一步表示为

$f=\mathrm{sgn}[\underset{i=1}{\overset{l_{\mathrm{NSV}}}{\Sigma }}{\lambda }_i{y}_{i}K(x_i,x)+b]$

尽管通过非线性函数将样本数据映射到具有高维的特征空间，并在特征空间中构造最优分类超平面，但在求解最优化问题和计算判别函数时并不需要显式计算该非线性函数，而只需计算核函数，从而避免特征空间维数灾难问题。核函数的选择必须满足Merce条件。

常用的核函数有：

线性核函数    K(x_i，x)＝x_i·x       (18)

多项式核函数K(x_i，x)＝(x_i·x+1)^d    (19)

径向基核函数$K(x_i,x)=\exp (-\frac{{\left|\right|x_i-x\left|\right|}^2}{{\sigma }^2})---\left(20\right)$

径向基核函数式(20)具有较好的分类性能，比多项式核函数和线性核函数性能要好，因此本发明选用径向机核函数式(20)构造分类器。

4利用步骤3建立的最佳分类判别函数，重复步骤1、2对涡结构样本进行涡结构识别。

输入为未知涡结构样本的特征量x，经过最佳分类判别函数f进行判别，输出为1的为大涡模式，输出为-1的为小涡模式。

5利用小波反变换进行大、小尺度涡重建。

在重构中需对低频带的小波系数矩阵进行特征量的提取，大尺度涡结构的重构也只是利用低频分量进行重构。高频分量重构后得到小尺度涡结构。根据Daubiches小波预先存储的有限序列{p_n，m}和{q_n，m}，这两个序列可以通过查找有关小波分析的书得到。c_k+1；n，m′d_k+1；n，m′分别为重构的大尺度涡结构和小尺度涡结构，其他具体的符号定义同步骤1。

$c_{k+1;n,m}^{\prime }=\underset{j={-L}_1+n}{\overset{n}{\Sigma }}\underset{l={-L}_2+m}{\overset{m}{\Sigma }}{p}_{n-2l,m-2j}{c}_{k;l,j}---\left(21\right)$

$d_{k+1;n,m}^{\prime }=\underset{j=L_1+n}{\overset{n}{\Sigma }}\underset{l={-L}_2+m}{\overset{m}{\Sigma }}{q}_{n-2l,m-2j}{d^p}_{k;l,j}$

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。