CN102542296A - 一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法 - Google Patents

Abstract

一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法，属于图像处理领域，本发明为解决现有二维经验模态分解方法受端点效应影响而无法有效提取图像本征信息的问题。本发明采用多变量灰色模型对给定图像进行边界延拓，再对延拓后的图像进行传统的二维经验模态分解，得到延拓的二维本征模态函数和残差，提取与原图像对应位置的二维本征模态函数和残差，作为最终分解结果。本发明充分发挥了灰色理论所需数据量少、短期预测精度高、计算速度快、对原始数据分布无特殊要求的优势，易于实现与推广应用。

Description

一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法

技术领域

本发明涉及一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法，属于图像处理领域。

背景技术

近年来，非线性非平稳数据分析方法在信号处理领域得到了广泛关注。由诺顿·黄等人在1998年提出的经验模态分解(Empirical Mode Decomposition，EMD)，以其强大的信号处理能力而成为继Fourier变换、小波变换、Wigner-Ville变换等方法之后，又一新的研究热点。2003年，Nunes等人将EMD扩展到二维情形，提出了二维经验模态分解(Bi-dimensional Empirical Mode Decomposition，BEMD)，将给定图像分解为有限个频率由高到低的二维本征模态函数(Bi-dimensional Intrinsic Mode Function，BIMF)以及一个趋势图像的和。BIMF满足两个约束条件：1)极大值点都为正，极小值点都为负(没有骑波)，且极大值点和极小值点的个数大体相等；2)上、下包络面的局部均值为零。

作为一种新兴的图像处理技术，BEMD一经提出便得到了广泛的应用，如图像压缩、纹理分析、图像融合等。同时也存在一系列公开问题，而端点效应是BEMD面临的重大挑战之一。在BEMD迭代过程中，当以极值点为节点做二维样条插值构造上、下包络面时，无法确保边界附近存在极值点(即使存在，也可能很稀疏)，这就导致由插值方法计算得到的包络面可能无法覆盖整个图像区域，使得包络面在端点处的插值精度很差，容易发生“过冲”或“欠冲”现象。更为严重的是，随着迭代次数的增加，这种“过冲”或“欠冲”现象会“扩散”到图像内部，而对整幅图像都造成“污染”，导致分解结果严重失真。

BEMD的端点效应阻碍了其在工程领域的进一步应用，虽然很多研究者已经认识到了该问题的严重性，但目前尚无一套完整公认的解决办法，这更加促使端点效应成为BEMD亟待解决的问题。已有的端点效应抑制算法不多，主要是基于镜像延拓的方法和基于纹理合成的延拓方法。其中，镜像延拓方法将最靠近图像边界的N行(或列)像素对称地映射到图像外部，然后对延拓后的图像进行传统的BEMD，得到一些延拓的二维本征模态函数BIMFs和一个延拓残差，最后提取出与原始图像对应的BIMFs和最终残差作为分解结果；基于纹理合成的延拓算法则以一种非采样的纹理合成方法一圈圈地将原始图像向外扩展，接下来也是对扩展的图像作BEMD，得到扩展的BIMFs和残差，最后提取出对应的BIMFs和残差作为最终结果。以上两种方法都有各自的缺点，前者速度较快，但精度偏低；后者精度适中，但速度较慢。

发明内容

本发明目的是为了解决现有二维经验模态分解(Bi-dimensional Empirical ModeDecomposition，BEMD)方法受端点效应影响而无法有效提取图像本征信息的问题，提供了一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法。

本发明所述一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法，该方法包括以下步骤：

步骤一、将原始图像IM_original(x，y)进行延拓迭代，初始化

IM_process(x，y)＝IM_original(x，y)，

其中，原始图像IM_original(x，y)的大小为P×P，

IM_process(x，y)为迭代过程中生成的迭代图像；

步骤二、判断迭代图像IM_process(x，y)是否已经被延拓成(P+2Q)×(P+2Q)的图像，如果是，执行步骤八；否则，执行步骤三；

步骤三、以光栅扫描的方式查看整个原始图像IM_original(x，y)，计算出待填充小块边界与原始图像IM_original(x，y)中小块边界之间的误差，并在容错误差为ε₁的范围内，从原始图像IM_original(x，y)中搜索出l个大小为T×T的小块B＝{B₁，B₂，...，B_l}，

其中，容错误差ε₁＝0.05～0.2，T＜Q；

步骤四、从B＝{B₁，B₂，...，B_l}中随机选取一个小块B_random，random＝1，2，...，l，令B₀＝B_random，B₀作为预填充块；

步骤五、以预填充块B₀边界值为多变量灰色模型MGM(1，3)的输入值，预测出预填充块B₀内部各点的像素值，并计算预测误差

步骤六、判断预测误差是否小于误差阈值ε₂，误差阈值ε₂＝0.1～0.4，如果是，执行步骤七；否则，令B₀＝B_random，并执行步骤七；

步骤七、计算预填充块B₀与迭代图像IM_process(x，y)中待填充区域重叠位置的误差面S，以最小误差路径来缝合预填充块B₀与待填充区域边界，并将预填充块B₀内部各点拼贴到相应的待填充区域，然后返回步骤二；

步骤八、令IM_extended(x，y)＝IM_process(x，y)，

其中，IM_extended(x，y)为最终延拓图像；

步骤九、对最终延拓图像IM_extended(x，y)进行的二维经验模态分解，得到I个延拓了的二维本征模态分量

和一个残差r_extended，I个延拓了的二维本征模态分量

和一个残差r_extended的尺寸都是(P+2Q)×(P+2Q)；

步骤十、提取I个延拓了的二维本征模态分量

和一个残差r_extended图像正中间P×P的块作为原始图像IM_original(x，y)的图像特征，

即I个延拓了的二维本征模态分量

正中间的P×P的块作为原始图像IM_original(x，y)的二维本征模态分量：BIMFⁱ；残差r_extended图像正中间P×P的块作为原始图像IM_original(x，y)的残差r。

本发明的优点：

1)本发明采用多变量灰色模型MGM(1，3)方法，对给定图像IM_original(x，y)进行恰当的延拓，充分发挥了灰色理论在进行数据预测时，所需数据量少、短期预测精度高、计算速度快、对原始数据分布无特殊要求的优势，有效抑制端点效应，准确提取IM_original(x，y)的本征组分BIMF。

2)本发明仅在算法开始对原始图像IM_original(x，y)边界进行延拓，在算法结束前提取出与IM_original(x，y)对应的BIMF和残差，没有改变原有的BEMD算法大框架，便于本发明的实现，有利于本发明的推广应用。

附图说明

图1是本发明所述一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法的流程图；

图2是实施方式四的流程图；

图3是实施方式五的流程图；

图4原始高光谱图像IM_original(x，y)，大小为145×145；

图5用本发明方法中，对IM_original(x，y)分布向四周预测16个像素值的延拓图像；

图6至图10为采用传统BEMD方法得到的IM_original(x，y)分解结果图；

图11至图15为采用本发明方法得到的IM_original(x，y)分解结果图。

具体实施方式

具体实施方式一：下面结合图1说明本实施方式，

本实施方式是通过以下技术方案实现的：对于给定图像IM_original(x，y)，设其原始大小和延拓后的大小分别为P×P和(P+2Q)×(P+2Q)，它向四周延拓，四周的延拓尺寸相同，四个方向的延拓区域分别为A，B，C和D四大块。从IM_original(x，y)中选取T×T(T＜Q)的最佳匹配块B₀，用多变量灰色模型MGM(1，3)(其中，“M”代表multivariate，“G”代表gray，“M”代表model，“1”代表模型阶数为1，“3”代表模型变量个数为3)预测B₀内部像素值，将其填充到A，B，C和D四大块的对应位置，各小块之间通过最小误差路径缝合。

本实施方式所述一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法，该方法包括以下步骤：

步骤一、将原始图像IM_original(x，y)进行延拓迭代，初始化IM_process(x，y)＝IM_original(x，y)，其中，原始图像IM_original(x，y)的大小为P×P，IM_process(x，y)为迭代过程中生成的迭代图像；

步骤二、判断迭代图像IM_process(x，y)是否已经被延拓成(P+2Q)×(P+2Q)的图像，如果是，执行步骤八；否则，执行步骤三；

步骤三、以光栅扫描的方式查看整个原始图像IM_original(x，y)，计算出待填充小块边界与原始图像IM_original(x，y)中小块边界之间的误差，并在容错误差为ε₁的范围内，从原始图像IM_original(x，y)中搜索出l个大小为T×T的小块B＝{B₁，B₂，...，B_l}，

其中，ε₁＝0.05～0.2，T＜Q；

步骤四、从B＝{B₁，B₂，...，B_l}中随机选取一个小块B_random，random＝1，2，...，l，令B₀＝B_random，B₀作为预填充块；

步骤五、以预填充块B₀边界值为多变量灰色模型MGM(1，3)的输入值，预测出预填充块B₀内部各点的像素值，并计算预测误差

步骤六、判断预测误差是否小于误差阈值ε₂，误差阈值ε₂＝0.1～0.4，如果是，执行

步骤七；否则，令B₀＝B_random，并执行步骤七；

步骤七、计算预填充块B₀与迭代图像IM_process(x，y)中待填充区域重叠位置的误差面S，以最小误差路径来缝合预填充块B₀与待填充区域边界，并将预填充块B₀内部各点拼贴到相应的待填充区域，然后返回步骤二；

步骤八、令IM_extended(x，y)＝IM_process(x，y)，其中，IM_extended(x，y)为最终延拓图像；

步骤九、对最终延拓图像IM_extended(x，y)进行的二维经验模态分解，得到I个延拓了的二维本征模态分量

和一个残差r_extended，I个延拓了的二维本征模态分量

和一个残差r_extended的尺寸都是(P+2Q)×(P+2Q)；

步骤十、提取I个延拓了的二维本征模态分量

和一个残差r_extended图像正中间P×P的块作为原始图像IM_original(x，y)的图像特征，即进行MGM-BEMD的结果，MGM-BEMD是基于多变量灰色模型的二维经验模态分解的缩写，

即I个延拓了的二维本征模态分量正中间的P×P的块作为原始图像IM_original(x，y)的二维本征模态分量：BIMFⁱ，i＝1，2，...，I；残差r_extended图像正中间P×P的块作为原始图像IM_original(x，y)的残差r。

具体实施方式二：本实施方式对实施方式一作进一步说明，步骤三中容错误差ε₁＝0.1。

具体实施方式三：本实施方式对实施方式一作进一步说明，步骤六中误差阈值ε₂＝0.2。

具体实施方式四：下面结合图2说明本实施方式，本实施方式对实施方式一作进一步说明，步骤五中所述以预填充块B₀边界值为多变量灰色模型的输入值，预测出预填充块B₀内部各点的像素值，并计算预测误差

的过程为：

步骤51、设定输入的预填充块B₀边界像素大小序列为

并作为多变量灰色模型MGM(1，3)的特征数据序列，

所述特征数据序列

对应的横坐标序列为

所述特征数据序列对应的纵坐标序列为

和

作为多变量灰色模型MGM(1，3)的相关因素序列；

步骤52、根据步骤51的特征数据序列和相关因素序列建立多变量灰色模型 $\mathrm{MGM}\left(\mathrm{1,3}\right)\text{=}{x}_1^{\left(0\right)}\left(k\right)+a{z}_1^{\left(1\right)}\left(k\right)=\underset{i=2}{\overset{3}{\mathrm{\Σ}}}{b}_i{x}_i^{\left(1\right)}\left(k\right),$

式中，k＝1，2，...，n，i＝1，2，3，

为背景值序列，且按公式

$z_1^{\left(1\right)}\left(k\right)=\frac{x_1^{\left(1\right)}(k-1)+x_1^{\left(1\right)}\left(k\right)}{2}$

获取；

为

的一阶累加生成序列，且按公式

$x_i^{\left(1\right)}\left(k\right)=\underset{r=1}{\overset{k}{\mathrm{\Σ}}}{x}_i^{\left(0\right)}\left(r\right)$

获取；

式中的系数a，b₂，b₃按最小二乘法

进行估计，

其中，

$B=\left[\begin{array}{ccc}-z_1^{\left(1\right)}\left(2\right)& x_2^{\left(1\right)}\left(2\right)& x_3^{\left(1\right)}\left(2\right)\\ -z_1^{\left(1\right)}\left(3\right)& x_2^{\left(1\right)}\left(3\right)& x_3^{\left(1\right)}\left(3\right)\\ \·& \·& \·\\ \·& \·& \·\\ \·& \·& \·\\ -z_1^{\left(1\right)}\left(n\right)& x_2^{\left(1\right)}\left(n\right)& x_3^{\left(1\right)}\left(n\right)\end{array}\right]$

$Y_n={(x_1^{\left(0\right)}\left(2\right),x_1^{\left(0\right)}\left(3\right),...,x_1^{\left(0\right)}\left(n\right))}^T,n\∈N^{+}$

步骤53、根据公式

$\left\{\begin{array}{c}{\hat{x}}_1^{\left(1\right)}\left(k\right)=[x_1^{\left(1\right)}\left(0\right)-\frac{1}{a}\underset{i=2}{\overset{3}{\mathrm{\Σ}}}{b}_i{x}_i^{\left(1\right)}\left(k\right)]e^{-a(k-1)}+\frac{1}{a}\underset{i=2}{\overset{3}{\mathrm{\Σ}}}{b}_i{x}_i^{\left(1\right)}\left(k\right)\\ {\hat{x}}_1^{\left(0\right)}\left(k\right)={\hat{x}}_1^{\left(1\right)}\left(k\right)-{\hat{x}}_1^{\left(1\right)}(k-1),k=\mathrm{2,3},...,n\end{array}\right.$

获取

的预测值

步骤54、令

为预填充块B₀内部各点横坐标值，令

为预填充块B₀内部各点纵坐标值，m代表被预测的点数，

则由

$\left\{\begin{array}{c}{\hat{x}}_1^{\left(1\right)}\left(k\right)=[x_1^{\left(1\right)}\left(0\right)-\frac{1}{a}\underset{i=2}{\overset{3}{\mathrm{\Σ}}}{b}_i{x}_i^{\left(1\right)}\left(k\right)]e^{-a(k-1)}+\frac{1}{a}\underset{i=2}{\overset{3}{\mathrm{\Σ}}}{b}_i{x}_i^{\left(1\right)}\left(k\right),k=n+1,n+2,...,n+m\\ {\hat{x}}_1^{\left(0\right)}\left(k\right)={\hat{x}}_1^{\left(1\right)}\left(k\right)-{\hat{x}}_1^{\left(1\right)}(k-1),k=n+2,n+3,...,n+m\end{array}\right.$

预测出预填充块B₀内部各点的像素值；

步骤55、由

计算

的残差序列ε(k)，并经

$\stackrel{\‾}{\mathrm{\Λ}}=\frac{1}{n}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}\left|\frac{\mathrm{\ϵ}\left(k\right)}{x_1^{\left(0\right)}\left(k\right)}\right|$

得出MGM(1，3)的预测误差。

计算的是平均相对误差。

经过以上五个步骤，便由MGM(1，3)预测出B₀内部各点，同时还给出了预测误差

具体实施方式五：下面结合图3说明本实施方式，本实施方式对实施方式一作进一步说明，步骤九中对最终延拓图像IM_extended(x，y)进行的二维经验模态分解，得到I个延拓了的二维本征模态分量

和一个残差r_extended的过程为：

步骤A、初始化r₁＝IM_extended(x，y)；i＝1；j＝0；SD＝1000；h_i，j＝r₁；c_i＝r₁，

其中，最终延拓图像IM_extended(x，y)的大小为(P+2Q)×(P+2Q)，

r1为第1次二维经验模态分解后的残差，

SD为终止迭代阈值，

h_i，j为第i次二维经验模态分解中经过第j次筛选后的剩余函数，

c_i为延拓了的二维本征模态分量

步骤B、令j＝j+1；h_i，(j-1)＝r_i，并通过与相邻像素对比，找出h_i，(j-1)的极大值和极小值，

r_i为第i次二维经验模态分解后的残差，

h_i，(j-1)为第i次二维经验模态分解中经过第j-1次筛选后的剩余函数；

步骤C、对步骤B检测出的第i次二维经验模态分解中经过第j-1次筛选后的剩余函数h_i，(j-1)的极大值和极小值，构造Denaulay三角形网格，由双三次样条插值法求取h_i，(j-1)的上包络e_max和下包络e_min，并计算包络均值

步骤D、计算h_ij＝h_i(j-1)-Ave和终止迭代阈值SD：

$\mathrm{SD}=\underset{x=1}{\overset{P+2Q}{\mathrm{\Σ}}}\underset{y=1}{\overset{P+2Q}{\mathrm{\Σ}}}\left[\frac{{|h_{(i-1)j}(x,y)-h_{\mathrm{ij}}(x,y)|}^2}{h_{(i-1)j}^2(x,y)}\right];$

步骤E、判断SD＜ε₀是否成立，ε₀＝0.2～0.3，如果是，执行步骤F；否则，返回步骤B；

步骤F、提取第i个延拓了的二维本征模态分量

c_i＝h_ij；和第i+1个残差r_i+1＝r_i-c_i(t)；

步骤G、判断第i+1个残差r_i+1是否为单调函数，

判断结果为否，令j＝0；i＝i+1，并返回步骤B；判断结果为是，令r_extended＝r_i+1，I＝i+1，则获取I个延拓了的二维本征模态分量

和一个残差r_extended。

从步骤E和步骤G能看出，传统BEMD有两个循环，步骤E所示的外循环控制BIMF的个数，步骤G所示的内循环以SD＜ε₀为准则通过迭代找出每个BIMF。通过内外两个循环，最终分解出图像IM_extended(x，y)的本征信息。其中，ε₀取值一般在0.2-0.3之间。

具体实施方式六：下面结合图1至图15说明本实施方式，本实施方式结合92AV3C高光谱图像的BEMD实施例阐述本发明的具体实施方式：

高光谱图像是由高光谱成像光谱仪从可见光到近红外的几百个连续波段内获取的一种三维立体图像，不但能反映地物空间信息，而且增加了一维光谱信息，其谱分辨率比全光谱图像明显提高。本实施例选取的92AV3C高光谱图像来自于AVIRIS(AirborneVisible/Infrared Imaging Spectrometer)传感器所采集的美国印地安那州西北部某农业地区的遥感观测。该图像集合包含224个连续波，从0.40μm到2.45μm大约每隔10nm一个波段。特别地，本实施例随机选取第167个波段作为原始图像IM_original(x，y)，如图4所示，其大小为145×145。

执行步骤一、初始化IM_process(x，y)＝IM_original(x，y)。其中，IM_original(x，y)为原始高光谱第167个图像，大小为145×145，此处P＝145，IM_process(x，y)为迭代过程中生成的迭代图像；

执行步骤二、判断迭代图像IM_process(x，y)是否已经被延拓成(145+2×16)×(145+2×16)的图像，此处，Q＝16，如果是，执行步骤八；否则，执行步骤三；

执行步骤三、以光栅扫描的方式查看整个原始图像IM_original(x，y)，计算出待填充小块边界与IM_original(x，y)中小块边界之间的误差，并在容错误差为ε₁＝0.1的范围内，从IM_original(x，y)中搜索出l个大小为5×5的小块B＝{B₁，B₂，...，B_l}；

执行步骤四、从B＝{B₁，B₂，...，B_l}中随机选取一个小块B_random，令B₀＝B_random为预填充块；

执行步骤五、以B₀边界值为多变量灰色模型的输入值，预测出B₀内部各点的像素值，并计算预测误差

执行步骤六、判断预测误差

是否小于识误差阈值ε₂＝0.2，如果是，执行步骤七；否则，B₀＝B_random，并执行步骤七；

执行步骤七、计算B₀与IM_process(x，y)中待填充区域重叠位置的误差面S，以最小误差路径来缝合B₀与待填充区域边界，并将B₀内部各点拼贴到相应的待填充区域，返回步骤二；

执行步骤八、IM_extended(x，y)＝IM_process(x，y)，其中，IM_extended(x，y)即为最终的延拓结果；

图5显示了对IM_original(x，y)分别向四周延拓16个像素点的图像IM_extended(x，y)，其中，粗线框内为原始图像，粗线框外为延拓部分。可以看出，IM_extended(x，y)较为准确地反映了IM_original(x，y)的变化趋势，延拓结果较为理想。

执行步骤九、对IM_extended(x，y)进行传统的BEMD，得到4个延拓了的BIMF和1个残差，分别表示为

i＝1，2，...，4和r_extended，它们的尺寸都为(145+2×16)×(145+2×16)；

执行步骤十、分别提取

i＝1，2，...，4和r_extended图像中间145×145的块作为原始图像IM_original(x，y)进行MGM-BEMD的结果，其中，BIMF表示为BIMFⁱ，i＝1，2，...，4，最终残差表示为r。

图6至图15对比了本发明方法与传统BEMD方法计算得到的BIMF组分和残差。其中，图6至图10为传统BEMD方法的分解结果，图11至图15为本发明MGM-BEMD方法的结果。可以看出，图6至图10所示的BIMF和残差存在明显的端点效应，而图11至图15的结果得到了很大改善。

总之，图6至图15充分说明了本发明方法能有效抑制传统BEMD的端点效应，准确提取给定图像的本征信息，是对传统BEMD方法的巧妙改进。

Claims (5)

1.一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法，其特征在于，该方法包括以下步骤：

步骤一、将原始图像IM_original(x，y)进行延拓迭代，初始化

IM_process(x，y)＝IM_original(x，y)，

其中，原始图像IM_original(x，y)的大小为P×P，

IM_process(x，y)为迭代过程中生成的迭代图像；

步骤二、判断迭代图像IM_process(x，y)是否已经被延拓成(P+2Q)×(P+2Q)的图像，如果是，执行步骤八；否则，执行步骤三；

步骤三、以光栅扫描的方式查看整个原始图像IM_original(x，y)，计算出待填充小块边界与原始图像IM_original(x，y)中小块边界之间的误差，并在容错误差为ε₁的范围内，从原始图像IM_original(x，y)中搜索出l个大小为T×T的小块B＝{B₁，B₂，...，B_l}，

其中，容错误差ε₁＝0.05～0.2，T＜Q；

步骤四、从B＝{B₁，B₂，...，B_l}中随机选取一个小块B_random，random＝1，2，...，l，令B₀＝B_random，B₀作为预填充块；

步骤五、以预填充块B₀边界值为多变量灰色模型MGM(1，3)的输入值，预测出预填充块B₀内部各点的像素值，并计算预测误差

步骤六、判断预测误差

是否小于误差阈值ε₂，误差阈值ε₂＝0.1～0.4，如果是，执行

步骤七；否则，令B₀＝B_random，并执行步骤七；

步骤七、计算预填充块B₀与迭代图像IM_process(x，y)中待填充区域重叠位置的误差面S，以最小误差路径来缝合预填充块B₀与待填充区域边界，并将预填充块B₀内部各点拼贴到相应的待填充区域，然后返回步骤二；

步骤八、令IM_extended(x，y)＝IM_process(x，y)，

其中，IM_extended(x，y)为最终延拓图像；

步骤九、对最终延拓图像IM_extended(x，y)进行的二维经验模态分解，得到I个延拓了的二维本征模态分量

和一个残差r_extended，I个延拓了的二维本征模态分量

和一个残差r_extended的尺寸都是(P+2Q)×(P+2Q)；

步骤十、提取I个延拓了的二维本征模态分量

和一个残差r_extended图像正中间P×P的块作为原始图像IM_original(x，y)的图像特征，

即I个延拓了的二维本征模态分量

正中间的P×P的块作为原始图像IM_original(x，y)的二维本征模态分量：BIMFⁱ；残差r_extended图像正中间P×P的块作为原始图像IM_original(x，y)的残差r。

2.根据权利要求1所述的一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法，其特征在于，步骤三中容错误差ε₁＝0.1。

3.根据权利要求1所述的一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法，其特征在于，步骤六中误差阈值ε₂＝0.2。

4.根据权利要求1所述的一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法，其特征在于，步骤五中所述以预填充块B₀边界值为多变量灰色模型的输入值，预测出预填充块B₀内部各点的像素值，并计算预测误差

的过程为：

步骤51、设定输入的预填充块B₀边界像素大小序列为

并作为多变量灰色模型MGM(1，3)的特征数据序列，

所述特征数据序列

对应的横坐标序列为所述特征数据序列

对应的纵坐标序列为

和

作为多变量灰色模型MGM(1，3)的相关因素序列；

步骤52、根据步骤51的特征数据序列和相关因素序列建立多变量灰色模型 $\mathrm{MGM}\left(\mathrm{1,3}\right)\text{=}{x}_1^{\left(0\right)}\left(k\right)+a{z}_1^{\left(1\right)}\left(k\right)=\underset{i=2}{\overset{3}{\mathrm{\Σ}}}{b}_i{x}_i^{\left(1\right)}\left(k\right),$

式中，k＝1，2，...，n，i＝1，2，3，

为背景值序列，且按公式

$z_1^{\left(1\right)}\left(k\right)=\frac{x_1^{\left(1\right)}(k-1)+x_1^{\left(1\right)}\left(k\right)}{2}$

获取；

为

的一阶累加生成序列，且按公式

$x_i^{\left(1\right)}\left(k\right)=\underset{r=1}{\overset{k}{\mathrm{\Σ}}}{x}_i^{\left(0\right)}\left(r\right)$

获取；

式中的系数a，b₂，b₃按最小二乘法

进行估计，

其中，

$B=\left[\begin{array}{ccc}-z_1^{\left(1\right)}\left(2\right)& x_2^{\left(1\right)}\left(2\right)& x_3^{\left(1\right)}\left(2\right)\\ -z_1^{\left(1\right)}\left(3\right)& x_2^{\left(1\right)}\left(3\right)& x_3^{\left(1\right)}\left(3\right)\\ \·& \·& \·\\ \·& \·& \·\\ \·& \·& \·\\ -z_1^{\left(1\right)}\left(n\right)& x_2^{\left(1\right)}\left(n\right)& x_3^{\left(1\right)}\left(n\right)\end{array}\right]$

$Y_n={(x_1^{\left(0\right)}\left(2\right),x_1^{\left(0\right)}\left(3\right),...,x_1^{\left(0\right)}\left(n\right))}^T,n\∈N^{+}$

步骤53、根据公式

$\left\{\begin{array}{c}{\hat{x}}_1^{\left(1\right)}\left(k\right)=[x_1^{\left(1\right)}\left(0\right)-\frac{1}{a}\underset{i=2}{\overset{3}{\mathrm{\Σ}}}{b}_i{x}_i^{\left(1\right)}\left(k\right)]e^{-a(k-1)}+\frac{1}{a}\underset{i=2}{\overset{3}{\mathrm{\Σ}}}{b}_i{x}_i^{\left(1\right)}\left(k\right)\\ {\hat{x}}_1^{\left(0\right)}\left(k\right)={\hat{x}}_1^{\left(1\right)}\left(k\right)-{\hat{x}}_1^{\left(1\right)}(k-1),k=\mathrm{2,3},...,n\end{array}\right.$

获取

的预测值

步骤54、令

为预填充块B₀内部各点横坐标值，令

为预填充块B₀内部各点纵坐标值，m代表被预测的点数，

则由

$\left\{\begin{array}{c}{\hat{x}}_1^{\left(1\right)}\left(k\right)=[x_1^{\left(1\right)}\left(0\right)-\frac{1}{a}\underset{i=2}{\overset{3}{\mathrm{\Σ}}}{b}_i{x}_i^{\left(1\right)}\left(k\right)]e^{-a(k-1)}+\frac{1}{a}\underset{i=2}{\overset{3}{\mathrm{\Σ}}}{b}_i{x}_i^{\left(1\right)}\left(k\right),k=n+1,n+2,...,n+m\\ {\hat{x}}_1^{\left(0\right)}\left(k\right)={\hat{x}}_1^{\left(1\right)}\left(k\right)-{\hat{x}}_1^{\left(1\right)}(k-1),k=n+2,n+3,...,n+m\end{array}\right.$

预测出预填充块B₀内部各点的像素值；

步骤55、由

计算

的残差序列ε(k)，并经

$\stackrel{\‾}{\mathrm{\Λ}}=\frac{1}{n}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}\left|\frac{\mathrm{\ϵ}\left(k\right)}{x_1^{\left(0\right)}\left(k\right)}\right|$

得出MGM(1，3)的预测误差。

5.根据权利要求1所述的一种采用基于多变量灰色模型的二维经验模态分解提取图像特征的方法，其特征在于，步骤九中对最终延拓图像IM_extended(x，y)进行的二维经验模态分解，得到I个延拓了的二维本征模态分量

和一个残差r_extended的过程为：

步骤A、初始化r₁＝IM_extended(x，y)；i＝1；j＝0；SD＝1000；h_i，j＝r₁；c_i＝r₁，

其中，最终延拓图像IM_extended(x，y)的大小为(P+2Q)×(P+2Q)，

r1为第1次二维经验模态分解后的残差，

SD为终止迭代阈值，

h_i，j为第i次二维经验模态分解中经过第j次筛选后的剩余函数，

c_i为延拓了的二维本征模态分量

步骤B、令j＝j+1；h_i，(j-1)＝r_i，并通过与相邻像素对比，找出h_i，(j-1)的极大值和极小值，

r_i为第i次二维经验模态分解后的残差，

h_i，(j-1)为第i次二维经验模态分解中经过第j-1次筛选后的剩余函数；

步骤C、对步骤B检测出的第i次二维经验模态分解中经过第j-1次筛选后的剩余函数h_i，(j-1)的极大值和极小值，构造Denaulay三角形网格，由双三次样条插值法求取h_i，(j-1)的上包络e_max和下包络e_min，并计算包络均值

步骤D、计算h_ij＝h_i(j-1)-Ave和终止迭代阈值SD：

$\mathrm{SD}=\underset{x=1}{\overset{P+2Q}{\mathrm{\Σ}}}\underset{y=1}{\overset{P+2Q}{\mathrm{\Σ}}}\left[\frac{{|h_{(i-1)j}(x,y)-h_{\mathrm{ij}}(x,y)|}^2}{h_{(i-1)j}^2(x,y)}\right];$

步骤E、判断SD＜ε₀是否成立，ε₀＝0.2～0.3，如果是，执行步骤F；否则，返回步骤B；

步骤F、提取第i个延拓了的二维本征模态分量

c_i＝h_ij；和第i+1个残差r_i+1＝r_i-c_i(t)；

步骤G、判断第i+1个残差r_i+1是否为单调函数，

判断结果为否，令j＝0；i＝i+1，并返回步骤B；

判断结果为是，令r_extended＝r_i+1，I＝i+1，则获取I个延拓了的二维本征模态分量

和一个残差r_extended。

Images