CN103116761A - 一种针对图像序列的动态纹理识别方法 - Google Patents

Abstract

一种针对图像序列的动态纹理识别方法，通过三维HMT模型参数的比较进行动态纹理识别，包括对图像序列进行Surfacelet变换得到系数矩阵，即通过塔形分解实现针对图像序列的多尺度分解，通过两个2D-DFB串联的3D-DFB将三维信号分解到不同方向，通过三维矩阵存储Surfacelet变换后得到的子带数据并由提取的特征向量生成系数矩阵；然后针对系数矩阵建立三维HMT模型，即使用高斯混合模型实现系数的分布建模，使用三维HMT模型实现系数的尺度间连续性建模，使用EM算法求解HMT模型的参数。本发明还对于所处理数据包含多种类型动态纹理的情况，提出了相应的扩展方案。本发明实现起来更加简单，适用性更强，且识别效果更好。

Description

一种针对图像序列的动态纹理识别方法

技术领域

本发明涉及针对图像序列的动态纹理识别领域，尤其涉及一种结合Surfacelet变换和三维HMT模型的动态纹理识别方法。

背景技术

随着计算机技术的不断发展，信息的表现形式也越来越丰富，在所有的信息表现形式中视频占据了很大的比例，并且这个比例还将不断扩大。动态纹理是指描述某种动态景观的具有时间相关重复特征的图像序列，是一类比较特殊的视频数据。动态纹理具备空间上的重复性和时间上的稳定性，在自然界中广泛存在，比如海浪、瀑布、飘扬的旗帜、飞翔的鸟群等等。动态纹理不同于静态纹理，静态纹理只是场景某一时刻的纹理图像，而动态纹理则可以表示出场景中的纹理随时间变化的情况。由于数据信息越来越庞杂，如何对数据进行有效的识别和分类就成为一个非常重要的课题，而动态纹理识别技术就是一种识别和区分出不同种类的动态纹理的技术。目前，动态纹理识别技术的应用已经越来越广泛，例如可以用于火灾监控视频中的烟火识别，也可以用于医学图像序列中的目标定位。

经过学者们多年的研究，已经有若干动态纹理识别方法被提出来。现有的动态纹理识别方法大致可以分为四种：基于光流的方法、基于时空域几何属性的方法、基于局部时空滤波的方法、基于模型参数估计和特征提取的方法。基于光流的方法把纹理信息和运动信息组合成特征向量，计算复杂度较低而且便于实现，但是其适用性并不是很好；基于时空域几何属性的方法先找到视频数据中的一组运动轮廓，然后跟踪这些轮廓，这种动态纹理识别方法的识别率比较高，但是找到图像序列中的运动轮廓往往是比较困难的，这一弱点制约了该方法的应用；基于局部时空滤波的方法可以更为精细地对动态纹理进行划分，比如可以将动态纹理分为平稳的、连续的、不连续的、闪烁的、发光的等不同种类，但是很多动态纹理往往同时具备上述特性中的几种特性，所以这种方法的实用性并不是很强；而基于模型参数估计和特征提取的方法是目前最为成功的动态纹理识别方法，该方法从模型的角度来观察动态纹理数据，并提供了一整套动态纹理分析和识别的算法框架。

已有的动态纹理识别方法多是在视频数据的时空域上进行的，而很多图像处理方法则是在频域上进行的。这主要因为图像数据是二维数据，可以看成是二维信号，已经有许多方法可对二维信号的频谱进行分析，例如小波变换、Contourlet（轮廓波）变换等。相比之下，三维或者更高维度的信号频谱分析方法则较少，并且不够成熟。而Surfacelet（表面波）变换是目前比较成熟的三维信号分析工具，Surfacelet可以认为是将Contourlet变换拓展到三维的情况。Surfacelet变换综合使用3D-DFB（3D Directional Filter Banks，三维方向滤波器组）和多尺度分解来实现三维信号的变换和分解。Surfacelet变换首先利用多尺度分解技术将三维信号的高频成分和低频成分分离开，接着对信号的高频成分进行3D-DFB分解，从而将三维信号的高频成分分解到不同的方向中去，这样就完成了单层Surfacelet变换。然后可以对低频成分进行进一步分解，这样就使Surfacelet变换能够迭代地进行下去，直到将信号分解至满足要求的层数。Surfacelet变换可以对三维信号的频域按照不同方向和不同层次进行划分，因此可以选择Surfacelet变换对三维视频数据信号进行频域处理，然后分析经过Surfacelet变换后得到的系数，进而实现针对动态纹理的识别方法。

Surfacelet变换的结果是一系列子带，及其对应的系数矩阵，这样对动态纹理的识别就转变成了对系数矩阵的识别。最基本的方法是将子带数据的均值和方差组合为特征向量，然后通过距离函数来评价样本间的相似性，但这种处理方式过于粗糙。考虑到Surfacelet变换后的系数分布与广义高斯分布比较接近，可以针对系数建立广义高斯模型，然后通过比较高斯函数的尺度参数和形状参数来进行动态纹理的分类，但这种方法并没有考虑到Surfacelet分解后的各个层次系数之间的关系。小波域的HMT（Hidden Markov Tree，隐马尔可夫树）模型考虑到了小波不同层次系数间的关系，因此将二维的HMT模型扩展到三维，针对Surfacelet变换后的系数进行建模，然后使用极大似然准则作为样本间的相似性判别方法，通过比较图像序列的三维HMT模型参数来区分动态纹理数据，将会得到更好的分类结果。但目前尚未有相关研究出现在本领域。

发明内容

针对上述存在的技术问题，本发明目的是提供一种结合Surfacelet变换与三维HMT模型的动态纹理识别方法，以实现更高精度的针对图像序列的动态纹理识别。

为达到上述目的，本发明采用如下的技术方案：

一种针对图像序列的动态纹理识别方法，包括以下步骤：

步骤1，建立标准动态纹理模型库TL，包括以下子步骤，

步骤1.1，设动态纹理种类数目为K，采用K段包含相应标准动态纹理的图像序列为样本，构成标准动态纹理样本库DL；设第t段图像序列为D_t，t的取值为1,2…K，标准动态纹理样本库表示为DL={D₁,D₂,…D_K}；

步骤1.2，对标准动态纹理样本库DL中的每段图像序列进行Surfacelet变换与三维HMT建模，设图像序列D_t对应的模型参数为T_t，得到标准动态纹理模型库TL={T₁,T₂,…T_K}；

步骤2，基于步骤1所得标准动态纹理模型库TL，对待识别的图像序列进行分类判断，包括以下子步骤，

步骤2.1，对待识别的图像序列进行Surfacelet变换与三维HMT建模，得到对应的模型参数T_c；

步骤2.2，将步骤2.1所得模型参数T_c与标准动态纹理模型库TL中的K个模型参数逐一比较，计算待识别的图像序列的模型参数T_c基于各样本的模型参数T_t的极大似然函数值f(T_c|T₁),f(T_c|T₂),…,f(T_c|T_K)；

步骤2.3，从步骤2.3所得极大似然函数值f(T_c|T₁),f(T_c|T₂),…,f(T_c|T_K)中寻找最大值，最大值对应的模型参数T_w所属动态纹理种类为待识别的图像序列所属的动态纹理类别。

而且，步骤1.2和步骤2.1中，对一个图像序列进行Surfacelet变换的实现过程包括以下步骤，

步骤a1，通过塔形分解实现针对图像序列的多尺度分解，包括以下子步骤，

步骤a1.1，将图像序列分别通过高通滤波器和低通滤波器，分解得到对应的信号高频成分和信号低频成分；

步骤a1.2，判断信号是否已经分解至用户预先设定的层数，是则得到多尺度分解后的三维信号；否则转至步骤a1.3；

步骤a1.3，对分解得到的信号低频成分进行2倍的上采样操作，然后通过反混叠滤波器减弱上采样操作带来的频谱混叠现象；

步骤a1.4，进行3倍的下采样操作；

步骤a1.5，步骤a1.4下采样后所得的信号作为新的图像序列，转至步骤a1.1迭代执行；

步骤a2，通过两个二维方向滤波器组串联的三维方向滤波器组，对任一尺度的三维高频信号X均执行以下子步骤，

步骤a2.1，对三维高频信号X首先使用一个二维方向滤波器组沿某一方向进行分解，得到分解后的信号，此步骤分解后的信号为楔形的分解子带；

步骤a2.2，对步骤a2.1分解后的信号使用另一个二维方向滤波器组作另外方向的分解，得到由两个楔形区域重合而成的锥形分解子带，实现将三维高频信号X分解到不同的方向，得到相应尺度各方向的高频子带；

步骤a3，通过系数矩阵存储Surfacelet变换结果，包括将步骤a1所得低频子带和各尺度各方向的高频子带的数据分别存储于一个三维矩阵中，设有N个三维矩阵，从N个三维矩阵提取特征向量，由所有特征向量组成的矩阵就是存储Surfacelet变换结果的系数矩阵。

而且，步骤1.2和步骤2.1中，针对一个图像序列进行Surfacelet变换后所得系数矩阵建立三维HMT模型实现过程包括以下步骤，

步骤b1，使用高斯混合模型实现系数的分布建模，

设x为系数矩阵中的任一系数，M为高斯混合模型中标准高斯模型的数目，则系数x的概率密度函数p(x)定义为

$p\left(x\right)=\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_{j}g(x;{\mathrm{\μ}}_j,{\mathrm{\Σ}}_j)$

$g(x;{\mathrm{\μ}}_j,{\mathrm{\Σ}}_j)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_j}\exp \{-\frac{{(x-{\mathrm{\μ}}_j)}^2}{{2\mathrm{\π}}_j^2}\}$

其中∑_j，μ_j，σ_j，α_j分别表示第j个标准高斯模型对应的系数矩阵、矩阵均值、方差以及相应的权重参数；

每个系数对应一个隐状态S，取值为0或1；

步骤b2，使用三维HMT模型实现系数的尺度间连续性建模，包括以下子步骤，

步骤b2.1，基于子带对应的系数，相应地建立一个系数八叉树；

步骤b2.2，在系数八叉树的基础之上，将每个系数的隐状态与它对应的8个子节点的隐状态连接起来，形成一个隐马尔可夫概率树；

步骤b2.3，设隐马尔可夫概率树中某系数xi的状态只由父节点p_i的状态决定，而状态转移概率表示当p_i的隐状态为m时x_i的隐状态为m'的概率，其中m与m'取值为0或1；定义如下的状态转移矩阵

$\left[\begin{array}{cc}{\mathrm{\ϵ}}_{x_i,0}^{p_i,0}& {\mathrm{\ϵ}}_{x_i,0}^{p_i,0}\\ {\mathrm{\ϵ}}_{x_i,0}^{p_i,0}& {\mathrm{\ϵ}}_{x_i,0}^{p_i,0}\end{array}\right]\→\left[\begin{array}{cc}{\mathrm{\ϵ}}_{x_i,0}^{p_i,0}& 1-{\mathrm{\ϵ}}_{x_i,0}^{p_i,0}\\ 1-{\mathrm{\ϵ}}_{x_i,1}^{p_i,1}& {\mathrm{\ϵ}}_{x_i,1}^{p_i,1}\end{array}\right]$

步骤b2.4，系数x_i表示为参数向量

${\mathrm{\Θ}}_i=\{P_S,{\mathrm{\μ}}_{i,S},{\mathrm{\σ}}_{i,S},{\mathrm{\ϵ}}_{x_i,m^{\′}}^{p_i,m}\}$

其中P_S表示根节点的状态概率分布，μ_i,S，σ_i,S分别表示系数x_i在隐状态S时高斯混合模型对应的矩阵均值与方差；在给出根节点的系数状态后，把上一层结果传至下一层，作为下一层的初始状态，得到初始化的三维HMT模型参数；

步骤b2.5，设三维HMT模型T包含N个分解子带，每个子带对应的系数都有自己的参数向量Θ，则三维HMT模型表示为

T={Θ₁,Θ₂,…Θ_N}

其中，Θ₁,Θ₂,…Θ_N表示第1,2…N个子带分别对应系数的参数向量；

步骤b3，使用期望最大化算法求解三维HMT模型的参数，包括以下子步骤，

步骤b3.1，初始化，包括基于步骤b1与步骤b2得到的一个初始化的参数向量Θ⁰，并且设置迭代次数l为0；

步骤b3.2，根据系数x_i和目前的参数向量Θ^l，计算联合后验概率分布P(S|x_i,Θ^l)，其中S为系数x_i的取值为0或1的隐状态，则针对参数向量Θ^l+1的对数似然函数的期望表示为

Q(Θ^l+1|Θ^l)=E_s[lnp(x_i,S|Θ^l+1)|x_i,Θ^l]

其中p为概率密度函数，E_s为期望函数；

步骤b3.3，计算能够使期望值最大化的三维HMT模型的新的参数向量Θ^l+1，该过程的数学表达式为

Θ^l+1=argmax{Q(Θ^l+1|Θ^l)}

过程对均值μ_i，S，方差σ_i，S，状态转移概率

进行更新；

步骤b3.4，收敛判断，如果已经达到了收敛的条件则停止，否则设置l=l+1并转到步骤b3.2迭代执行。

与现有的动态纹理识别算法不同，本发明提出了结合Surfacelet变换和三维HMT模型的新方法，将图像序列看成三维信号并用Surfacelet变换进行分解处理，然后使用三维HMT模型对变换后的系数进行建模，最后基于模型参数进行动态纹理的分类与识别。本发明具有以下优点和积极效果：

（1）本发明提出的动态纹理识别方法的实现过程更加简单，因为该方法把整个图像序列作为一个整体进行处理，且基于变换系数的模型参数进行识别，从而避免了图像序列的逐帧计算与大量系数矩阵的直接比较。

（2）本发明提出的动态纹理识别方法的适用性更强，可应用于更多种类动态纹理的图像序列，且性能稳定；相比之下，现有的很多动态纹理识别方法只适用于某一种类或某些种类的动态纹理的图像序列。

（3）本发明提出的动态纹理识别方法的识别效果更好，与其它方法相比，该方法所采用的Surfacelet变换与三维HMT模型，是对图像序列的动态纹理的更加深入与准确的描述，因此更加接近动态纹理的本质特征。

具体实施方式

本发明的技术方案可由本领域技术人员采用计算机软件手段实现自动运行流程。下面以具体实施例对本发明作进一步说明。

实施例通过三维HMT模型参数进行动态纹理识别，包括以下步骤：

步骤1，针对各种标准动态纹理数据进行3D-HMT建模以构造样本模型库

与一般的模式识别算法框架类似，动态纹理的分类与识别也需要先建立一个标准的动态纹理库，具体实现过程描述如下：

步骤1.1，若支持识别的动态纹理种类数目为K，选取K段包含相应标准动态纹理的图像序列构造标准动态纹理样本库DL，若其中第t段图像序列为D_t（t的取值为1,2…K），则标准动态纹理样本库可表示为DL={D₁,D₂,…D_K}；

步骤1.2，对标准动态纹理样本库中的每段图像序列进行Surfacelet变换与三维HMT建模，若图像序列D_t对应的模型参数为T_t，则相应的标准动态纹理模型库可表示为TL={T₁,T₂,…T_K}。

具体而言，标准动态纹理数据的选择应遵循这样的准则：作为标准样本数据的图像序列必须具有代表性，每一种类的样本数据应尽量全面覆盖该类别的动态纹理的所有形态，就是说标准样本数据对该类别的动态纹理的反映要全面。

在对标准图像序列样本数据进行Surfacelet变换时，由于HMT建模与建立标准模型库的需要，Surfacelet分解的层数与方向数要保持一致。基于实验经验，分解层数的选择以2-4层为宜，而分解方向数的选择以4-6为宜。

步骤2，针对待识别的图像序列基于极大似然准则并借助标准动态纹理模型库TL进行分类判断

标准动态纹理模型库TL预先建立后，需要时再对于当前输入的图像序列进行识别，具体实现过程描述如下：

步骤2.1，对当前输入的待识别的图像序列进行Surfacelet变换与三维HMT建模，得到其对应的模型参数T_c；Surfacelet变换与三维HMT建模具体实现方式与步骤1中对标准动态纹理样本库中的每段图像序列的处理一致；

步骤2.2，与标准动态纹理模型库TL中K个样本的模型参数逐一比较，计算待识别的图像序列的模型参数T_c基于各样本的模型参数T_t的极大似然函数值f(T_c|T₁),f(T_c|T₂),…,f(T_c|T_K)；

步骤2.3，寻找最大的极大似然函数值f(T_c|T_w)=Max{f(T_c|T₁),f(T_c|T₂),…,f(T_c|T_K)}，则最大值对应的模型参数T_w所属动态纹理种类为待识别的图像序列所属的动态纹理类别。

上面的动态纹理识别方法的描述针对所处理图像序列只对应一种动态纹理的情况，若所处理图像序列包含多种类型的动态纹理，例如先后出现水流、旗帜等不同的动态场景，实施过程可相应地做如下调整：

（1）采用现有的视频聚类技术，把待识别的图像序列按照相似度测量准则进行分割与聚集，得到原始视频的若干图像序列分段；

（2）将每个视频分段作为当前输入的待识别图像序列，分别采用本发明的动态纹理识别方法进行动态纹理的类别判断。

为便于实施参考起见，进一步提供实施例的Surfacelet变换具体实现说明，包括以下步骤：

步骤a1，通过塔形分解实现针对图像序列的多尺度分解

实施例使用三种滤波器，即高通滤波器、低通滤波器、反混叠滤波器，对图像序列进行塔形分解，具体实现过程描述如下：

（1）将图像序列分别通过高通滤波器和低通滤波器，分解得到对应的信号高频成分和信号低频成分；

（2）判断信号是否已经分解至用户预先设定的层数，

是则得到多尺度分解后的三维信号，步骤a1完成；

否则转至步骤（3）作进一步操作；

（3）对分解得到的信号低频成分进行2倍的上采样操作，然后通过反混叠滤波器来减弱上采样操作带来的频谱混叠现象；

（4）接着进行3倍的下采样操作，这样小数下采样因子1.5（上采样2下采样3）可通过仅仅保持方向滤波器组的部分响应来减少混频；

（5）把采样操作之后步骤（4）所得的信号作为新的图像序列，转至步骤（1）加以迭代执行。

多尺度分解后的三维信号中包括每轮执行步骤（1）得到的信号高频成分，即多个尺度的三维高频信号，以及最后一轮执行（1）时得到的信号低频成分，即一个三维低频信号。

步骤a2，通过两个二维方向滤波器组串联的三维方向滤波器组将信号分解到不同方向

串联两个2D-DFB得到一个3D-DFB，并用于多方向分解，具体实现方式描述如下：

3D-DFB只能处理信号的高频部分，因此可依次选取步骤a1分解所得的不同尺度的三维高频信号，并逐一处理，对任一尺度的三维高频信号X均执行以下步骤：

（1）对三维高频信号X首先使用一个2D-DFB沿某一方向进行分解，得到分解后的信号，此步骤分解后的信号为楔形的分解子带；

（2）对步骤（1）分解后的信号使用另一个2D-DFB作另外方向的分解，得到由两个楔形区域重合而成的锥形分解子带，从而将三维高频信号分解到不同的方向；

若某一尺度的三维高频信号所用第一个2D-DFB为l₁级二叉树结构，所用第二个2D-DFB为l₂级二叉树结构，则该三维高频信号被划分为

个不同方向的锥形分解子带。

处理完所有尺度的三维高频信号后，得到每个尺度信号在各方向的高频子带。

步骤a3，通过系数矩阵存储Surfacelet变换后的结果

经过步骤a1执行多尺度分解与步骤a2执行多方向分解完成Surfacelet变换后，图像序列被分解为一系列子带，其中只有一个是步骤a1所得低频子带，剩下的则是步骤a2所得各尺度各方向的高频子带。这些子带中的数据存储于三维矩阵中，存储后一个子带就对应着一个三维矩阵。

Surfacelet变换后的各三维矩阵中的数据量非常庞大，无法直接将它们合成为一个特征向量。因此需要提取能够反映三维矩阵中数据整体分布状况的主要特征，例如矩阵的均值、方差、熵等，然后由矩阵特征构成特征向量。若某图像序列分解后得到N个子带，数据分别存储于N个三维矩阵。

例如当每个三维矩阵都以矩阵均值μ作为其主要特征，则由N个矩阵均值组成的一个N维向量，就是该图像序列的动态纹理的特征向量

θ=(μ₁,μ₂,…,μ_N)              （1）

同理，若每个子带对应的三维矩阵同时以矩阵均值μ和矩阵方差σ作为其主要特征，则可得到该图像序列的一个2N维特征向量

θ=(μ₁,σ₁,μ₂,σ₂,…,μ_N,σ_N)                （2）

由所有特征向量组成的矩阵，就是存储Surfacelet变换结果的系数矩阵。

为便于实施参考起见，进一步提供实施例的针对系数矩阵建立三维HMT模型具体实现说明，包括以下子步骤：

步骤b1，使用高斯混合模型实现系数的分布建模

Surfacelet变换后所得系数矩阵中的系数具有能量紧支性，即由大量的小系数和少量的大系数组成，因此使用高斯混合模型来进行系数分布的建模。

假设x为系数矩阵中的任一系数，M为高斯混合模型中标准高斯模型的数目，则系数x的概率密度函数p(x)可定义为

$p\left(x\right)=\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_{j}g(x;{\mathrm{\μ}}_j,{\mathrm{\Σ}}_j)---\left(3\right)$

$g(x;{\mathrm{\μ}}_j,{\mathrm{\Σ}}_j)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_j}\exp \{-\frac{{(x-{\mathrm{\μ}}_j)}^2}{{2\mathrm{\π}}_j^2}\}---\left(4\right)$

其中∑_j，μ_j，σ_j，α_j分别表示第j个标准高斯模型对应的系数矩阵，矩阵均值与方差，以及相应的权重参数。

实施例中设置M=2，即使用两个基本的标准高斯模型。则基于高斯混合模型的系数分布建模，每个系数x对应一个隐状态S，该变量取值为0或1。可采用状态0对应于零均值、小方差的高斯分布（即比较平低的标准高斯分布）；状态1对应于零均值、大方差的高斯分布（即比较尖锐的标准高斯分布）。

这样，Surfacelet的系数分布问题可归结为隐状态的确定问题，当隐状态确定后每个系数的分布也随之确定。

步骤b2，使用三维HMT模型实现系数的尺度间连续性建模

基于图像序列三维信号Surfacelet分解的八叉树结构特点，以及某系数状态只由其父节点状态决定的关系特点，针对动态纹理不同层次系数间的连续性通过3D-HMT建模，具体实现过程描述如下：

（1）Surfacelet变换不同尺度间的子带具有连续性，在每一层分解后，下一层子带的尺寸大小是上一层子带尺寸的1/8，即给定尺度上相邻的8个子带与上层较低分辨率尺度上的1个子带对应，基于子带对应的系数，相应地可以逐层建立一个系数八叉树；

（2）在系数八叉树的基础之上，将每个系数的隐状态与它对应的8个子节点的隐状态连接起来，进而形成一个隐马尔可夫概率树；

（3）假设隐马尔可夫概率树中某系数x_i的状态只由其父节点p_i的状态决定，而状态转移概率

表示当p_i的隐状态为m时x_i的隐状态为m'的概率，其中m与m'取值为0或1，则对应所有的m与m'组合（即00，01，10，11），定义如下的状态转移矩阵

$\left[\begin{array}{cc}{\mathrm{\ϵ}}_{x_i,0}^{p_i,0}& {\mathrm{\ϵ}}_{x_i,0}^{p_i,0}\\ {\mathrm{\ϵ}}_{x_i,0}^{p_i,0}& {\mathrm{\ϵ}}_{x_i,0}^{p_i,0}\end{array}\right]\→\left[\begin{array}{cc}{\mathrm{\ϵ}}_{x_i,0}^{p_i,0}& 1-{\mathrm{\ϵ}}_{x_i,0}^{p_i,0}\\ 1-{\mathrm{\ϵ}}_{x_i,1}^{p_i,1}& {\mathrm{\ϵ}}_{x_i,1}^{p_i,1}\end{array}\right]---\left(5\right)$

（4）Surfacelet变换后的系数x_i可表示为参数向量

${\mathrm{\Θ}}_i=\{P_S,{\mathrm{\μ}}_{i,S},{\mathrm{\σ}}_{i,S},{\mathrm{\ϵ}}_{x_i,m^{\′}}^{p_i,m}\}---\left(6\right)$

其中P_S表示根节点的状态概率分布，μ_i,S，σ_i,S分别表示系数x_i在隐状态S时高斯混合模型对应的矩阵均值与方差。因此在给出根节点的系数状态后，可把上一层结果传至下一层，作为下一层的初始状态，从而得到初始化的HMT模型参数；

（5）若一个完整的三维HMT模型T包含N个分解子带，每个子带对应的系数都有自己的参数向量Θ，则完整的三维HMT模型可以表示为

T={Θ₁,Θ₂,…Θ_N}               （7）

其中，Θ₁,Θ₂,…Θ_N表示第1,2…N个子带分别对应系数的参数向量。

步骤b3，使用期望最大化算法求解三维HMT模型的参数

在使用3D-HMT描述Surfacelet系数矩阵之前，必须先训练这个模型，通过训练来决定能够最好表征系数的三维HMT模型参数。由于一个参数向量Θ对应着一个分解子带的数据，实施例针对模型T的每个参数向量Θ，使用EM（Expectation Maximization，期望最大）算法逐个求解，具体实现过程描述如下：

（1）初始化，首先输入基于步骤b1与步骤b2得到的一个初始化的参数向量Θ⁰（即采用公式6初始化后的结果），并且设置迭代次数l为0；

（2）E步骤，根据Surfacelet变换后的系数x_i和目前的参数向量Θ^l，计算联合后验概率分布P(S|x_i,Θ^l)，其中S为系数x_i的取值为0或1的隐状态，则针对参数向量Θ^l+1的对数似然函数的期望可表示为

Q(Θ^l+1|Θ^l)=E_s[lnp(x_i,S|Θ^l+1)|x_i,Θ^l]            （8）其中p为概率密度函数，E_s为期望函数；

（3）M步骤，基于极大似然估计单调增加的特点，计算能够使期望值最大化的HMT模型的新参数向量估计Θ^l+1，该过程的数学表达式为

Θ^l+1=argmax{Q(Θ^l+1|Θ^l)}            （9）

由于系数x_i是给定隐状态的条件高斯分布，M步骤包括对均值μ_i,S，方差σ_i,S，状态转移概率的更新；

（4）收敛判断，如果已经达到了EM算法收敛的预设条件（例如Θ^l+1与Θ^l之间的差值小于用户预先设定的阈值），则算法停止，否则设置l=l+1，转到步骤（2）迭代执行。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (3)

1.一种针对图像序列的动态纹理识别方法，其特征在于，包括以下步骤：

步骤1，建立标准动态纹理模型库TL，包括以下子步骤，

步骤1.1，设动态纹理种类数目为K，采用K段包含相应标准动态纹理的图像序列为样本，构成标准动态纹理样本库DL；设第t段图像序列为D_t，t的取值为1,2…K，标准动态纹理样本库表示为DL={D₁,D₂,…D_K}；

步骤1.2，对标准动态纹理样本库DL中的每段图像序列进行Surfacelet变换与三维HMT建模，设图像序列D_t对应的模型参数为T_t，得到标准动态纹理模型库TL={T₁,T₂,…T_K}；

步骤2，基于步骤1所得标准动态纹理模型库TL，对待识别的图像序列进行分类判断，包括以下子步骤，

步骤2.1，对待识别的图像序列进行Surfacelet变换与三维HMT建模，得到对应的模型参数T_c；步骤2.2，将步骤2.1所得模型参数T_c与标准动态纹理模型库TL中的K个模型参数逐一比较，计算待识别的图像序列的模型参数T_c基于各样本的模型参数T_t的极大似然函数值f(T_c|T₁),f(T_c|T₂),…,f(T_c|T_K)；

步骤2.3，从步骤2.3所得极大似然函数值f(T_c|T₁),f(T_c|T₂),…,f(T_c|T_K)中寻找最大值，最大值对应的模型参数T_w所属动态纹理种类为待识别的图像序列所属的动态纹理类别。

2.根据权利要求1所述针对图像序列的动态纹理识别方法，其特征在于：步骤1.2和步骤2.1中，对一个图像序列进行Surfacelet变换的实现过程包括以下步骤，

步骤a1，通过塔形分解实现针对图像序列的多尺度分解，包括以下子步骤，

步骤a1.1，将图像序列分别通过高通滤波器和低通滤波器，分解得到对应的信号高频成分和信号低频成分；

步骤a1.2，判断信号是否已经分解至用户预先设定的层数，是则得到多尺度分解后的三维信号；否则转至步骤a1.3；

步骤a1.3，对分解得到的信号低频成分进行2倍的上采样操作，然后通过反混叠滤波器减弱上采样操作带来的频谱混叠现象；

步骤a1.4，进行3倍的下采样操作；

步骤a1.5，步骤a1.4下采样后所得的信号作为新的图像序列，转至步骤a1.1迭代执行；

步骤a2，通过两个二维方向滤波器组串联的三维方向滤波器组，对任一尺度的三维高频信号X均执行以下子步骤，

步骤a2.1，对三维高频信号X首先使用一个二维方向滤波器组沿某一方向进行分解，得到分解后的信号，此步骤分解后的信号为楔形的分解子带；

步骤a2.2，对步骤a2.1分解后的信号使用另一个二维方向滤波器组作另外方向的分解，得到由两个楔形区域重合而成的锥形分解子带，实现将三维高频信号X分解到不同的方向，得到相应尺度各方向的高频子带；

步骤a3，通过系数矩阵存储Surfacelet变换结果，包括将步骤a1所得低频子带和各尺度各方向的高频子带的数据分别存储于一个三维矩阵中，设有N个三维矩阵，从N个三维矩阵提取特征向量，由所有特征向量组成的矩阵就是存储Surfacelet变换结果的系数矩阵。

3.根据权利要求2所述针对图像序列的动态纹理识别方法，其特征在于：步骤1.2和步骤2.1中，针对一个图像序列进行Surfacelet变换后所得系数矩阵建立三维HMT模型实现过程包括以下步骤，

步骤b1，使用高斯混合模型实现系数的分布建模，

设x为系数矩阵中的任一系数，M为高斯混合模型中标准高斯模型的数目，则系数x的概率密度函数p(x)定义为

$p\left(x\right)=\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_{j}g(x;{\mathrm{\μ}}_j,{\mathrm{\Σ}}_j)$

$g(x;{\mathrm{\μ}}_j,{\mathrm{\Σ}}_j)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_j}\exp \{-\frac{{(x-{\mathrm{\μ}}_j)}^2}{{2\mathrm{\π}}_j^2}\}$

其中∑_j，μ_j，σ_j，α_j分别表示第j个标准高斯模型对应的系数矩阵、矩阵均值、方差以及相应的权重参数；

每个系数对应一个隐状态S，取值为0或1；

步骤b2，使用三维HMT模型实现系数的尺度间连续性建模，包括以下子步骤，

步骤b2.1，基于子带对应的系数，相应地建立一个系数八叉树；

步骤b2.2，在系数八叉树的基础之上，将每个系数的隐状态与它对应的8个子节点的隐状态连接起来，形成一个隐马尔可夫概率树；

步骤b2.3，设隐马尔可夫概率树中某系数x_i的状态只由父节点p_i的状态决定，而状态转移概率

表示当p_i的隐状态为m时x_i的隐状态为m'的概率，其中m与m'取值为0或1；定义如下的状态转移矩阵

$\left[\begin{array}{cc}{\mathrm{\ϵ}}_{x_i,0}^{p_i,0}& {\mathrm{\ϵ}}_{x_i,0}^{p_i,0}\\ {\mathrm{\ϵ}}_{x_i,0}^{p_i,0}& {\mathrm{\ϵ}}_{x_i,0}^{p_i,0}\end{array}\right]\→\left[\begin{array}{cc}{\mathrm{\ϵ}}_{x_i,0}^{p_i,0}& 1-{\mathrm{\ϵ}}_{x_i,0}^{p_i,0}\\ 1-{\mathrm{\ϵ}}_{x_i,1}^{p_i,1}& {\mathrm{\ϵ}}_{x_i,1}^{p_i,1}\end{array}\right]$

步骤b2.4，系数x_i表示为参数向量

${\mathrm{\Θ}}_i=\{P_S,{\mathrm{\μ}}_{i,S},{\mathrm{\σ}}_{i,S},{\mathrm{\ϵ}}_{x_i,m^{\′}}^{p_i,m}\}$

其中P_S表示根节点的状态概率分布，μ_i,S，σ_i,S分别表示系数x_i在隐状态S时高斯混合模型对应的矩阵均值与方差；在给出根节点的系数状态后，把上一层结果传至下一层，作为下一层的初始状态，得到初始化的三维HMT模型参数；

步骤b2.5，设三维HMT模型T包含N个分解子带，每个子带对应的系数都有自己的参数向量Θ，则三维HMT模型表示为

T={Θ₁,Θ₂,…Θ_N}

其中，Θ₁,Θ₂,…Θ_N表示第1,2…N个子带分别对应系数的参数向量；

步骤b3，使用期望最大化算法求解三维HMT模型的参数，包括以下子步骤，

步骤b3.1，初始化，包括基于步骤b1与步骤b2得到的一个初始化的参数向量Θ⁰，并且设置迭代次数l为0；

步骤b3.2，根据系数x_i和目前的参数向量Θ^l，计算联合后验概率分布P(S|x_i,Θ^l)，其中S为系数xi的取值为0或1的隐状态，则针对参数向量Θ^l+1的对数似然函数的期望表示为

Q(Θ^l+1|Θ^l)=E_s[lnp(x_i,S|Θ^l+1)|x_i,Θ^l]

其中p为概率密度函数，E_s为期望函数；

步骤b3.3，计算能够使期望值最大化的三维HMT模型的新的参数向量Θ^l+1，该过程的数学表达式为

Θ^l+1=argmax{Q(Θ^l+1|Θ^l)}

过程对均值μ_i,S，方差σ_i,S，状态转移概率进行更新；

步骤b3.4，收敛判断，如果已经达到了收敛的条件则停止，否则设置l=l+1并转到步骤b3.2迭代执行。