CN103456020B - 基于treelet特征融合的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于treelet特征融合的遥感图像变化检测方法，主要解决现有技术中的单一类型差异图检测精度低，适用范围窄的问题。其实现步骤为：(1)输入两时相遥感图像X₁和X₂；(2)计算两幅图像的差值图，对数比值图和均值比图；(3)获取由三幅不同差异图中像素的邻域向量组成的样本矩阵；(4)用treelet方法对样本矩阵进行特征提取，得到特征向量矩阵；(5)采用Kmeans方法将特征向量矩阵聚为两类，根据聚类结果得到变化检测结果。本发明具有操作简单，抗噪性能好和检测精度高的优点，可普遍应用于环境监测、农业研究、灾害评估等领域。

Description

基于treelet特征融合的遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，涉及遥感图像变化检测方法，可用于城区规划、自然灾害评估、土地利用和土地覆盖的动态监测等领域。

背景技术

遥感图像变化检测是指通过对同一地区不同时期的两幅或多幅图像进行比较分析，根据图像之间的差异来获取地物的变化信息。遥感图像变化检测技术已成功地应用于众多领域，如环境监测、土地利用和土地覆盖的动态监测、森林或植被的变化分析、灾害评估、农业调查、城镇变化研究及在军事中的人造目标监测和地面武装部署分析。

在变化检测研究中，学者们将现有的多种变化检测方法归为不同的类别，其中最常见的是基于差异图像分析的变化检测方法。该方法通常包括3个关键的步骤：(1)图像预处理；(2)差异图的构造；(3)变化信息的提取。其中：差异图构造的好坏对变化检测的结果有着重要的影响，差值法和比值法作为两种最基本的方法，分别通过对校正后的两时相遥感图像逐像素相减和相除运算来获得差异图。光谱图像广泛采用差值法构造差异图，而SAR图像中由于相干斑噪声的影响，多采用比值法构造差异图，且比值法对校准误差不敏感。目前在用比值法构造差异图时主要采用对数比或均值比的形式。Dekker提出使用对数比的方法构造差异图，“SpecklefilteringinsatelliteSARchangedetectionimagery”Int.J.RemoteSens.,vol.19,pp.1133–1146,1998.对数比法就是对比值法取对数，它除了可以将乘性噪声转化为加性噪声外，还对比值图像的变化范围进行了压缩，因此削弱了高像素值区域的变化。均值比法就是取对应像素的邻域均值再求比值，它对噪声具有更强的鲁棒性，但是对背景区域也就是未变化区域带来过多的伪变化信息。因此，使用单一类型的差异图进行变化检测会存在检测精度低以及适用范围窄等问题。

对于变化信息的提取，传统的直接对差异图的灰度信息进行聚类或阈值分割的方法，由于没有利用空间邻域信息，抗噪性能差，变化检测精度低。

发明内容

本发明目的在于提出一种基于treelet特征融合的遥感图像变化检测方法，以解决现有技术抗噪性能差，单一类型差异图检测精度低、适用范围窄的问题，更好地检测出变化区域。

为实现上述目的，本发明的实现步骤包括如下：

(1)读入在不同时刻从同一地区获取的两幅遥感图像X₁和X₂；

(2)获取两幅遥感图像的差异图：

(2a)计算两幅遥感图像对应像素灰度值的差，并进行归一化，得到差值图Xd₁；

(2b)计算两幅遥感图像对应像素灰度值的商，并对商取对数后，再进行归一化，得到对数比值图Xd₂；

(2c)计算两幅遥感图像的均值比，并归一化，得到均值比图Xd₃；

(3)将差值图Xd₁，对数比值图Xd₂和均值比图Xd₃中每个像素所在的h×h邻域小块作为该像素的邻域向量，将位于该3幅图像相同位置处像素的邻域向量合为一个样本向量，用所有样本向量组成样本矩阵；

(4)用treelet变换对样本矩阵进行特征提取，得到特征向量矩阵；

(5)从差值图Xd₁，对数比值图Xd₂，均值比图Xd₃中任取一幅，用K-means方法，根据特征向量矩阵将所取出差异图聚成不同的两类，分别计算这两个不同类别的均值，定义均值较大的那一类为变化类，均值较小的那一类为非变化类，由此得到最终的变化检测结果。

本发明与现有技术相比具有以下优点：

1.本发明由于采用treelet变换进行特征提取，因而操作简单，正确率高，抗噪性能好；

2.本发明方法由于利用了不同差异图的有效信息和空间邻域信息进行变化检测，进一步提高了抗噪性能和变化检测精度；

3.仿真实验表明，本发明对SAR图像和光谱图像都可以得到满意的变化检测结果，鲁棒性好。

附图说明

图1是本发明的流程图；

图2是模拟遥感图像数据集和该数据集的标准变化检测图；

图3是用本发明和对比方法对模拟遥感图像数据集的变化检测结果图；

图4是墨西哥地区遥感图像数据集和和该数据集的标准变化检测图；

图5是本发明和对比方法对墨西哥郊区遥感图像数据集的变化检测结果图；

图6是Ottawa地区SAR图像数据集和该数据集的标准变化检测图；

图7是用本发明和对比方法对Ottawa地区SAR图像数据集的变化检测结果图。

具体实施方式

参照图1，本发明的具体步骤如下：

步骤1：读入在不同时刻从同一地区获取的两幅遥感图像X₁和X₂。

步骤2：获取两幅遥感图像的差异图。

2a)对遥感图像X₁在坐标(i,j)处的灰度值X₁(i,j)与遥感图像X₂在坐标(i,j)处的灰度值X₂(i,j)进行差值运算，得到差矩阵D在坐标(i,j)处的值D(i,j)：

D(i,j)＝|X₁(i,j)-X₂(i,j)|，进而得到差矩阵D＝{D(i,j)}；

2b)归一化差矩阵，得到差值图Xd₁在坐标(i,j)处的值Xd₁(i,j)：

进而得到差值图：Xd₁＝{Xd₁(i,j)}；

2c)对遥感图像X₁在坐标(i,j)处的灰度值X₁(i,j)与遥感图像X₂在坐标(i,j)处的灰度值X₂(i,j)进行对数比运算，得到对数比矩阵R在坐标(i,j)处的值R(i,j)：

进而得到对数比矩阵R＝{R(i,j)}；

2d)归一化对数比矩阵，得到对数比值图Xd₂在坐标(i,j)处的值Xd₂(i,j)：

进而得到对数比值图：Xd₂＝{Xd₂(i,j)}；

2e)对遥感图像X₁在坐标(i,j)处的灰度值X₁(i,j)与遥感图像X₂在坐标(i,j)处的灰度值X₂(i,j)进行均值比运算，得到均值比矩阵U在坐标(i,j)处的值U(i,j)：

$U(i,j)=1-min(\frac{{\mathrm{\μ}}_1(i,j)}{{\mathrm{\μ}}_2(i,j)},\frac{{\mathrm{\μ}}_2(i,j)}{{\mathrm{\μ}}_1(i,j)}),$ 进而得到均值比矩阵U＝{U(i,j)}，

其中，μ₁(i,j)和μ₂(i,j)分别表示图像X₁和X₂中以坐标(i,j)为中心的3×3邻域窗口内所有像素灰度值的平均值；

2f)归一化均值比矩阵，得到均值比图Xd₃在坐标(i,j)处的值Xd₃(i,j)：

进而得到均值比图：Xd₃＝{Xd₃(i,j)}。

步骤3：将差值图Xd₁，对数比值图Xd₂和均值比图Xd₃中每个像素所在的h×h邻域小块作为该像素的邻域向量，将位于该3幅图像相同位置处像素的邻域向量合为一个样本向量。

3a)将差值图Xd₁中每个像素所在的h×h邻域小块转化为列向量，得到该像素的邻域向量λ＝1,2,...,H×W，其中，h≥2，H表示差值图的长，W表示差值图的宽；

3b)将对数比值图Xd₂中每个像素所在的h×h邻域小块转化为列向量，得到该像素的邻域向量

3c)将均值比图Xd₃中每个像素所在的h×h邻域小块转化为列向量，得到该像素的邻域向量

3d)根据三幅差异图Xd₁，Xd₂和Xd₃大小相同，且该三幅图像中每个像素都可用它的邻域向量来表示的特性，将三幅差异图中相同位置处像素的邻域向量合为一个，得到样本向量 $y^{\mathrm{\λ}}=\left[\begin{array}{c}{x}_1^{\mathrm{\λ}}\\ x_2^{\mathrm{\λ}}\\ x_3^{\mathrm{\λ}}\end{array}\right].$

步骤4：根据步骤3中的样本向量，得到大小为3h²×HW的样本矩阵P＝[y¹,y²,...,y^H×W]，该样本矩阵中有3h²个变量，且变量之间存在冗余，可进行特征提取和维度约减。

步骤5：用treelet变换对样本矩阵P进行特征提取，得到特征向量矩阵。

5a)计算均值向量得到大小为3h²×3h²的协方差矩阵：其中，y^λ为样本矩阵P中的第λ个列向量，T表示转置；

5b)根据协方差矩阵C第η行第τ列的元素值C(η,τ)，计算相似度矩阵M中第η行第τ列的元素值，得到大小为3h²×3h²的相似度矩阵M＝{M(η,τ)}，其中，C(η,η)表示协方差矩阵C中第η行第η列的元素值，C(τ,τ)表示协方差矩阵C中第τ行第τ列的元素值，该相似度矩阵M表明了3h²个变量中两两之间的相似性；

5c)定义treelet变换的聚类层数为l＝0,1,2,...,L-1，其中L＝3h²；定义在第l＝0层，初始化基矩阵B⁽⁰⁾为3h²×3h²的单位矩阵，相似度矩阵M⁽⁰⁾＝M；定义变量下标集合为δ＝{1,2,....,3h²}；

5d)当l≠0时，从l-1层的相似度矩阵M^(l-1)中找出最大值，并将最大值的行坐标和列坐标分别记为α和β，即：

$(\mathrm{\α},\mathrm{\β})=\arg \underset{\mathrm{\η},\mathrm{\τ}\∈\mathrm{\δ}}{max}{M}^{(l-1)}(\mathrm{\η},\mathrm{\τ}),$

其中，M^(l-1)(η,τ)表示第l-1层的相似度矩阵M^(l-1)在坐标(η,τ)处的值，η<τ，η和τ均是变量下标集合δ中的元素；

5e)计算雅克比旋转矩阵J，

其中，c＝cos(θ_l)，s＝sin(θ_l)；θ_l为旋转角，

${\mathrm{\&theta;}}_l={\tan }^{-1}\&lsqb;\frac{C_{\mathrm{\&alpha;}\mathrm{\&alpha;}}^{\left(l\right)}-C_{\mathrm{\&beta;}\mathrm{\&beta;}}^{\left(l\right)}\&PlusMinus;\sqrt{{(C_{\mathrm{\&alpha;}\mathrm{\&alpha;}}^{\left(l\right)}-C_{\mathrm{\&beta;}\mathrm{\&beta;}}^{\left(l\right)})}^2+4{\left(C_{\mathrm{\&alpha;}\mathrm{\&beta;}}^{\left(l\right)}\right)}^2}}{2C_{\mathrm{\&alpha;}\mathrm{\&beta;}}^{\left(l\right)}}\&rsqb;,$ 且|θ_l|≤π/4，

式中，表示协方差矩阵C^(l)在坐标(α,α)处的值，表示协方差矩阵C^(l)在坐标(β,β)处的值，表示协方差矩阵C^(l)在坐标(α,β)处的值；

5f)根据雅克比旋转矩阵J，得到第l层的下三个矩阵，即基矩阵B^(l)＝B^(l-1)J，协方差矩阵C^(l)＝J^TC^(l-1)J，相似度矩阵M^(l)＝J^TM^(l-1)J，其中，T表示转置；

5g)根据协方差矩阵C^(l)在坐标(α,α)处的值C^(l)(α,α)和在坐标(β,β)处的值C^(l)(β,β)，定义尺度函数φ_l和细节函数ψ_l：

如果C^(l)(α,α)≥C^(l)(β,β)，定义尺度函数φ_l和细节函数ψ_l分别为基矩阵B^(l)的第α列和第β列，否则定义尺度函数φ_l和细节函数ψ_l分别为基矩阵B^(l)的第β列和第α列；

将β从变量下标集合δ中去除，即δ＝δ\{β}，其中，\表示删除操作；

5h)重复步骤4d)至4g)直至达到聚类的最高层l＝L-1，得到treelet变换的正交基矩阵B＝[φ_l,ψ_l,....,ψ₂,ψ₁]，其中，φ_l为第l层treelet变换的尺度函数，ψ_k,k＝1,2,....,l为第k层treelet变换的细节函数；

5i)取正交基矩阵B的前S个列向量，构成矩阵B_S，1≤S≤3h²，将样本矩阵中的每一个样本向量y^λ沿B_S方向投影，得到S×1维的投影向量该投影向量包含不同差异图的信息和邻域信息；

5j)根据步骤5i)中的投影向量，得到大小为S×HW的特征向量矩阵Q＝[z¹,z²,...,z^HW]。

步骤6：从差值图Xd₁，对数比值图Xd₂，均值比图Xd₃中任取一幅，用K-means方法，根据特征向量矩阵将所取出差异图聚成不同的两类。

6a)将特征向量矩阵Q中的每一列看作一个样本，用K-means方法将所有样本聚为两个不同的类别A和B；

6b)从差值图Xd₁，对数比值图Xd₂，均值比图Xd₃中任取一幅，由于所取出差异图中的每个像素对应于Q中的每一列，因此该差异图被划分为与Q的聚类结果相对应的两个不同类别A'和B'；

确定差异图中的每个像素属于A'类还是B'类完全取决于Q的聚类结果，而与差异图本身的类型无关。

步骤7：分别计算A'类和B'类的均值，定义均值较大的那一类为变化类，均值较小的那一类为非变化类，由此得到最终的变化检测结果。

本发明的效果可以通过以下仿真结果进一步说明：

1.实验条件

实验环境为：windowsXP，SPI，CPUPentium(R)4，基本频率2.4GHZ，软件平台为MatlabR2010a。

第一个数据集为模拟遥感图像数据集，如图2(a)和2(b)所示，其中图2(a)是原始图像，它是一幅ATM(AirborneThematicMapper)图像的第3波段光谱图像，图2(b)是模拟变化图像，它通过模拟地球的天气变化和电磁波的辐射特性等因素影响并人工地嵌入一些变化区域得到，图像大小均为470×335像素，灰度级为256。图2(c)为该数据集的标准变化检测图，包括153214个非变化像素和4236个变化像素。

第二个数据集是墨西哥地区遥感图像数据集，它由两幅在不同时刻拍摄的Landsat7ETM+第4波段光谱图像组成，变化区域主要是由火灾破坏了大量的植被所导致的，如图4(a)和4(b)所示，其中图4(a)为2000年4月的图像，图4(b)为2002年5月的图像，图像大小均为512×512像素，灰度级为256。图4(c)为该数据集的标准变化检测图，包括236545个非变化像素和25599个变化像素。

第三个数据集是Ottawa地区SAR图像数据集，它由两幅在不同时刻拍摄的RadarsatSAR图像组成，变化区域主要是由洪水灾情引起的，如图6(a)和6(b)所示，其中图6(a)为1997年5月图像，图6(b)为1997年8月的图像，图像大小均为290×350像素，灰度级为256。图6(c)为该数据集的标准变化检测图，包括85451个非变化像素和16049个变化像素。

2.实验评价指标

实验使用的评价指标是漏检数、误检数、总错误数和正确检测率PCC。其中，漏检数为没有检测出来的实际发生了变化的像素的总和，误检数为实际没有发生变化但被当作变化检测出来的像素的总和，总错误数是漏检数和误检数之和，正确检测率＝1-总错误数/总像素数。

3.实验内容和实验结果

用本发明方法和现有的4种变化检测方法对3个图像数据集进行变化检测，该4种对比方法分别为：采用差值运算构造差异图然后用K-means聚类的方法记作差值法；采用对数比运算构造差异图然后用K-means聚类的方法记做对数比值法；采用均值比运算构造差异图然后用K-means聚类的方法记作均值比法；CelikT在文章“Unsupervisedchangedetectioninsatelliteimagesusingprincipalcomponentanalysisandk-meansclustering”提出的方法，记作PCA-Kmeans法。其中PCA-Kmeans方法在构造差异图时，第一个和第二个数据集用差值运算构造差异图，第三个SAR图像数据集用对数比运算构造差异图。

实验1.

用本发明方法和现有的4种变化检测方法对第一个数据集进行变化检测，实验结果如图3，其中3(a)是差值法的变化检测结果图，3(b)是对数比值法的变化检测结果图，3(c)是均值比法的变化检测结果图，3(d)是PCA-Kmeans方法的变化检测结果图，3(e)是本发明的变化检测结果图。

实验2.

用本发明方法和现有的4种变化检测方法对第二个数据集进行变化检测，实验结果如图5，其中5(a)是差值法的变化检测结果图，5(b)是对数比值法的变化检测结果图，5(c)是均值比法的变化检测结果图，5(d)是PCA-Kmeans方法的变化检测结果图，5(e)是本发明的变化检测结果图。

实验3.

用本发明方法和现有的4种变化检测方法对第三个数据集进行变化检测，实验结果如图7，其中7(a)是差值法的变化检测结果图，7(b)是对数比值法的变化检测结果图，7(c)是均值比法的变化检测结果图，7(d)是PCA-Kmeans方法的变化检测结果图，7(e)是本发明的变化检测结果图。

对本发明方法和现有的4种变化检测方法对3个数据集进行变化检测的数值统计，结果如表1所示。

表1本发明与现有四种方法的变化检测数值比较

由表1可看出，本发明的PCC均高于其他对比方法。光谱图像一般用差值运算构造差异图，从图3，和表1中第一个数据集的结果可以看出，差值法相较于对数比值法和均值比法，有较少的错误数和较好的直观效果。对于第二个墨西哥数据集，从图5及表1中可看出，同第一个数据集的结果类似，差值法的效果比对数比值法的效果好，PCA-Kmeans方法使用PCA方法从邻域小块集合中提取特征具有一定的抗噪能力，比差值法的效果好，PCC比差值法高了0.51％。本发明有效地利用了不同差异图的信息和空间邻域信息，拥有最少的总错误数和最好的直观效果。

对于第三个Ottawa地区SAR图像数据集，比值运算更适合于差异图的构造，从图7及表1可以看到，对数比值法和均值比法要比差值法拥有好的直观效果和统计结果。PCA-Kmeans方法由于采用考虑邻域信息的特征提取方法，进一步降低了对数比值法的总错误数。本发明有最少的错误数。

从图3，图5，图7及表1可看出，传统的基于单一类型差异图的构造，然后进行变化信息的提取，抗噪能力差，总错误数大，检测精度低，直观效果不好。PCA-Kmeans利用了邻域小块的数据分析方法，具有一定的抗噪能力，比单纯的K-means聚类方法的效果好，但只利用了单一类型差异图的信息，检测精度不高。本发明中为了利用空间邻域信息，使用差异图中每个像素的h×h邻域小块作为该像素的邻域向量，将三幅不同差异图中相同位置上的邻域向量合成一个样本向量，构成样本矩阵，并用treelet变换对样本矩阵进行特征提取，获得的特征向量既包含不同差异图的信息，又含有邻域信息，具有较好的抗噪能力和较高的检测精度，且操作简单。普遍适用于光谱图像和SAR图像。

Claims (2)

1.一种基于treelet特征融合的遥感图像变化检测方法，包括如下步骤：

(1)读入在不同时刻从同一地区获取的两幅遥感图像X₁和X₂；

(2)获取两幅遥感图像的差异图：

(2a)计算两幅遥感图像对应像素灰度值的差，并进行归一化，得到差值图Xd₁；

(2b)计算两幅遥感图像对应像素灰度值的商，并对商取对数后，再进行归一化，得到对数比值图Xd₂；

(2c)计算两幅遥感图像的均值比，并归一化，得到均值比图Xd₃；

(3)将差值图Xd₁，对数比值图Xd₂和均值比图Xd₃中每个像素所在的h×h邻域小块作为该像素的邻域向量，将位于该3幅图像相同位置处像素的邻域向量合为一个样本向量，用所有样本向量组成样本矩阵；

(4)用treelet变换对样本矩阵进行特征提取，得到特征向量矩阵；

(5)从差值图Xd₁，对数比值图Xd₂，均值比图Xd₃中任取一幅，用K-means方法，根据特征向量矩阵将所取出差异图聚成不同的两类，分别计算这两个不同类别的均值，定义均值较大的那一类为变化类，均值较小的那一类为非变化类，由此得到最终的变化检测结果。

2.根据权利要求1中所述的方法，其中步骤(4)所述的用treelet变换对样本矩阵进行特征提取，得到特征向量矩阵，按如下步骤进行：

4a)根据步骤(3)中得到的样本矩阵P计算均值向量得到大小为3h²×3h²的协方差矩阵：其中，h≥2，y^λ为样本矩阵P中的第λ个列向量，λ＝1,2,...,H×W，T表示转置，H表示差异图的长，W表示差异图的宽，三幅差异图Xd₁，Xd₂和Xd₃具有相同的大小；

4b)根据协方差矩阵C第η行第τ列的元素值C(η,τ)，计算相似度矩阵M中第η行第τ列的元素值，得到大小为3h²×3h²的相似度矩阵M＝{M(η,τ)}，

其中，C(η,η)表示协方差矩阵C中第η行第η列的元素值，C(τ,τ)表示协方差矩阵C中第τ行第τ列的元素值，该相似度矩阵M表明了3h²个变量中两两之间的相似性；

4c)定义treelet变换的聚类层数为l＝0,1,2,...,L-1，其中L＝3h²；定义在第l＝0层，初始化基矩阵B⁽⁰⁾为3h²×3h²的单位矩阵，相似度矩阵M⁽⁰⁾＝M；定义变量下标集合为δ＝{1,2,....,3h²}；

4d)当l≠0时，从l-1层的相似度矩阵M^(l-1)中找出最大值，并将最大值的行坐标和列坐标分别记为α和β，即：

$(\mathrm{\&alpha;},\mathrm{\&beta;})=\arg \underset{\mathrm{\&eta;},\mathrm{\&tau;}\&Element;\mathrm{\&delta;}}{max}{M}^{(l-1)}(\mathrm{\&eta;},\mathrm{\&tau;})$

其中，M^(l-1)(η,τ)表示第l-1层的相似度矩阵M^(l-1)在坐标(η,τ)处的值，η<τ，η和τ均是变量下标集合δ中的元素；

4e)计算雅克比旋转矩阵J，

其中，c＝cos(θ_l)，s＝sin(θ_l)；θ_l为旋转角，

${\mathrm{\&theta;}}_l={\tan }^{-1}\&lsqb;\frac{C_{\mathrm{\&alpha;}\mathrm{\&alpha;}}^{\left(l\right)}-C_{\mathrm{\&beta;}\mathrm{\&beta;}}^{\left(l\right)}\&PlusMinus;\sqrt{{(C_{\mathrm{\&alpha;}\mathrm{\&alpha;}}^{\left(l\right)}-C_{\mathrm{\&beta;}\mathrm{\&beta;}}^{\left(l\right)})}^2+4{\left(C_{\mathrm{\&alpha;}\mathrm{\&beta;}}^{\left(l\right)}\right)}^2}}{2C_{\mathrm{\&alpha;}\mathrm{\&beta;}}^{\left(l\right)}}\&rsqb;,$ 且|θ_l|≤π/4

式中，表示协方差矩阵C^(l)在坐标(α,α)处的值，表示协方差矩阵C^(l)在坐标(β,β)处的值，表示协方差矩阵C^(l)在坐标(α,β)处的值；

4f)根据雅克比旋转矩阵J，得到第l层的下三个矩阵，即基矩阵B^(l)＝B^(l-1)J，协方差矩阵C^(l)＝J^TC^(l-1)J，相似度矩阵M^(l)＝J^TM^(l-1)J，其中，T表示转置；

4g)根据协方差矩阵C^(l)在坐标(α,α)处的值C^(l)(α,α)和在坐标(β,β)处的值C^(l)(β,β)，定义尺度函数φ_l和细节函数ψ_l：

如果C^(l)(α,α)≥C^(l)(β,β)，定义尺度函数φ_l和细节函数ψ_l分别为基矩阵B^(l)的第α列和第β列，否则定义尺度函数φ_l和细节函数ψ_l分别为基矩阵B^(l)的第β列和第α列；

将β从变量下标集合δ中去除，即δ＝δ\{β}，其中\表示删除操作；

4h)重复步骤4d)至4g)直至达到聚类的最高层l＝L-1，得到treelet变换的正交基矩阵B＝[φ_l,ψ_l,....,ψ₂,ψ₁]，其中，φ_l为第l层treelet变换的尺度函数，ψ_k,k＝1,2,....,l为第k层treelet变换的细节函数；

4i)取正交基矩阵B的前S个列向量，构成矩阵B_S，1≤S≤3h²，将样本矩阵中的每一个样本向量y^λ沿B_S方向投影，得到S×1维的投影向量

4j)根据2i)中的投影向量，得到大小为S×HW的特征向量矩阵Q＝[z¹,z²,...,z^HW]。