CN102622608B - 基于非线性时间序列分析的高光谱图像分类方法 - Google Patents

Abstract

基于非线性时间序列分析的新型高光谱图像分类方法，即通过一种非线性时间序列分析处理高光谱反射率曲线、从而对高光谱图像中不同像元进行特征构造、然后依据构造的特征完成分类的方法。该方法如下：第一，通过人机交互接口模块获得待处理的高光谱数据；第二，通过高光谱特征构造模块获得用于地物分类的特征组合；第三，通过高光谱地物分类模块，使用特征构造结果对实例进行地物分类；第四，通过分类结果输出模块，输出分类结果。本发明的优点在于鲁棒性强，空间复杂度小，而时间复杂度与样本点数保持线性关系，分类精度高，适用范围广。

Description

基于非线性时间序列分析的高光谱图像分类方法

技术领域

本发明涉及一种基于非线性时间序列分析的新型高光谱图像分类方法，适用于高光谱图像处理系统中，属于高光谱图像处理领域。

背景技术

20世纪80年代遥感技术的最大成就之一是高光谱遥感技术的兴起，由于兼具成像和光谱探测的优点，该技术已经在商业、军事和民间领域得到广泛应用。随着高光谱影像光谱分辨率的提高，在常规遥感中无法探测的物质在高光谱遥感中可以被探测，这为后期进行地物详细分类提供了前提条件。尽管高光谱可以提供丰富的地面信息，但其大量的波段数导致了呈指数级增长的数据量和信息冗余，从而增加了后期数据处理的复杂度。

在地物识别或分类时，并非使用的波段数越多越好。这是因为：

1)不能把光谱通道数简单等同于信息通道数。因为高光谱图像相邻的波段之间具有较大的相关性，如果不加分析地使用所有谱段进行识别或分类，不仅不会提高分类精度，反而影响识别或分类结果。

2)并不是所有的波段对于后续处理都有着同等的重要性，通过线性或者非线性变换得到新的图像空间并保留不同目标地物类别间有显著差异的特征，在不损失重要信息的条件下可以代表其它波段的信息。

3)选择谱段过多，不仅增加运算量，影响计算速度，而且还需要大量的训练样本，否则难以获得预期的识别或分类结果。针对高光谱数据高维、海量信息需在合理的数学模型指导下，按照一定的准则选取最佳波段组合，进行数据降维，是提升高光谱数据处理系统处理速度与精度的有效途径。

综上所述，高光谱遥感图像的优越性是以其大的数据量以及高的数据维为代价的，因此传统的遥感图像处理技术在这里将不再适用，一些针对传统遥感数据的图像处理方法和技术面临挑战，需要研究适合于高光谱图像的处理方法。

时间序列分析是一种动态数据处理的经典方法。该方法基于随机过程理论和数理统计学方法，研究随机数据序列所遵从的统计规律，以用于解决实际问题。狭义上讲，时间序列是按时间顺序的一组数字序列；广义上讲，时间序列包括按照一定顺序排列的数字序列。时间序列分析就是利用这组数列，应用数理统计方法加以处理，以预测未来事物的发展。它一般采用曲线拟合和参数估计方法，其优点是既承认了观测数据的依存性，又考虑了随机波动的不定性。时间序列预测一般反映三种实际变化规律：趋势性变化、周期性变化、随机性变化。时间序列分析常用在国民经济宏观控制、区域综合发展规划、企业经营管理、市场潜量预测、工业自动化生产、气象预报、水文预报、地震前兆预报、农作物病虫灾害预报、环境污染控制、生态平衡、天文学和海洋学等方面。

然而，传统的时间序列分析方法存在如下假设：

1)时间序列变量是从某个随机过程中随机抽取并按一定顺序排列而成的，因而存在一个稳定趋势。

2)时间序列变量的波动率或方差是固定的，即等于一个常数。

这些假设显然不符合实际，极大限制了传统时间序列分析方法在实际中的应用，而混沌时间序列分析则可以有效规避这些限制。混沌现象是非周期的具有渐近自相似有序性的现象，是确定性非线性系统由于系统内部非线性相互作用所产生的貌似无规则的、类似随机运动的现象。但混沌不是简单的无序，虽然它没有明显的周期和对称，却具有丰富的内部层次，蕴涵着深层次的结构和秩序，并在一定条件下以某种方式显示出来。混沌理论试图建立一个很小或几乎没有随机误差的非线性复杂模型，在分析时主要考虑模型本身对时间序列未来值的影响，从而能够预测时间序列微小起伏波动的精细程度。混沌时间序列分析方法已经在众多领域取得成果，本技术发明把该分析方法引入高光谱图像处理中并加以改进。

在高光谱数据处理系统中，出于计算复杂度要求，分类算法要具有小的计算量，而到达此要求的传统算法分类精度差。出于分类精度的要求，分类算法要具有对不同分类场景的良好的健壮性，而到达此要求的算法计算复杂度高。因此，需要寻找一种综合考虑计算复杂度与分类精度这两方面平衡的算法，使其具有低的时间复杂度及良好的鲁棒性。

发明内容

本发明的目的在于提供一种基于非线性时间序列分析的新型高光谱图像分类方法，发明了一种运用非线性时间序列分析处理高光谱反射率曲线、从而对高光谱图像中不同像元进行特征构造、进而依据构造的特征完成分类的方法。本方法保持计算量小的同时又保持了分类精度。此方法应用到高光谱数据处理系统中具有强的鲁棒性。

本发明方法基于仿真原型系统，该系统具有人机交互接口模块、高光谱特征构造模块、高光谱地物分类模块、分类结果输出模块等四个功能模块。其中，人机交互接口模块完成高光谱数据的读入与相关参数的设定；高光谱特征构造模块完成运用混沌时间序列分析对高光谱数据进行维数约简；高光谱地物分类模块完成应用特征构造结果对高光谱数据的不同地物进行分类；分类结果输出模块完成输出对高光谱数据不同地物的分类结果。

本发明所涉及的方法流程包括以下步骤：（1）获得初始数据及相关初始化操作；（2）计算每个纯像元高光谱反射率曲线的Hurst指数；（3）计算每个纯像元高光谱反射率曲线的最大Lyapunov指数；（4）分类；（5）结果输出。其中，人机交互接口模块对应（1），高光谱特征构造模块对应（2）～（3），高光谱地物分类模块对应（4），分类结果输出模块对应（5）。

下面对该方法流程各步骤进行详细说明：

对任一高光谱像元X，记为X=(x₁,…,x_i,…,x_K)，其中K是高光谱数据的波段总数，x_i,i=1,…,K代表第i个波段的响应数值。Hurst指数记为

，最大Lyapunov指数记为

。

通过人机交互接口模块完成步骤一。

步骤一获得初始数据及相关初始化操作

通过人机交互接口模块获得待处理的高光谱数据。

通过高光谱特征构造模块完成步骤二～步骤三。

步骤二计算每个纯像元高光谱反射率曲线的Hurst指数

将X=(x₁,…,x_i,…,x_K)按照长度n划分为D个子区间，第d个子区间X_d如公式(1)所示：

X_d={x_(d-1)×n+1,…,x_d×n}                  (1)

其中，d=1,…,D，n是每个子区间的长度，是介于lnK到

之间的正整数。

对第d个子区间X_d的各个元素按照公式(2)重命名为g_d,i：

g_d,i=x_(d-1)×n+i                       (2)

其中，i=1,…,n，d=1,…,D。

按照公式(3)计算第d个子区间X_d的平均值m_d和标准差S_d：

$\left\{\begin{array}{c}{m}_d=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{g}_{d,i}\\ S_d=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(g_{d,i}-m_d)}^2}\end{array}\right.---\left(3\right)$

其中，i=1,…,n，d=1,…,D。对第d个子区间X_d的各个元素g_d,i减去第d个子区间X_d的平均值m_d，即g*_d,i=g_d,i-m_d，其中，i=1,…,n，d=1,…,D。

按照公式(4)构造g*_d,i的累积和序列{c_d,i}：

$c_{d,i}=\underset{j=1}{\overset{i}{\mathrm{\Σ}}}g{*}_{d,j}---\left(4\right)$

其中，i=1,…,n，d=1,…,D。按照公式(5)得到第d个子区间X_d关于累积和序列{c_d,i}的极差R_d：

$R_d=\underset{1\≤i\≤n}{\max }\left\{c_{d,i}\right\}-\underset{1\≤i\≤n}{\min }\left\{c_{d,i}\right\}---\left(5\right)$

其中，d=1,…,D。

按照公式(6)计算所有长度为n的子区间的重标极差R_d/S_d的平均值(R/S)_n：

${(R/S)}_n=\frac{1}{D}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{R}_d/S_d---\left(6\right)$

对所有取值从lnK取到

的n重复以上操作。对ln(R/S)_n作关于lnn的最小二乘回归，回归直线的斜率就是

的计算值，其中n是介于lnK到

之间的正整数。

步骤三计算每个纯像元高光谱反射率曲线的最大Lyapunov指数

给定平均周期p、延迟时间τ和嵌入维数m，将X=(x₁,…,x_i,…,x_K)重构，以得到相空间{Y_j}。其中，Y_j=(x_j,x_j+τ,…,x_j+(m-1)τ)，j=1,…,M，M=K-(m-1)τ。M是相空间中的点数，即Y_j的个数。

对相空间中的每一个点Y_j，j=1,…,M，寻找满足以下两个条件的临近点Y_j'：

1)|j'-j|>p，即限制短暂分离；

2)Y_j'与Y_j的欧式距离最小。

根据相空间中点的下标的先后次序，计算临近点对Y_j'与Y_j的演化距离d_j(i)=||Y_j'+i-Y_j+i||，其中i=1,…,min(M-j',M-j)是演化步长。

记所有临近点对的最大演化步长是I₀，依次给定演化步长i=1,…,I₀，设当前演化步长i下非空且非零的d_j(i)有q个，将这q个非空且非零的d_j(i)记为

计算 $y\left(i\right)=\frac{1}{q}\underset{j=1}{\overset{q}{\mathrm{\Σ}}}\ln {\stackrel{\‾}{d}}_j\left(i\right),$ i＝1，…，I₀。

对{y(i)}作关于i=1,…,I₀的最小二乘回归，回归直线的斜率就是

的计算值。

任一高光谱像元的维数约简结果由其自身的Hurst指数和最大Lyapunov指数

组成。

通过高光谱地物分类模块完成步骤四。

步骤四分类

使用维数约简结果对高光谱数据进行地物分类。

通过分类结果输出模块完成步骤五。

步骤五结果输出

输出高光谱图像地物分类结果。

本发明是一种基于混沌时间序列分析的高光谱图像分类方法，其优点在于：用于高光谱数据处理系统中，特征构造结果不受其它像元的干扰，鲁棒性强，空间复杂度小，而时间复杂度与样本点数保持线性关系，分类精度高，适用范围广。

具体实施方式

下面结合实例进一步说明本发明的应用方法。

基于本发明开发了仿真原型系统，该系统包括：人机交互接口模块、高光谱特征构造模块、高光谱地物分类模块、分类结果输出模块等四个功能模块。

第一，通过人机交互接口模块获得待处理的高光谱数据。本实例采用Washington D.C.Mall高光谱数据，大小为1280×307像素，波长范围为0.4～2.4mm，除去水汽吸收波段和低信噪比波段后，保留191个波段，并截取其中一子图大小为562×307像素，子图共包含7类地物，分别是：屋顶、草地、树木、小路、街道、水、阴影。

第二，通过高光谱特征构造模块获得用于地物分类的特征组合，即对实例的每个像元依次进行如下处理。

1)计算每个纯像元高光谱反射率曲线的Hurst指数

设定子区间长度n的范围介于ln191与

之间的正整数，即n=6,7,…,12,13。

n从6开始，将像元的191维向量按照每6个一组划分为个子区间，其中

是向下取整函数，第d个子区间记为X_d。

计算X_d的平均值m_d和标准差S_d，将X_d的每个元素减去平均值m_d，再构造对应的累积和序列{c_d,i}，从而得到累积和序列{c_d,i}的极差R_d，最后计算第d个子区间X_d的重标极差R_d/S_d。计算当n=6时所有子区间的重标极差，再求其平均，得到(R/S)_n。

重复以上操作，直至遍历所有合适的n。对ln(R/S)_n作关于lnn的最小二乘回归，回归直线的斜率就是

的计算值，即的第一个元素，其中

2)计算每个纯像元高光谱反射率曲线的最大Lyapunov指数

取定平均周期p、延迟时间τ和嵌入维数m分别为p=21,τ=3,m=8。则相空间一共有M=191-(8-1)×3=170个点，其中Y₁=(x₁,x₄,x₇,…,x₂₂)，Y₂=(x₂,x₅,x₈,…,x₂₃)，……，Y₁₇₀=(x₁₇₀,x₁₇₃,x₁₇₆,…,x₁₉₁)。

假设当前像元对应Y₁的满足限制短暂分离的临近点是Y₈₀，则此时演化步长i的最大允许值是min(170-80,170-1)=90，否则将造成相空间溢出。d₁(1)=||Y₈₁-Y₁||，d₁(2)=||Y₈₂-Y₂||，d₁(3)=||Y₈₃-Y₃||，……，d₁(90)=||Y₁₉₁-Y₉₁||，而d₁(91),d₁(92),…均无意义。对相空间中的170个点重复上述操作。

记所有临近点对的最大演化步长是I₀=100。设当前演化步长i=1下非空且非零的d_j(1)有q=20个，将这20个非空且非零的d_j(i)记为计算

依次给定演化步长i=1,…,100重复上述操作。

对{y(i)}作关于i的最小二乘回归，回归直线的斜率就是

的计算值，即

的第一个元素，其中

第三，通过高光谱地物分类模块，使用特征构造结果对实例进行地物分类。

第四，通过分类结果输出模块，输出对实例的分类结果。

本发明方法是一种基于混沌时间序列分析的高光谱图像分类方法，经过仿真系统的具体实施，能够有效避免维数灾难和豪斯现象的干扰，在保留目标地物主要特征的条件下完成高光谱图像的监督分类。用于高光谱数据处理系统中，特征构造结果不受其它像元的干扰，鲁棒性强，空间复杂度小，而时间复杂度与样本点数保持线性关系，分类精度高，适用范围广。

Claims (1)

1.基于非线性时间序列分析的高光谱图像分类方法，其特征在于：基于仿真原型系统，具有人机交互接口模块、高光谱特征构造模块、高光谱地物分类模块、分类结果输出模块这四个功能模块，包括如下步骤：

假设X为任一高光谱像元，记为X=(x₁,…,x_i,…,x_K)，其中K是高光谱数据的波段总数，x_i,i=1,…,K代表第i个波段的响应数值，Hurst指数记为

，最大Lyapunov指数记为

；

通过人机交互接口模块完成步骤一

步骤一获得初始数据及相关初始化操作

通过人机交互接口模块获得待处理的高光谱数据；

通过高光谱特征构造模块完成步骤二～步骤三

步骤二计算每个纯像元高光谱反射率曲线的Hurst指数

将X=(x₁,…,x_i,…,x_K)按照长度n划分为D个子区间，第d个子区间X_d如公式(1)所示：

X_d={x_(d-1)×n+1,…,x_d×n}                     (1)

其中，d=1,…,D，n是每个子区间的长度，是介于从lnK到之间的正整数；

对第d个子区间X_d的各个元素按照公式(2)重命名为g_d,i：

g_d,i=x_(d-1)×n+i                             (2)

其中，i=1,…,n，d=1,…,D；

按照公式(3)计算第d个子区间X_d的平均值m_d和标准差S_d：

$\left\{\begin{array}{c}{m}_d=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{g}_{d,i}\\ S_d=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(g_{d,i}-m_d)}^2}\end{array}\right.---\left(3\right)$

其中，i=1,…,n，d=1,…,D，对第d个子区间X_d的各个元素g_d,i减去第d个子区间X_d的平均值m_d，即g*_d,i=g_d,i-m_d，其中，i=1,…,n，d=1,…,D；

按照公式(4)构造g*_d,i的累积和序列{c_d,i}：

$c_{d,i}=\underset{j=1}{\overset{i}{\mathrm{\Σ}}}g{*}_{d,j}---\left(4\right)$

其中，i=1,…,n，d=1,…,D；

按照公式(5)得到第d个子区间X_d关于累积和序列{c_d,i}的极差R_d：

$R_d=\underset{1\≤i\≤n}{\max }\left\{c_{d,i}\right\}-\underset{1\≤i\≤n}{\min }\left\{c_{d,i}\right\}---\left(5\right)$

其中，d=1,…,D；

按照公式(6)计算所有长度为n的子区间的重标极差R_d/S_d的平均值(R/S)_n：

${(R/S)}_n=\frac{1}{D}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{R}_d/S_d---\left(6\right)$

对所有取值从lnK取到的n重复以上操作，对ln(R/S)_n作关于lnn的最小二乘回归，回归直线的斜率就是的计算值，其中n是介于lnK到之间的正整数；

步骤三计算每个纯像元高光谱反射率曲线的最大Lyapunov指数

给定平均周期p、延迟时间τ和嵌入维数m，将X=(x₁,…,x_i,…,x_K)重构，以得到相空间{Y_j},其中，Y_j=(x_j,x_j+τ,…,x_j+(m-1)τ)，j=1,…,M，M=K-(m-1)τ，M是相空间中的点数，即Y_j的个数；

对相空间中的每一个点Y_j，j=1,…,M，寻找满足以下两个条件的临近点Y_j'：

1)|j'-j|>p，即限制短暂分离；

2)Y_j'与Y_j的欧式距离最小；

根据相空间中点的下标的先后次序，计算临近点对Y_j'与Y_j的演化距离d_j(i)=||Y_j'+i-Y_j+i||，其中i=1,…,min(M-j',M-j)是演化步长；

记所有临近点对的最大演化步长是I₀，依次给定演化步长i=1,…,I₀，设当前演化步长i下非空且非零的d_j(i)有q个，将这q个非空且非零的d_j(i)记为计算 $y\left(i\right)=\frac{1}{q}\underset{j=1}{\overset{q}{\mathrm{\Σ}}}\ln {\stackrel{\‾}{d}}_j\left(i\right),$ i=1,…,I₀；

对{y(i)}作关于i=1,…,I₀的最小二乘回归，回归直线的斜率就是

的计算值；

任一高光谱像元的维数约简结果由其自身的Hurst指数

和最大Lyapunov指数

组成；

通过高光谱地物分类模块完成步骤四

步骤四分类

使用维数约简结果对高光谱数据进行地物分类；

通过分类结果输出模块完成步骤五

步骤五结果输出

输出高光谱图像地物分类结果。