CN105761272A - 一种成像光谱混合像元中纯物质数量确定方法 - Google Patents

Abstract

本发明公开一种成像光谱混合像元中纯物质数量确定方法，通过采用成像光谱仪获取混合像元光谱数据；对信噪比较低的原始光谱进行平滑预处理；对获取像元光谱进行相似度判别，将相似度高的原始光谱进行均值中心化处理；通过简化特征值分解将光谱矩阵投影至正交子空间；将获取的子空间特征值降序归一化为权重序列；建立权重序列二维坐标系，通过一种几何距离自动判别方式，确定混合像元纯物质个数。采用本发明提供的方法，可以在无需任何先验信息的情况下，高效准确地检测出混合光谱纯物质数目。本发明不仅适用于近距离成像光谱(例如：显微成像光谱)和遥感高光谱，也可以用于非成像光谱(例如近红外光谱分析)黑色体系组分数的确定。

Description

一种成像光谱混合像元中纯物质数量确定方法

技术领域

本发明涉及成像光谱分析技术领域，尤其涉及一种成像光谱混合像元中纯物质数量确定方法。

背景技术

成像光谱技术是20世纪90年代以来发展最快、最引人注目的光谱分析技术，利用成像光谱仪对混合物质化学成分性质或含量进行分析，具有速度快、不破坏样品、操作简单、稳定性好、效率高等特点，可广泛应用于遥感、农业和食品检测等领域。成像光谱数据立方体往往包含成千上万个像元点的光谱，由于成像光谱仪空间分辨率有限或者由于被测物为混合化合物，每个像元中通常包含的不仅仅是一种纯物质，所以形成混合像元，这些混合像元的光谱由若干种纯物质光谱混合而成，因此高效准确地确定混合像元纯物质数量有助于混合物分析，并对研究混合像元解混有着非常重要的意义。

在对混合物质进行光谱测量的过程中，由于光谱仪所获取的混合物光谱信息往往是由不同纯物质的光谱信号相混合所致，基于此，结合著名的朗伯比尔定律，一种线性混合模型被广泛用于研究混合像元光谱解析。线性混合模型的基本思想是，所研究的混合像元光谱(矢量)可被看做是由纯物质光谱与纯物质丰度值相乘并且累加的结果，近年来，基于光谱线性混合模型，很多算法被开发出来用于纯物质数量估计，可将这些算法大致分为如下三类：

第一类算法是基于主成分分析(PCA)的算法。PCA的基本原理是在原始高维数据空间找到一组合适的正交基向量，其维数低于原始数据的维数，但是可以保留原数据的主要信息。PCA可以有效地估计出纯物质的数目，并对噪声信号有很好的抵抗性，但不足之处在于，当向量空间维数较大时，算法中数据协方差矩阵和矩阵特征值分解的计算复杂度和内存存储消耗都将迅速增加。

第二类方法是基于虚拟维度(VD)的算法，通过计算样本特征值检验信号存在的显著性，将最终检测出的信号源的个数作为纯物质个数，算法可以准确的确定纯物质个数，但不足之处在于对于不同的光谱数据需要设置特定的虚警率，运行效率不高。

第三类方法是一种基于最小误差的方法，其中具有代表性的是HySime方法。这是一种无监督的算法，算法通过将原始信号向子空间投影来寻找能够最小化投影误差的子空间，就以子空间的维度作为纯物质个数。算法可以完全自动获取纯物质个数，并能有效的克服噪声干扰，但不足在于需要先估计噪声矩阵，噪声矩阵估计的结果会对纯物质数目的估计产生很大影响，且随着光谱波段数目的增多，算法的估计精度会明显下降。

以上方法都能较为准确地估计出纯物质个数，但也存在如下一些问题：1)大部分算法方法需要人为设定阈值，而阈值的选择在纯像元与噪声的特征值相似时很难确定。2)一些算法计算复杂度较高，需要消耗大量内存，不具备纯物质数量检测的实时性。3)一些算法在针对一般高光谱能够很好的识别，但对于波段数目较多时的高光谱算法精度下降。

发明内容

有鉴于此，本发明的主要目的在于提供一种纯物质个数自动确定的方法用于成像光谱分析，本方法能够解决现有纯物质个数确定方法需要人为设置阈值参数、计算复杂度高、运行效率不高、光谱波段数增多检测精度下降等技术问题。

为达到上述目的，本发明的技术方案是这样实现的：一种成像光谱混合像元中纯物质数量确定方法，如图1所示，该方法步骤如下：

步骤A、通过采用成像光谱仪获取混合像元光谱数据；

步骤B、对信噪比较低的原始光谱进行平滑预处理；

步骤C、对获取像元光谱进行相似度判别，将相似度高的预处理光谱进行均值中心化处理；

步骤D、通过简化特征值分解将光谱矩阵投影至正交子空间；

步骤E、将获取的正交子空间特征值降序归一化为权重序列；

步骤F、建立权重序列二维坐标系，通过一种几何距离自动判别方式，确定混合像元纯物质数量p。

其中，所述的步骤A中通过采用成像光谱仪获取混合像元光谱数据具体为：获取所有像元点光谱反射率或吸光度得到原始光谱矩阵Y_m×n，其中m是所有像元的数目，n是光谱的波段数目。

其中，所述的步骤B中对信噪比较低的原始光谱进行平滑预处理具体为：

为降低高频噪声对光谱的影响，将原始光谱矩阵Y_m×n中的每条光谱y_1×n行预处理，采用移动窗平滑滤波处理：

$X_i^{*}=\underset{j=-t}{\overset{t}{\Σ}}{X}_{i+j}{W}_j/\underset{i=-t}{\overset{t}{\Σ}}{W}_j---\left(1\right)$

其中X_i分别是平滑后和平滑前光谱向量y_1×n中每个波段的反射率或吸光度值，W_j是移动窗口平滑中的权重因子，为简化平滑过程，取W_j＝1，2t+1是平滑点数，一般t取1或者2。

需要指出的是，本发明的预处理方法并不仅仅局限于上述方法，其他任何去噪、去背景等无用信息消除的预处理方法都应属于本发明的保护范围。

其中，所述的步骤C中对获取像元光谱进行相似度判别具体为：

C1、求取预处理后的成像光谱矩阵X_m×n的平均光谱

$\stackrel{\‾}{X}=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}x\[i\]---\left(2\right)$

其中x[i]为第i条预处理后的光谱；

C2、将预处理后的光谱矩阵X_m×n逐一选出每条光谱与均值光谱进行光谱角匹配SAM相似度判别：

$SAM(x_{il},\stackrel{\‾}{X})={\cos }^{-1}\left(\frac{\underset{l=1}{\overset{n}{\mathrm{\Σ}}}{x}_{il}\stackrel{\‾}{X}}{\sqrt{\left(\underset{l=1}{\overset{n}{\mathrm{\Σ}}}{x}_{il}^2\right)\left(\underset{l=1}{\overset{n}{\mathrm{\Σ}}}{\stackrel{\‾}{X}}^2\right)}}\right)---\left(3\right)$

其中x_il为从X_m×n中选出的第i条光谱第l个波长处的反射率值或吸光度值；

C3、如果X_m×n矩阵的每条光谱间相似度SAM值小于0.01，对X_m×n进行均值中心化处理，得到矩阵M，均值中心化处理如下：

$M=\[x_1^T-\stackrel{\‾}{X};x_2^T-\stackrel{\‾}{X};...x_m^T-\stackrel{\‾}{X}\]---\left(4\right)$

如果SAM值不小于0.01，则跳过步骤C3，并令M＝X_m×n。

其中，所述的步骤D具体包括如下步骤：

步骤D1、对M矩阵进行奇异值分解：

M＝SΛD^T(5)

式中：S为m×r阶行正交矩阵S^-1＝S^T；D为n×r阶列正交矩阵D^-1＝D^T，而Λ为r×r阶对角阵；

步骤D2、令C_n＝M^TM，对C_n矩阵进行特征值分解：

C_n＝M^TM＝DΛS^TSΛD^T＝DΛ²D^T(6)

由于D为n×r阶列正交矩阵，所以可以将其作为一组列标准正交基，将C_n矩阵投影至D正交基上，可以得到：

D^TC_n＝(D^TD)Λ²D^T＝Λ²D^T(7)。

其中，步骤E、权重序列几何判别，包括以下步骤：

所述权重序列几何距离判别方法是对正交投影后的特征值进行重要因子和次要因子分类判别，该步骤实施过程如图2所示。

所述的步骤E中将获取的子空间特征值降序归一化为权重序列具体为：所述权重序列几何距离判别方法是对正交投影后的特征值进行重要因子和次要因子分类判别，所述权重序列几何距离判别方法包括以下步骤：

步骤E1、对获取特征值矩阵Λ进行降序归一化处理：

$\left[\begin{array}{cccc}{\mathrm{\λ}}_1^2& {\mathrm{\λ}}_2^2& ...& {\mathrm{\λ}}_r^2\end{array}\right]=sort\left(\mathrm{\Λ}\right)---\left(8\right)$

$\begin{array}{c}W=1/(10^INT({\log }_{10}\left(max\right({\mathrm{\λ}}_i^i)\left)\right))\×\left[\begin{array}{cccc}{\mathrm{\λ}}_1^i& {\mathrm{\λ}}_2^i& ...& {\mathrm{\λ}}_r^i\end{array}\right]\\ s.t.i=1,2,3;W=\[w_1,w_2...w_r\]\end{array}---\left(9\right)$

其中sort()代表降序排列运算，[λ₁λ₂…λ_r]为特征值矩阵Λ主对角线上的数值，W是权重矩阵，w₁,w₂…w_r代表权重值；

步骤E2、建立权重序列坐标[A₁,A₂…A_r]＝[(S₁,w₁),(S₂,w₂)…(S_r,w_r)]，将A₁与Ar之间的连线L₁作为边界直线，该直线公式如下：

L₁:a₁x+b₁y+c₁＝0(10)

其中(S₁,S₂…S_r)为任意一组升序排列数组，如(1,2…r)；a₁、b₁、c₁代表直线方程L₁的系数；

步骤E3、计算权重序列点[A₁,A₂…A_r]到L₁的距离定义如下：

${\mathrm{DT}}_{L_1}\left(i\right)=\frac{|a_{1}r+b_1{w}_i+c_1|}{\sqrt{{a_1}^2+{b_1}^2}},i=1,2,3...r---\left(11\right)$

步骤E4、搜寻边界距离的极大值：

$p=S\left(\underset{i}{max}\right({\mathrm{DT}}_{L1}\left(i\right)\left)\right)-1,(i=1,2...r)---\left(12\right)$

步骤E5、对于第一次判别p大于4的序列，重新将A₂与Ar之间的连线L₂作为边界直线，则

iffloor(p/2)>2

L₂:a₂x+b₂y+c₂＝0(13)

其中floor()代表向下取整运算；a₂、b₂、c₂代表直线方程L₂的系数；

对于第一次判别p小于等于4的序列，跳过E5-E7；

步骤E6、计算权重序列点[A₂,A₃…A_r]到L₂的距离

${\mathrm{DT}}_{L_2}\left(i\right)=\frac{|a_{2}r+b_2{w}_i+c_2|}{\sqrt{{a_2}^2+{b_2}^2}},i=2,3...r---\left(14\right)$

步骤E7、搜寻边界距离的极大值，找出对应的序列值：

$p=S\left(\underset{i}{max}\right({\mathrm{DT}}_{L2}\left(i\right)\left)\right)-1,(i=1,2...r)---\left(15\right).$

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

(1)本发明开创性地将光谱矩阵正交投影后的权重序列点到边界直线的几何距离作为纯物质数目判定依据，该方法解决了目前大多数算法需要根据经验人为设置阈值参数的缺点，相对于常规算法具有更高的自动化程度。

(2)本发明对相似度高的光谱矩阵会先进行均值中心化处理，在一定程度上增强了微量物质的光谱信息强度，相较于常规算法能够更加准确地判定混合物质中是否含有微量物质或相似物质。

(3)本方法利用简化特征值获得权重序列，通过一种准确高效的几何距离判别方式确定混合物光谱中纯物质个数，相较于常规方法需求出协方差和相关系数矩阵，本方法计算量更小，具有更高的运行效率。

(4)本方法对原始光谱矩阵进行统一正交投影，解决了部分算法对于光谱波段数增多纯物质个数估计精度下降的问题。

附图说明

图1为本发明方法的流程图；

图2权重序列几何分离判别实施过程示意图；

图3为实验获得样本不同浓度的均值光谱图；

图4为100％三聚氰胺以及20％三聚氰胺检测权重序列图；

图5为20％三聚氰胺权重序列点的几何距离判别及其局部放大示意图；

图6本发明提出方法与Hysime、VD算法运行时间比较结果。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下举实施例并参照附图，对本发明进一步详细说明。

为验证本发明的实际使用效果，实验针对日益严重的饲料原料掺假问题进行研究。不法商家为了伪造出饲料原料蛋白质含量高的假象，常在饲料原料如豆粕中添加一些非法含高蛋白有机物如三聚氰胺等，由于豆粕外观呈现浅黄色或黄褐色，三聚氰胺为白色晶体状，豆粕的颜色会掩盖三聚氰胺的存在，而且掺杂的三聚氰胺含量很低，因此通过表观特征识别豆粕掺杂的难度较大。本发明采用显微成像光谱技术对其纯物质种类进行检测。

实验选取饲料原料豆粕和99％的三聚氰胺，光谱检测采用美国PerkinElmer公司的FITRMicroscopeSpotlight400近红外显微成像光谱仪。实验配备三聚氰胺与纯豆粕质量比为100％、20％、1.5％的样本，掺假样本采用离心搅拌机高速混合五分钟，取出少量样本用压片机进行压片，其中每个浓度配备两个压片样本。实验首先预热光谱仪30分钟，将液氮加入冷却检测器用来提高仪器灵敏度，设置实验参数为：成像模式下采集反射光谱、空间分辨率为25×25μm、光谱分辨率为4cm^-1、扫描速度为1.0cm/s、扫描波段共761个波段，所有实验均采取相同的扫描方法，并在相同的实验参数下进行。

利用显微成像光谱对掺杂样本进行测量，实验均选择同样区域，测得三聚氰胺混合浓度为100％、20％、1.5％的光谱，每个浓度梯度样本像元点个数为1024个。

步骤A、获取成像光谱数据：

获取成像光谱吸光度数据得到原始光谱矩阵Y_1024×761，其中1024是所有像元的数目，761是光谱的波段数目。图2为实验获得样本100％、20％、1.5％三个浓度下的均值光谱图；

步骤B、成像光谱反射率数据预处理

为降低高频噪声对光谱的影响，将原始光谱矩阵Y_1024×761进行5点平滑滤波处理，处理后生成新的光谱矩阵X_1024×761。

步骤C、光谱矩阵相似度判别

C1、分别求取不同浓度的样本成像光谱吸光度数据X_1024×761的平均光谱，如图3所示。

$\stackrel{\‾}{X}=\frac{1}{1024}\underset{i=1}{\overset{1024}{\Σ}}x\[i\]$

C2、将预处理后的光谱矩阵Y_1024×761随机选出一条光谱均值光谱进行SAM相似度判别：

$SAM(x_{il},\stackrel{\‾}{X})={\cos }^{-1}\left(\frac{\underset{l=1}{\overset{761}{\mathrm{\Σ}}}{x}_{il}\stackrel{\‾}{X}}{\sqrt{\left(\underset{l=1}{\overset{n}{\mathrm{\Σ}}}{x}_{il}^2\right)\left(\underset{l=1}{\overset{n}{\mathrm{\Σ}}}{\stackrel{\‾}{X}}^2\right)}}\right)$

C3、判别发现100％、20％、1.5％三聚氰胺样本SAM值均大于0.01，跳过均值中心化处理过程。

步骤D、光谱矩阵正交投影

步骤D1、分别对不同浓度的X_1024×761矩阵进行奇异值分解。

步骤D2、对不同浓度的C_n矩阵进行特征值分解，得到不同主对角特征值Λ。

步骤E、权重序列几何判别

步骤E1、对获取特征值矩阵Λ进行sort(Λ)降序处理，并将其进行归一化。

步骤E2、建立权重序列坐标[A₁,A₂…A_r]＝[(S₁,w₁),(S₂,w₂)…(S_r,w_r)]，图4所示分别为100％、20％三聚氰胺样本检测权重序列图，计算边界直线L₁。

步骤E3、计算权重序列点[A₁,A₂…A_r]到L₁的距离图5所示为20％三聚氰胺样本到直线L₁几何距离示意图。

步骤E4、搜寻边界距离的极大值对应的序列值。

步骤E5、四个浓度混合样本第一次判别p分别为1、2、2均小于4，跳过步骤E5-E7。最终检测结果如下表所示。

表格1本发明提出方法WSGSDM与Hysime、虚拟维度VD算法运行结果

表格所示为本发明的方法与目前较为流行的Hysime算法、虚拟维度算法VD比较后的结果，表中可以看出本发明提出的方法检测结果明显好于Hysime、VD算法运行结果，由于Hysime算法需要先对原始光谱的每一个波段进行噪声估计，而实际测量中噪声由测量误差和仪器误差等因素产生，单一的噪声类型无法正确判定其噪声端元，因而Hysime所估计出的端元数目包含了多余的噪声端元和变异信息，不能准确估计出端元数目，VD算法能比较准确的估计出端元数目，但在浓度较低时容易过多估计端元数目，图6所示为三种算法的运行时间，图中可以看出本发明提出的算法运行效率明显高于其他两种算法。

提供以上实施例仅仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。

Claims (6)

1.一种成像光谱混合像元中纯物质数量确定方法，其特征在于：该方法步骤如下：

步骤A、通过采用成像光谱仪获取混合像元光谱数据；

步骤B、对信噪比较低的原始光谱进行平滑预处理；

步骤C、对获取像元光谱进行相似度判别，将相似度高的预处理光谱进行均值中心化处理；

步骤D、通过简化特征值分解将光谱矩阵投影至正交子空间；

步骤E、将获取的正交子空间特征值降序归一化为权重序列；

步骤F、建立权重序列二维坐标系，通过一种几何距离自动判别方式，确定混合像元纯物质数量p。

2.根据权利要求1所述的一种成像光谱混合像元中纯物质数量确定方法，其特征在于：所述的步骤A中通过采用成像光谱仪获取混合像元光谱数据具体为：获取所有像元点光谱反射率或吸光度得到原始光谱矩阵Y_m×n，其中m是所有像元的数目，n是光谱的波段数目。

3.根据权利要求2所述的一种成像光谱混合像元中纯物质数量确定方法，其特征在于：所述的步骤B中对信噪比较低的原始光谱进行平滑预处理具体为：

为降低高频噪声对光谱的影响，将原始光谱矩阵Y_m×n中的每条光谱y_1×n行预处理，采用移动窗平滑滤波处理：

$X_i^{*}=\underset{j=-t}{\overset{t}{\Σ}}{X}_{i+j}{W}_j\underset{i=-t}{\overset{t}{\Σ}}{W}_j---\left(1\right)$

其中X_i分别是平滑后和平滑前光谱向量y_1×n中每个波段的反射率或吸光度值，W_j是移动窗口平滑中的权重因子，为简化平滑过程，取W_j＝1，2t+1是平滑点数，一般t取1或者2。

4.根据权利要求2所述的一种成像光谱混合像元中纯物质数量确定方法，其特征在于：所述的步骤C中对获取像元光谱进行相似度判别具体为：

C1、求取预处理后的成像光谱矩阵X_m×n的平均光谱

$\stackrel{\‾}{X}=\frac{1}{m}\underset{i=1}{\overset{m}{\Σ}}x\[i\]---\left(2\right)$

其中x[i]为第i条预处理后的光谱；

C2、将预处理后的光谱矩阵X_m×n逐一选出每条光谱与均值光谱进行光谱角匹配SAM相似度判别：

$SAM(x_{il},\stackrel{\‾}{X})={\cos }^{-1}\left(\frac{\underset{l=1}{\overset{n}{\Σ}}{x}_{il}\stackrel{\‾}{X}}{\sqrt{\left(\underset{l=1}{\overset{n}{\Σ}}{x}_{il}^2\right)\left(\underset{l=1}{\overset{n}{\Σ}}{\stackrel{\‾}{X}}^2\right)}}\right)---\left(3\right)$

其中x_il为从X_m×n中选出的第i条光谱第l个波长处的反射率值或吸光度值；

C3、如果X_m×n矩阵的每条光谱间相似度SAM值小于0.01，对X_m×n进行均值中心化处理，得到矩阵M，均值中心化处理如下：

$M=\[x_1^T-\stackrel{\‾}{X};x_2^T-\stackrel{\‾}{X};...x_m^T-\stackrel{\‾}{X}\]---\left(4\right)$

如果SAM值不小于0.01，则跳过步骤C3，并令M＝X_m×n。

5.根据权利要求1所述的一种成像光谱混合像元中纯物质数量确定方法，其特征在于：所述的步骤D具体包括如下步骤：

步骤D1、对M矩阵进行奇异值分解：

M＝SΛD^T(5)

式中：S为m×r阶行正交矩阵S^-1＝S^T；D为n×r阶列正交矩阵D^-1＝D^T，而Λ为r×r阶对角阵；

步骤D2、令C_n＝M^TM，对C_n矩阵进行特征值分解：

C_n＝M^TM＝DΛS^TSΛD^T＝DΛ²D^T(6)

由于D为n×r阶列正交矩阵，所以可以将其作为一组列标准正交基，将C_n矩阵投影至D正交基上，可以得到：

D^TC_n＝(D^TD)Λ²D^T＝Λ²D^T(7)。

6.根据权利要求1所述的一种成像光谱混合像元中纯物质数量确定方法，其特征在于：所述的步骤E中将获取的子空间特征值降序归一化为权重序列具体为：所述权重序列几何距离判别方法是对正交投影后的特征值进行重要因子和次要因子分类判别，所述权重序列几何距离判别方法包括以下步骤：

步骤E1、对获取特征值矩阵Λ进行降序归一化处理：

$\left[\begin{array}{cccc}{\mathrm{\λ}}_1^2& {\mathrm{\λ}}_2^2& \mathrm{...}& {\mathrm{\λ}}_r^2\end{array}\right]=sort\left(\mathrm{\Λ}\right)---\left(8\right)$

$W=1/(10^INT({\log }_{10}\left(max\left({\mathrm{\λ}}_i^i\right)\right)\left)\right)\×\left[\begin{array}{cccc}{\mathrm{\λ}}_1^i& {\mathrm{\λ}}_2^i& \mathrm{...}& {\mathrm{\λ}}_r^i\end{array}\right]---\left(9\right)$

s.t.i＝1,2,3；W＝[w₁,w₂…w_r]

其中sort()代表降序排列运算，[λ₁λ₂…λ_r]为特征值矩阵Λ主对角线上的数值，W是权重矩阵，w₁,w₂…w_r代表权重值；

步骤E2、建立权重序列坐标[A₁,A₂…A_r]＝[(S₁,w₁),(S₂,w₂)…(S_r,w_r)]，将A₁与Ar之间的连线L₁作为边界直线，该直线公式如下：

L₁:a₁x+b₁y+c₁＝0(10)

其中(S₁,S₂…S_r)为任意一组升序排列数组，如(1,2…r)；a₁、b₁、c₁代表直线方程L₁的系数；

步骤E3、计算权重序列点[A₁,A₂…A_r]到L₁的距离定义如下：

${\mathrm{DT}}_{L_1}\left(i\right)=\frac{\left|a_{1}r+b_1{w}_i+c_1\right|}{\sqrt{{a_1}^2+{b_1}^2}},i=1,2,3...r---\left(11\right)$

步骤E4、搜寻边界距离的极大值：

$p=S\left(\underset{i}{max}\right({\mathrm{DT}}_{L1}\left(i\right)\left)\right)-1,(i=1,2...r)---\left(12\right)$

步骤E5、对于第一次判别p大于4的序列，重新将A₂与Ar之间的连线L₂作为边界直线，则

iffloor(p/2)>2

L₂:a₂x+b₂y+c₂＝0(13)

其中floor()代表向下取整运算；a₂、b₂、c₂代表直线方程L₂的系数；

对于第一次判别p小于等于4的序列，跳过E5-E7；

步骤E6、计算权重序列点[A₂,A₃…A_r]到L₂的距离：

${\mathrm{DT}}_{L_2}\left(i\right)=\frac{\left|a_{2}r+b_2{w}_i+c_2\right|}{\sqrt{{a_2}^2+{b_2}^2}},i=2,3...r---\left(14\right)$

步骤E7、搜寻边界距离的极大值，找出对应的序列值：

$p=S\left(\underset{i}{max}\right({\mathrm{DT}}_{L2}\left(i\right)\left)\right)-1,(i=2,3...r)---\left(15\right).$