一种基于分形维数的数字图像分析方法

摘要

本发明公开了一种基于分形维数的数字图像分析方法，其特征是通过图像信息矩阵U提取数字图像中所包含的全部颜色属性信息，并利用图像信息矩阵U计算出像素点之间的距离，根据图像信息矩阵U中元素进行增量判断，根据增量判断结果选用增量流形算法或非增量流形算法获得与图像信息矩阵U的低维像素点空间，利用低维像素点空间计算低维像素点空间中每任意两个像素点间的距离，通过所获得的像素点间的距离计算出数字图像的分形维数，从而实现对不同的数字图像进行分类。本发明能有效地从待处理的图像中得到分形维数，且拟合误差相对最小，对数字图像具有良好的分类能力。

说明书

技术领域

本发明涉及一种图像分析方法，具体是说是一种基于分形维数的数字图像分析方法，属 于信息系统与信息管理领域。

背景技术

图像分析是指利用数学模型和图像处理的技术来分析图像的低层特征和上层结构以获得 具有一定价值的图像信息的处理过程。20世纪60年代以来，在图像分析方面已有许多研究 成果，从针对具体问题和应用的图像分析技术逐渐向建立一般理论的方向发展。图像分析主 要分为输入、分割、识别和解释等四个过程。针对图像分析的不同过程，研究人员提出了一 系列的图像分析方法，主要方法有：统计几何特征方法、随机建模方法、傅里叶分析、小波 分析、偏微分方程法和分形维数分析方法等。图像分析技术现已广泛应用于故障诊断、目标 识别和专家系统等方面中。

分形维数概念是由哈佛大学数学系教授Mandelbrot为描述一类具有某种相似性的复杂几 何对象而提出的。在此基础上，Pentland提出可利用物体表面图像的分形模型获得物体形状 和纹理的信息。此后，各国学者提出很多基于分形维数的数字图像分析方法。

目前，分形维数是运用于分形图像处理中的其它技术的主要度量工具。计算图像的分形 维数本质上是对图像表面复杂度分布的变化的一种测量。分形维数不仅反映了灰度幅值的变 化，而且还顾及了表面在不同尺度下的变化情况。

分形维数在图像上的应用是以两点为基础的：

(l)自然界中不同种类的形态物质一般具有不同的维数；

(2)自然界中的分形与图像的灰度表示之间有着一定的对应关系。

由于分形维数反映了人们对物体表面粗糙程度的感受，同时是独立于图像一定范围内分 辨率比例的，独立于视角的，而稳定存在的物质的表示的量，因而该参数在图像分析中备受 青睐。

传统基于分形维数的数字图像分析方法存在的问题包括：

①图像分析能力不足。随着待处理图像的颜色信息的增多，其颜色属性值分布空间的复 杂化增加了对图像分析的难度，传统分形维数计算方法多针对灰度图像进行维数求解，忽视 了对真彩图像中多种颜色信息的捕获，导致传统计算方法得出的结果很难全面反映不同颜色 属性对图像分析结果的影响。

②计算时间长。传统分形维数计算方法对图像的处理是以像素点为基本单位，即对图像 中的所有像素点进行计算。当需要对每个像素点反复比较时，计算方法的复杂度将大幅增加， 导致计算时间过长。此外，部分改进的分形维数计算方法虽然可以获得符合人类视觉特性的 分形维数，但是需要对图像进行分解处理，导致计算时间也大幅增加。

拓扑流形思想是在21世纪初，由Roweis S.T、Wong W.K等人提出来的，该思想认为高 维复杂数据存在着低维流形结构，每一个数据点都是整个空间的关联结构组成部分。目前该 思想还处于不断完善和发展过程中。

发明内容

本发明是为了克服现有技术存在的不足之处，提供一种基于分形维数的数字图像分析方 法，利用拓扑流形改进分形维数对图像的分析和分类能力，并能有效地从待处理的图像中得 到分形维数，且拟合误差相对最小；改善了传统方法对图像分析能力的不足和计算时间较长 等问题，从而对数字图像具有良好的分类能力。

本发明为解决技术问题采用如下技术方案：

本发明一种基于分形维数的数字图像分析方法，所述数字图像是长宽相等的正方形图像， 其特点是按如下过程进行：

A．图像信息的提取

通过图像信息矩阵U提取数字图像中所包含的全部颜色属性信息，并利用图像信息矩阵 U中元素所表示的像素点在多维空间中所对应的坐标计算出像素点之间的距离，所述颜色属 性信息为：数字图像中像素点的位置和像素点的灰度值；或为数字图像中像素点的位置和素 点的三基色分量值；

B．图像信息的转换

根据图像信息矩阵U中元素进行增量判断，根据增量判断结果选用增量流形算法或非增 量流形算法获得与图像信息矩阵U同胚的低维像素点空间Y；

C．分形维数的计算

利用步骤B所获得的低维像素点空间Y计算出在所述低维像素点空间中每任意两个像素 点间的距离，定义所获得的像素点间的距离为低维像素点空间Y中的像素点相似度，利用所 述像素点相似度获得数字图像的分形维数，根据所获得的数字图像的分形维数对不同数字图 像进行区分，实现对不同数字图像的分类。

本发明基于分形维数的数字图像分析方法的特点也在于：

所述步骤A的实现方法为：

所述图像信息矩阵U是一个m×n的矩阵，其中m为像素点中包含颜色属性信息的个数， n为像素点的个数，由每个像素点i所包含的颜色属性信息组成所述图像信息矩阵U中一个 列向量u_i，以每个列向量u_i＝[u_i1,u_i2,…,u_im]^T作为多维空间中的坐标(u_i1,u_i2,…,u_im)，则像 素点i与其它像素点j的像素点距离d_ij为：

$d_{\mathrm{ij}}=\sqrt{\Sigma {(u_{\mathrm{iq}}-u_{\mathrm{jq}})}^2}$(q=1,2,...,m)            (1)

定义所述图像信息矩阵U中的列向量u_i为像素点i的颜色属性向量。

所述步骤B的实现方法为：

定义：所述数字图像增量前矩阵为U₀，数字图像的非增量矩阵为U₁，数字图像的增量矩 阵为U₂；

若增量判断的判断结果是所述图像信息矩阵U为非增量矩阵U₁，则按照步骤B1的非增 量流形算法获得所述非增量矩阵U₁的低维像素点空间Y₁；若增量判断的判断结果是所述图 像信息矩阵U为增量矩阵U₂，则按照步骤B2的增量流形算法获得所述增量矩阵U₂的低维 像素点空间Y₂；

B1、设ε为所述数字图像中各像素点之间的邻域范围尺度；比较每一个像素点i的颜色 属性向量u_i与所述像素点距离d_ij，当像素点距离d_ij＜邻域范围尺度ε时，则像素点i的颜色 属性向量u_j为所述颜色属性向量u_i的邻域点u_ij；否则像素点i的颜色属性向量u_j就不是像 素点颜色属性向量u_i的邻域点；由所述邻域点u_ij所组成的集合称为所述颜色属性向量u_i的 邻域集；

假设所述邻域集中的邻域点u_ij的个数为k，则像素点i的颜色属性向量u_i与所述颜色属 性向量u_i的k个邻域点的线性表示之间的损失函数Φ(W)为：

$\Phi \left(W\right)=\underset{i=1}{\overset{n}{\Sigma }}{|u_{\mathrm{ij}}-\underset{j=1}{\overset{k}{\Sigma }}(W_{\mathrm{ij}}·u_{\mathrm{ij}})|}^2---\left(2\right)$

式(2)中k为任意一个小于m的正整数；W为式(2)的权值矩阵，W_ij为式(2)的权重矩阵中 元素；

所述权重矩阵中元素W_ij的限制条件为：

$\underset{j=1}{\overset{k}{\Sigma }}{W}_{\mathrm{ij}}=1---\left(3\right)$

对式(2)进行公式变换得损失函数Φ(W)为：

$\Phi \left(W\right)=w_i^T·{(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))}^T·(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))·w_i---\left(4\right)$

式(4)中W_i为所述权值矩阵W的行向量；U’为k个所述颜色属性向量u_i所组成的矩阵；

根据拉格朗日乘子法，利用式(3)与式(4)构建函数L(W_i)：

$L\left(w_i\right)=w_i^T·{(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))}^T·(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))·w_i+\lambda (w_i^T·\overrightarrow{1}-1)---\left(5\right)$

式(5)中，λ为拉格朗日因子；

对式(5)左右两边求导后的结果恒等于零获得：

${(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))}^T·(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))·w_i=c·\overrightarrow{1}---\left(6\right)$

其中c=-2/λ且取值为1；

令局部协方差矩阵Vⁱ为：

${(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))}^T·(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))=V^i---\left(7\right)$

则所述局部协方差矩阵Vⁱ中元素v_jtⁱ为：

$V_{\mathrm{jt}}^i={(u_i-u_{\mathrm{ij}})}^T·(u_i-u_{\mathrm{it}})$(t=1,2,...,m)            (8)

将式(8)代入式(6)，获得所述权值矩阵W中元素W_ij为：

$w_{\mathrm{ij}}=\frac{\underset{t=1}{\overset{k}{\Sigma }}{\left(V^i\right)}_{\mathrm{jt}}^{-1}}{\underset{a=1}{\overset{k}{\Sigma }}\underset{b=1}{\overset{k}{\Sigma }}{\left(V^i\right)}_{\mathrm{ab}}^{-1}}---\left(9\right)$

构造关于低维像素点空间Y的损失函数Φ(Y)为：

$\Phi \left(Y\right)=\underset{i=1}{\overset{n}{\Sigma }}{|y_{\mathrm{ij}}-\underset{j=1}{\overset{k}{\Sigma }}(W_{\mathrm{ij}}·y_{\mathrm{ij}})|}^2---\left(10\right)$

式(10)中，y_ij为低维像素点空间Y中元素；

对式(10)进行公式变换得损失函数Φ(Y)为：

φ(Y)=trace(Y(I-W)(I-W)^TY^T)            (11)

式(11)中，trace表示对括号中的函数求解矩阵的迹；

令n×n的对称矩阵M为：

M=(I-W)^T·(I-W)                        (12)

将式(12)代入式(11)获得所述损失函数Φ(Y)的变换式：

$\Phi \left(Y\right)=\underset{i=1}{\overset{n}{\Sigma }}\underset{j=1}{\overset{n}{\Sigma }}(M_{\mathrm{ij}}·y_i^T·y_j)---\left(13\right)$

式(13)中的M_ij为对称矩阵M中元素；

对式(10)中的损失函数Φ(Y)的限制条件为：

$\left(\begin{array}{c}\underset{i=1}{\overset{n}{\Sigma }}{y}_i=0\\ \frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}(y_i·y_i^T)=I\end{array}\right)---\left(14\right)$

式(14)中，y_i为所述低维像素点空间中Y的行向量，I为d×d的单位矩阵；

根据拉格朗日乘子法，利用式(10)和式(14)构造函数F(Y)：

F(Y)=YMY^T+λ(YY^T-nI)                    (15)

根据拉格朗日乘子法，对式(15)左右两边求导后的结果恒等于零获得：

MY^T=λY^T                                (16)

将式(12)代入式(16)获得所述图像信息矩阵U同胚的低维像素点空间Y=[y₁,y₂,…,y_n]；

B2、利用式(9)计算获得增值前矩阵U₀的权值矩阵W₀，利用式(12)获得增值前矩阵U₀的对称矩阵M₀，利用步骤B1获得增值前矩阵U₀的低维像素点空间Y₀；利用式(12)获得所 述增量矩阵U₂的对称矩阵M₂，则矩阵直和

针对增值矩阵，新增一像素点n+1，所述像素点n+1的颜色属性向量u_n+1为所述数字图 像中第n+1像素点所对应的列向量；

当所述颜色属性向量u_n+1加入后，比较每一个像素点i的颜色属性向量u_i与所加入的颜 色属性向量u_n+1的像素点距离d_in+1，当像素点距离d_in+1＜邻域范围尺度ε，则像素点n+1的 颜色属性向量u_n+1为像素点i的颜色属性向量u_i的邻域点，以颜色属性向量u_in+1替换原来 第k个邻域点u_ik，并重新计算权重矩阵中元素W_in+1；否则像素点n+1的颜色属性向量u_n+1不为像素点i的颜色属性向量u_i邻域点，保持原像素点i的颜色属性向量u_i的邻域集，保持 原权重矩阵中元素W_ij；

比较所述颜色属性向量u_n+1与所述每一个像素点i的颜色属性向量u_i的像素点距离d_n+1i， 当像素点距离d_n+1i＜邻域范围尺度ε，则像素点i的颜色属性向量u_i为所述颜色属性向量u_n+1的邻域点，并获得所述颜色属性向量u_n+1的邻域集，利用式(9)获得所述颜色属性向量u_n+1的权重矩阵中元素W_n+1j，并得到增值矩阵U₂所对应的权值矩阵W₂，否则像素点i的颜色 属性向量u_i不为所述颜色属性向量u_n+1邻域点；

利用特征值关系式式(17)和式(18)分别获得增值前矩阵U₀的对称矩阵M₀的特征值λ₀和与 所述特征值λ₀对应的特征向量Y₀以及矩阵值和M₃的特征值λ₃与所述特征值λ₃对应的特征 向量为Y₃：

|M₀-E₀λ₀|=0和|M₃-E₃λ₃|=0            (17)

(M₀-E₀λ₀)Y₀^T＝0和(M₃-E₃-λ₃)Y₃^T=0      (18)

式(17)和式(18)中M₀和M₃是方阵且分别为n₀阶矩阵和n₃阶矩阵，E₀为n₀阶单位矩阵， E₃为n₃阶单位矩阵；

取出所述特征值λ₃中与向量[Y₀,0]^T相对应的特征值，则所述特征值λ₃中剩余特征值为特 征值λ_t；

根据特征值λ_t，利用式(19)获得增值矩阵U₂的低维像素点空间Y₂=[y₁,y₂,…,y_n,y_n+1]；

(M₂-Iλ_t)Y₂^T＝0                       (19)

所述步骤C的实现方法为：

C1、利用式(20)获得所述低维像素点相似度h_ij：

$h_{\mathrm{ij}}=\sqrt{{(y_i-y_j)}^2}---\left(20\right)$

将所述低维像素点相似度h_ij作为分量组成像素点的相似度矩阵H：

C2、将所述数字图像以正方形区域进行划分，所述正方形区域的边长为S=2^t，t为整数 且取值为1≤t≤log2(B)，B为所述数字图像的边长；

利用所述像素点相似度矩阵H比较出每个正方形区域内的像素点间的最大差异值C；

将所述每个正方形区域扩展为一列正方体的盒子，所述每列盒子覆盖所述低维像素点空 间Y中各个元素，利用所述最大差异值C获得每个正方形区域上所需要覆盖的盒子数Nr：

$\mathrm{Nr}(i,j)=\frac{C}{L}---\left(22\right)$

则每个正方形区域内所需要覆盖的总的盒子数Hr为：

$\mathrm{Hr}=\underset{i,j}{\Sigma }\mathrm{Nr}(i,j)---\left(23\right)$

则所述数字图像所对应的分形维数D为：

$D=\lim \frac{\log \left(\mathrm{Hr}\right)}{\log (1/r)}---\left(24\right)$

式(24)中，边长比例参数r=B/S，利用式(24)中的盒子数Hr和边长比例参数r运用最小二 乘法求出分形维数D。

与已有技术相比，本发明有益效果体现在：

1、本发明应用广泛，适用于对任何图像进行分形维数计算，包括传统分形维数方法无 法直接进行计算的彩色图，都具有良好的计算结果，且计算效率较高，拟合误差较小。

2、本发明对彩色图像进行分形维数计算时，避免了传统分形维数方法需要将彩色图像转 换成灰度图像的繁琐步骤，根据增量判断结果选用拓扑流形中的增量流形算法或非增量流行 算法直接针对彩色图像进行分形维数的计算，改进了分形维数的分析方法。

3、本发明具有良好的图像分类能力，特别是针对全色遥感图像进行计算时，可以将各幅 遥感图像有效地区分开来，且线性拟合相对较好。

4、本发明无论在灰度纹理图像方面，还是在彩色纹理图像方面进行分形维数计算时都具 有良好的时间性能，且针对图像中增量数据，利用增量流形算法获得原始图像的低维流形并 通过所得到的低维流形直接获得增量图像的分形维数，有效降低了计算时间。

附图说明：

图1为本发明数字图像分析方法过程示意图；

图2为本发明进行测试的灰度纹理图像；

图3为本发明进行测试的彩色纹理图像；

图4为本发明进行测试的灰度图像。

具体实施方式

本实施例中所采用的数字图像是长宽相等的正方形图像，基于该数字图像分析的分形维 数分析方法按照如下过程进行：

1．图像信息的提取

通过图像信息矩阵U提取数字图像中所包含的全部颜色属性信息，图像信息矩阵U是一 个m×n的矩阵，其中m为像素点中包含颜色属性信息的个数，n为像素点的个数；颜色属性 信息为：数字图像中像素点的位置和像素点的灰度值；或为数字图像中像素点的位置和素点 的三基色分量值；

利用图像信息矩阵U中元素所表示的像素点在多维空间中所对应的坐标计算出像素点之 间的距离；由每个像素点i所包含的颜色属性信息组成图像信息矩阵U中一个列向量u_i，以 每个列向量u_i＝[u_i1,u_i2,…,u_im]^T作为多维空间中的坐标(u_i1,u_i2,…,u_im)，则像素点i与其它 像素点j的像素点距离d_ij为：

$d_{\mathrm{ij}}=\sqrt{\Sigma {(u_{\mathrm{iq}}-u_{\mathrm{jq}})}^2}$(q=1,2,...,m)            (1)

一般公知的多维空间是由m个颜色属性信息为坐标轴，以m个颜色属性信息同时为零时 所确定的坐标系，其中像素点距离d_ij为多维空间中的坐标间的距离，像素点距离d_ij需要根 据数字图像而确定，如果是灰度图像，则颜色属性信息为数字图像中像素点的行数、像素点 的列数和像素点的灰度值，其坐标轴就是三维坐标轴，即以灰度图像中像素点的行数为第一 维坐标点、像素点的列数为第二维坐标点和像素点的灰度值为第三维坐标点所构成的三维坐 标轴来表示灰度图像中各像素点的坐标，像素点距离d_ij就是在三维坐标中各坐标所表示的像 素点之间的距离；如果是彩色图则颜色属性信息为数字图像中像素点的行数、像素点的列数 和像素点的三基色分量值，其坐标轴就是五维坐标轴，即以彩色图像中像素点的行数为第一 维坐标点、像素点的列数为第二维坐标点、像素点的红色分量值为第三维坐标点、像素点的 绿色分量值为第四维坐标点和像素点的蓝色分量值为第五维坐标来点所构成的五维坐标轴来 表示彩色图像中各像素点的坐标，并利用该坐标计算出像素点距离d_ij；定义图像信息矩阵U 中的列向量u_i为像素点i的颜色属性向量，其中，i、j=1,2,…,n。

2．图像信息的转换

图像转换是依据拓扑流形求解出与图像信息矩阵U同胚的低维像素点空间Y；根据图像 信息矩阵U中元素进行增量判断，根据增量判断结果选用不同图像信息转换方式获得与图像 信息矩阵U同胚的低维像素点空间Y；

由于传统分形维数算法是直接对数字图像中像素点所对应的灰度值进行分形维数计算， 这种方法会导致图像分析能力不足和计算时间较长等两个问题，因此本方法利用拓扑流形思 想分析得出，数字图像中每一个像素点都不是孤立的单元，而是与其他像素点密切相关，共 同组成图像纹理，而图像纹理又可以通过数字图像颜色属性信息在空间上的变化和重复来表 示；因此本方法利用图像中像素点中包含的颜色属性信息，精确地对图像进行分形维数的求 解。图像信息的转换正是利用拓扑流形对图像信息矩阵U进行计算，通过像素点之间的距离 求解出与图像信息矩阵U同胚的低维像素点空间Y。

定义：数字图像增量前矩阵为U₀，数字图像的非增量矩阵为U₁，数字图像的增量矩阵为 U₂；

增量判断的过程是逐次比较图像信息矩阵U中第i列列向量u_i与数字图像增量前矩阵U₀第i列列向量u_i中是否一致，若前n列列向量完全一致且没有出现第n+1列列向量，则认为 图像信息矩阵U为为非增量矩阵U₁；若前n列列向量都完全一致并且出现第n+1列列向量， 则图像信息矩阵U为增量矩阵U₂。

若增量判断的判断结果是图像信息矩阵U为非增量矩阵U₁，则按照步骤B1的非增量流 形算法计算非增量矩阵U₁的低维像素点空间Y₁；若增量判断的判断结果是图像信息矩阵U 为增量矩阵U₂，则按照步骤B2的增量流形算法计算增量矩阵U₂的低维像素点空间Y₂；

B1、设ε为数字图像中各像素点之间的邻域范围尺度；比较每一个像素点i的颜色属性向 量u_i与像素点距离d_ij，当像素点距离d_ij＜邻域范围尺度ε时，则像素点i的颜色属性向量u_j为颜色属性向量u_i的邻域点u_ij；否则像素点i的颜色属性向量u_j就不是像素点颜色属性向量 u_i的邻域点；由邻域点u_ij所组成的集合称为颜色属性向量u_i的邻域集；

假设邻域集中的邻域点u_ij的个数为k，则像素点i的颜色属性向量u_i与颜色属性向量u_i的k个邻域点的线性表示之间的损失函数Φ(W)为：

$\Phi \left(W\right)=\underset{i=1}{\overset{n}{\Sigma }}{|u_{\mathrm{ij}}-\underset{j=1}{\overset{k}{\Sigma }}(W_{\mathrm{ij}}·u_{\mathrm{ij}})|}^2---\left(2\right)$

式(2)中k为任意一个小于m的正整数；W为式(2)的权值矩阵，W_ij为式(2)的权重矩阵中 元素；

权重矩阵中元素W_ij的限制条件为：

$\underset{j=1}{\overset{k}{\Sigma }}{W}_{\mathrm{ij}}=1---\left(3\right)$

对式(2)进行公式变换得损失函数Φ(W)为：

$\Phi \left(W\right)=w_i^T·{(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))}^T·(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))·w_i---\left(4\right)$

式(4)中W_i为权值矩阵W的行向量；U’为k个颜色属性向量u_i所组成的矩阵；

根据拉格朗日乘子法，利用式(3)与式(4)构建函数L(W_i)：

$L\left(w_i\right)=w_i^T·{(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))}^T·(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))·w_i+\lambda (w_i^T·\overrightarrow{1}-1)---\left(5\right)$

式(5)中，λ为拉格朗日因子；

对式(5)左右两边求导后的结果恒等于零获得：

${(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))}^T·(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))·w_i=c·\overrightarrow{1}---\left(6\right)$

其中c=-2/λ且取值为1；

令局部协方差矩阵Vⁱ为：

${(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))}^T·(U^{\text{'}}-\underset{j=1}{\overset{k}{\Sigma }}\left(u_{\mathrm{ij}}\right))=V^i---\left(7\right)$

则局部协方差矩阵Vⁱ中元素V_jtⁱ为：

$V_{\mathrm{jt}}^i={(u_i-u_{\mathrm{ij}})}^T·(u_i-u_{\mathrm{it}})$(t=1,2,...,m)                       (8)

将式(8)代入式(6)，获得权值矩阵W中元素W_ij为：

$w_{\mathrm{ij}}=\frac{\underset{t=1}{\overset{k}{\Sigma }}{\left(V^i\right)}_{\mathrm{jt}}^{-1}}{\underset{a=1}{\overset{k}{\Sigma }}\underset{b=1}{\overset{k}{\Sigma }}{\left(V^i\right)}_{\mathrm{ab}}^{-1}}---\left(9\right)$

构造关于低维像素点空间Y的损失函数Φ(Y)为：

$\Phi \left(Y\right)=\underset{i=1}{\overset{n}{\Sigma }}{|y_{\mathrm{ij}}-\underset{j=1}{\overset{k}{\Sigma }}(W_{\mathrm{ij}}·y_{\mathrm{ij}})|}^2---\left(10\right)$

式(10)中，y_ij为低维像素点空间Y中元素；

对式(10)进行公式变换得损失函数Φ(Y)为：

φ(Y)=trace(Y(I-W)(I-W)^TY^T)    (11)

式(11)中，trace表示求矩阵的迹的函数；

令n×n的对称矩阵M为：

M=(I-W)^T·(I-W)                (12)

将式(12)代入式(11)获得损失函数Φ(Y)的变换式：

$\Phi \left(Y\right)=\underset{i=1}{\overset{n}{\Sigma }}\underset{j=1}{\overset{n}{\Sigma }}(M_{\mathrm{ij}}·y_i^T·y_j)---\left(13\right)$

式(13)中的M_ij为对称矩阵M中元素；

对式(10)中的损失函数Φ(Y)的限制条件为：

$\left(\begin{array}{c}\underset{i=1}{\overset{n}{\Sigma }}{y}_i=0\\ \frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}(y_i·y_i^T)=I\end{array}\right)---\left(14\right)$

式(14)中，y_i为低维像素点空间中Y的行向量，I为d×d的单位矩阵；

根据拉格朗日乘子法，利用式(10)和式(14)构造函数F(Y)：

F(Y)=YMY^T+λ(YY^T-nI)                (15)

根据拉格朗日乘子法，对式(15)左右两边求导后的结果恒等于零获得：

MY^T=λY^T                            (16)

将式(12)代入式(16)获得图像信息矩阵U同胚的低维像素点空间Y=[y₁,y₂,…,y_n]；低维 像素点空间Y为式(12)中对称矩阵M的前m个最小非零特征值所对应的特征向量；

B2、利用式(9)计算获得增值前矩阵U₀的权值矩阵W₀，利用式(12)获得增值前矩阵U₀的对称矩阵M₀，利用步骤B1获得增值前矩阵U₀的低维像素点空间Y₀；利用式(12)获得增 量矩阵U₂的对称矩阵M₂，则矩阵直和

针对增值矩阵，新增一像素点n+1，像素点n+1的颜色属性向量u_n+1为数字图像中第n+1 像素点所对应的列向量；

当颜色属性向量u_n+1加入后，比较每一个像素点i的颜色属性向量u_i与所加入的颜色属 性向量u_n+1的像素点距离d_in+1，当像素点距离d_in+1＜邻域范围尺度ε，则像素点n+1的颜色 属性向量u_n+1为像素点i的颜色属性向量u_i的邻域点，以颜色属性向量u_in+1替换原来第k 个邻域点u_ik，并重新计算权重矩阵中元素W_in+1；否则像素点n+1的颜色属性向量u_n+1不 为像素点i的颜色属性向量u_i邻域点，保持原像素点i的颜色属性向量u_i的邻域集，保持原 权重矩阵中元素W_ij；

比较颜色属性向量u_n+1与每一个像素点i的颜色属性向量u_i的像素点距离d_n+1i，当像素 点距离d_n+1i＜邻域范围尺度ε，则像素点i的颜色属性向量u_i为颜色属性向量u_n+1的邻域点， 并获得颜色属性向量u_n+1的邻域集，利用式(9)获得颜色属性向量u_n+1的权重矩阵中元素 W_n+1j，并得到增值矩阵U₂所对应的权值矩阵W₂，否则像素点i的颜色属性向量u_i不为颜 色属性向量u_n+1邻域点；

利用特征值关系式式(17)和式(18)分别获得增值前矩阵U₀的对称矩阵M₀的特征值λ₀和与 特征值λ₀对应的特征向量Y₀以及矩阵值和M₃的特征值λ₃与特征值λ₃对应的特征向量为Y₃：

|M₀-E₀λ₀|=0和|M₃-E₃λ₃|=0            (17)

(M₀-E₀λ₀)Y₀^T＝0和(M₃-E₃λ₃)Y₃^T=0     (18)

式(17)和式(18)中M₀和M₃是方阵且分别为n₀阶矩阵和n₃阶矩阵，E₀为n₀阶单位矩阵， E₃为n₃阶单位矩阵；

取出特征值λ₃中与向量[Y₀,0]^T相对应的特征值，则特征值λ₃中剩余特征值为特征值λ_t；

根据特征值λ_t，利用式(19)获得增值矩阵U₂的低维像素点空间Y₂=[y₁,y₂,…,y_n,y_n+1]， 增值矩阵U₂的低维像素点空间Y₂中各个元素分别为增值矩阵U₂的低维像素点空间Y₂中的 一维向量；

(M₂-1λ_t)Y₂^T＝0                       (19)

3．分形维数的计算

利用步骤2所获得的低维像素点空间Y计算出在低维像素点空间中每任意两个像素点间 的距离，定义所获得的像素点间的距离为低维像素点空间Y中的像素点相似度，利用像素点 相似度获得数字图像的分形维数，根据所获得的数字图像的分形维数对不同数字图像进行区 分，实现对不同数字图像的分类。

(1)、利用式(20)获得低维像素点相似度h_ij：

$h_{\mathrm{ij}}=\sqrt{{(y_i-y_j)}^2}---\left(20\right)$

式(20)低维像素点相似度即表征每个低维像素点数据y_i与其他低维像素点数据y_j的距离 h_ij；

将低维像素点相似度h_ij作为分量组成像素点的相似度矩阵H：

(2)、将数字图像以正方形区域进行划分，正方形区域的边长为S=2^t，t为整数且取值为 1≤t≤log2(B)，B为数字图像的边长；

利用像素点相似度矩阵H比较出每个正方形区域内的像素点间的最大差异值C；

将每个正方形区域扩展为一列正方体的盒子，每列盒子覆盖低维像素点空间Y中各个元 素，利用最大差异值C获得每个正方形区域上所需要覆盖的盒子数Nr：

$\mathrm{Nr}(i,j)=\frac{C}{L}---\left(22\right)$

则每个正方形区域内所需要覆盖的总的盒子数Hr为：

$\mathrm{Hr}=\underset{i,j}{\Sigma }\mathrm{Nr}(i,j)---\left(23\right)$

则数字图像所对应的分形维数D为：

$D=\lim \frac{\log \left(\mathrm{Hr}\right)}{\log (1/r)}---\left(24\right)$

式(24)中，边长比例参数r=B/S，利用式(24)中盒子数Hr和边长比例参数r运用最小二乘 法求出分形维数D。

由于数字图像的分形维数概念本身就是用来描述一类具有相似性的复杂几何对象，并且 已经证明对自然景物的表面纹理映射成的数字图像具有相同分形特征的分形表面，因此将每 一幅待测数字图像通过上述三个步骤计算出对应分形维数，根据数字图像的分形维数的不同， 对数字图像进行区分，从而实现对数字图像的分类。

下面通过三个实例作进一步的说明。

实例1：利用本方法对图2中的灰度纹理图像进行分形维数计算；

根据本方法的第1步骤，将灰度纹理图像中各像素点的颜色属性信息（即像素点位置、 灰度值或三基色分量值）视为该像素点的一组矢量属性，通过一个m×n矩阵U保存起来。

根据本方法的第2步骤，对矩阵U进行图像信息转换方式，进而求出与U同胚的低维像 素点空间Y。(由于地位像素空间Y中的数据过于庞大，仅取出五个数据显示其结构)

Y=[1.5761,1.5239,0.6454,…,-2.6334,-2.0134]；

根据本方法中第3步骤计算出相应的分形维数：

D=2.0154；

本方法的拟合误差为0.0002。

同时再利用Peleg方法、Pentland方法、Sarkar方法和Gangepain方法计算出的分形维数 和拟合误差分别为2.5481和0.2407；2.7533和0.0057；2.3721和0.0322；2.7834和0.0058。

虽然5种分形维数方法均可获取灰度纹理图像的分形维数，但本方法的拟合误差相对最 小，即线性拟合相对最好，这是由于本方法根据拓扑流形，将图像中每个像素点所对应的所 有颜色属性都作为分形维数计算中的要素，图像中信息得到了充分利用。而Peleg方法需要 对每层“厚度”进行计算，所产生的累计误差导致该方法的拟合误差相对最大；Pentland方 法需要从多个方向计算分形维数，但当方向个数有限时，则该方法无法全面地反映灰度曲面 的总体复杂程度，直接影响结果的可靠性；Sarkar方法计算出的分形维数相对偏大，这是因 为该方法常存在着过度估计盒子数的情况；Gangepain方法的拟合误差小于Peleg方法、 Pentland方法和Sarkar方法，但比本方法要大，这是由于窗口因素会使该方法损失局部性特 征，其分形维数的精确度存在不足。由此可见，本方法相对于其他方法在分析灰度纹理图像 上具有较为突出的优势。

实例2：利用本方法对图3中的彩色图像进行分形维数计算，其中图3的原图为彩色图像， 由于专利申请只能用黑白图像，故显示的图片为黑白图像；

根据本方法的第1步骤，将灰度纹理图像中各像素点的颜色属性信息（即像素点位置、 灰度值或三基色分量值）视为该像素点的一组矢量属性，通过一个m×n矩阵U保存起来。

根据本方法的第2步骤，对矩阵U进行图像信息转换方式，进而求出与U同胚的低维像 素点空间Y。(由于地位像素空间Y中的数据过于庞大，仅取出五个数据显示其结构)

Y=[1.1101,1.0069,1.1530,…,-1.6602,-0.9400]；

根据本方法中第3步骤计算出相应的分形维数：

D=2.004；

本方法的拟合误差为0.0001。

同时利用传统差分盒维法和多尺度多方向纹理特征提取方法计算出的分形维数和拟合误 差分别为2.4884和0.0040；2.1907、2.3643、2.3659、2.3646和0.0321、0.0910、0.1061、0.1106。

虽然以上3种算法均计算出彩色图像的分形维数，但本发明方法的拟合误差相对最小， 即线性拟合相对最好。其中，本方法可以直接计算出彩色图像的分形维数，并且拟合误差相 对最小，这是因为本方法根据拓扑流形，将与图像中每个像素点所对应的所有颜色属性同胚 的一维向量求解出来，并利用所求解出一维向量作为分形维数计算中的要素；传统差分盒维 法的拟合误差比本方法要大，这是由于该算法在色彩转换过程中丢失了原图的部分特征信息 并且在统计盒子数时有过度估计的情况；多尺度多方向纹理特征提取方法的拟合误差相对较 大，且计算所得分形维数个数较多，很难快速、准确地把握图像信息。由此可见，本方法相 对于其他算法在彩色图像上也具有较为突出的优势。

实例3：利用本方法对图4中的部分图像数据进行增量模拟测试；

根据本发明方法的第1步骤，将灰度纹理图像中各像素点的颜色属性信息（即像素点位 置、灰度值或三基色分量值）视为该像素点的一组矢量属性，通过一个m×n矩阵U保存起来。

根据本方法的第2步骤，对矩阵U进行图像信息转换方式，进而求出与U同胚的低维像 素点空间Y₁。

Y₁=[1.6309，-0.0771，-0.7982，0.4624，-1.2180]；

而根据LLE方法对图4中中的部分图像数据进行计算，结果Y₂为：

Y₂=[1.6309，-0.0771，-0.7982，0.4624，-1.2180]；

因为Y₁=Y₂，所以实验表明本能够将灰度图像中像素点间的关系进行有效保存，并映射 出与LLE方法相一致的低维流形结构。这两种方法的运行时间情况如表1所示，从表中的用 时数据可以看出，本方法能有效地降低计算时间，提高计算效率。

表1方法的运行时间

  数据集   本发明方法用时   LLE方法用时   第一组数据   2.6   3.3   第二组数据   2.9   3.5   第三组数据   2.4   3.1

时间单位：s