基于领域与l2,1范数的光谱图像分类方法及系统与流程

本发明涉及光谱图像分类领域，具体地说，本发明涉及一种基于领域与l_2,1范数的光谱图像分类技术方案。

背景技术：

在过去的几十年中，由于高光谱(通常为200波段左右)成像能够同时提供丰富的空间和光谱信息，高光谱成像已经广泛地应用在遥感领域中。而高光谱图像分类是高光谱成像很常见的应用之一，各种各样的分类技术已经应用于高光谱成像的实际应用中，如目标识别、土地利用分析和环境污染监测等。

在已经提出的分类方法中，最具代表性的方法之一是支持向量机(SVM)，SVM即使在少量训练数据的情况下也取得了令人满意的分类效果。近年来，基于稀疏表示的分类方法(SRC)已经在高光谱图像分类中得到了广泛地关注。每个测试像素由带标签的训练数据稀疏表示，测试数据的类标签由表示误差最小的类决定。在文献[1]中，为了将高光谱图像的空间信息考虑进去，陈等人提出了同时正交匹配追踪(SOMP)方法用于高光谱图像分类，每个测试像素的邻域类的所有像素由带标签的训练数据共同稀疏表示。在文献[2]中，李等人提出了联合鲁棒稀疏分类(JRSRC)方法将稀疏表示剩余(即异常)考虑进去，它能处理高光谱图像分类中的异常。在文献[3]中，陈等人将联合稀疏表示方法扩展到核稀疏表示方法，并且用到高光谱图像分类中。此外，张等人[4]提出了一种基于协作表示的分类框架(CRC)用于分类中，并且取得了很好的分类效果，且时间消耗要低于基于稀疏表示的分类方法。

尽管上述的基于稀疏表示和协作表示方法已经在高光谱图像分类中取得了不错的效果，但都没有同时将所有的测试像素考虑进去。

[1]Y.Chen,N.M.Nasrabadi,and T.D.Tran,“Hyperspectral image classification using dictionary-based sparse representation,”IEEE Trans.Geosci.Remote Sens.,vol.49,no.10,pp.3973–3985,2011.

[2]C.Li,Y.Ma,X.Mei,C.Liu,and J.Ma,“Hyperspectral image classification with robust sparse representation,”IEEE Geosci.Remote Sens.Lett.,vol.13,no.5,pp.641–645,2016.

[3]Y.Chen,N.M.Nasrabadi,and T.D.Tran,“Hyperspectral image classification via kernel sparse representation,”IEEE Trans.Geosci.Remote Sens.,vol.51,no.1,pp.217–231,Jan.2013.

[4]L.Zhang,M.Yang,and X.Feng,“Sparse representation or collaborative representation:Which helps face recognition？”in IEEE International Conference on Computer Vision,2011,pp.471–478.

技术实现要素：

为克服相应技术缺陷，本发明提出了一种基于领域与l_2,1范数的高光谱图像分类技术方案方案。

本发明技术方案提供一种基于领域与l_2,1范数的高光谱图像分类方法，包括以下步骤：

步骤1，令光谱图像为T∈R^K×P×B，其中K和P为图像在空间维度上的高度和宽度，B为光谱维度上的波段数，R为实数，T中每个像素对应的光谱记为T_i,j，i＝1,2…K，j＝1,2…P，对光谱图像的每个像素都进行领域平均操作，得到领域平均操作后的光谱图像Y∈R^K×P×B；

步骤2，建立光谱图像分类的数学模型，得到相应最优化问题，实现如下，

设有M种待分类目标类别，令A＝[a₁,a₂,...,a_M]∈R^B×M表示M个端元的光谱，则光谱图像Y中各像素对应的光谱Y_1,1,Y_1,2,…Y_1,P,Y_2,1,…,Y_K,P的数学模型如下，

C＝[Y_1,1,Y_1,2,…Y_1,P,Y_2,1,…,Y_K,P]＝AX+N

其中，X∈R^M×KP是丰度系数矩阵，N∈R^B×KP表示误差矩阵，C表示光谱图像Y中从左到右从上到下的像素对应光谱按列排列成的矩阵；

则光谱图像上各像素对应光谱的分类转化为如下最优化问题，

$\underset{X\≥0}{min}\left|\right|AX-C|{|}_{2,1}+\mathrm{\λ}|\left|X\right|{|}_{2,1}$

上式表示求令目标函数||AX-C||_2,1+λ||X||_2,1关于X的最小值；其中，min是最小化算子，λ＞0表示正则化参数，||·||_2,1表示计算l_2,1范数，定义如下，

$\left|\right|AX-C|{|}_{2,1}={\mathrm{\Σ}}_{b=1}^B\sqrt{{\mathrm{\Σ}}_{p=1}^{KP}{N}_{b,p}^2}$

$\left|\right|X|{|}_{2,1}={\mathrm{\Σ}}_{m=1}^M\sqrt{{\mathrm{\Σ}}_{p=1}^{KP}{X}_{m,p}^2}$

其中，N_b,p为误差矩阵N坐标(b,p)处数值，X_m,p为丰度系数矩阵X坐标(m,p)处数值；

步骤3，求解步骤2所得最优化问题，得到令目标函数|AX-C||_2,1+λ||X||_2,1取最小值的解X；

步骤4，根据步骤3求解到的X对光谱图像Y进行分类，

光谱图像Y中各像素Y_i,j的类别Class(Y_i,j)由最小类重构误差得到如下，

$Class\left(Y_{i,j}\right)=\arg \underset{m=1,...,M}{min}\left|\right|Y_{i,j}-a_m{X}_{m,i\×P+j}\left|\right|,i=1,2...K;j=1,2...P$

上式表示如果1到M中的某个值m令||Y_i,j-a_mX_m,i×P+j||取到最小值，则Y_i,j就归属到端元库中的第m个类别；X_m,i×P+j表示丰度系数矩阵X坐标(m,i×P+j)处数值。

而且，步骤1中，进行领域平均操作的实现方式如下，

在以T_i,j为中心且大小为D×D的窗口Q_i,j内进行领域平均操作，计算平均结果Y_i,j如下，

$Y_{i,j}=\frac{1}{D^2}\underset{T_Q\∈Q_{i,j}}{\Σ}{T}_Q$

其中，(D+1)/2≤i≤K-(D-1)/2，(D+1)/2≤j≤P-(D-1)/2；

对于光谱图像T边界(D-1)/2范围内的像素，采用边界扩展技术后再进行领域平均操作。

而且，步骤3中，采用交替方向乘子法求解最优化问题。

本发明提供一种基于领域与l_2,1范数的光谱图像分类系统，包括以下模块，

领域平均模块，用于令光谱图像为T∈R^K×P×B，其中K和P为图像在空间维度上的高度和宽度，B为光谱维度上的波段数，R为实数，T中每个像素对应的光谱记为T_i,j，i＝1,2…K，j＝1,2…P，对光谱图像的每个像素都进行领域平均操作，得到领域平均操作后的光谱图像Y∈R^K×P×B；模型构建模块，用于建立光谱图像分类的数学模型，得到相应最优化问题，实现如下，

设有M种待分类目标类别，令A＝[a₁,a₂,...,a_M]∈R^B×M表示M个端元的光谱，则光谱图像Y中各像素对应的光谱Y_1,1,Y_1,2,…Y_1,P,Y_2,1,…,Y_K,P的数学模型如下，

C＝[Y_1,1,Y_1,2,…Y_1,P,Y_2,1,…,Y_K,P]＝AX+N

其中，X∈R^M×KP是丰度系数矩阵，N∈R^B×KP表示误差矩阵，C表示光谱图像Y中从左到右从上到下的像素对应光谱按列排列成的矩阵；

则光谱图像上各像素对应光谱的分类转化为如下最优化问题，

$\underset{X\≥0}{min}\left|\right|AX-C|{|}_{2,1}+\mathrm{\λ}|\left|X\right|{|}_{2,1}$

上式表示求令目标函数||AX-C||_2,1+λ||X||_2,1关于X的最小值；其中，min是最小化算子，λ＞0表示正则化参数，||·||₂₁表示计算l_2,1范数，定义如下，

$\left|\right|AX-C|{|}_{2,1}={\mathrm{\Σ}}_{b=1}^B\sqrt{{\mathrm{\Σ}}_{p=1}^{KP}{N}_{b,p}^2}$

$\left|\right|X|{|}_{2,1}={\mathrm{\Σ}}_{m=1}^M\sqrt{{\mathrm{\Σ}}_{p=1}^{KP}{X}_{m,p}^2}$

其中，N_b,p为误差矩阵N坐标(b,p)处数值，X_m,p为丰度系数矩阵X坐标(m,p)处数值；

求解模块，用于求解最优化问题，得到令目标函数|AX-C||_2,1+λ||X||_2,1取最小值的解X；

分类模块，用于根据求解到的X对光谱图像Y进行分类，

光谱图像Y中各像素Y_i,j的类别Class(Y_i,j)由最小类重构误差得到如下，

$Class\left(Y_{i,j}\right)=\arg \underset{m=1,...,M}{min}\left|\right|Y_{i,j}-a_m{X}_{m,i\×P+j}\left|\right|,i=1,2...K;j=1,2...P$

上式表示如果1到M中的某个值m令||Y_i,j-a_mX_m,i×P+j||取到最小值，则Y_i,j就归属到端元库中的第m个类别；X_m,i×P+j表示丰度系数矩阵X坐标(m,i×P+j)处数值。

而且，领域平均模块中，进行领域平均操作的实现方式如下，

在以T_i,j为中心且大小为D×D的窗口Q_i,j内进行领域平均操作，计算平均结果Y_i,j如下，

$Y_{i,j}=\frac{1}{D^2}\underset{T_Q\∈Q_{i,j}}{\Σ}{T}_Q$

其中，(D+1)/2≤i≤K-(D-1)/2，(D+1)/2≤j≤P-(D-1)/2；

对于光谱图像T边界(D-1)/2范围内的像素，采用边界扩展技术后再进行领域平均操作。

而且，求解模块中，采用交替方向乘子法求解最优化问题。

本发明提出了一种基于联合领域与l_2,1范数的光谱图像分类技术方案用于高光谱图像分类，同时将所有的测试像素同时考虑进去，比单独地对每个测试像素进行稀疏恢复具有更好的效果。此外，还采用了l_2,1范数的损失函数来使算法对异常鲁棒，并通过ADMM进行求解。本发明技术方案具有分类准确度高的优点，也能适用于超光谱图像，市场价值高。

附图说明

图1是本发明实施例的流程图。

具体实施方式

下面结合附图和实施例对本发明进行进一步的说明。

参照附图1，本发明实施例主要由4个步骤组成：对光谱图像的每个像素都进行领域平均操作，建立光谱图像分类的数学模型，求解最优化模型，根据求解到的X对光谱图像进行分类。实施例选取的真实数据是Indian Pines高光谱数据集，是在1992年通过Airborne Visible/Infrared Imaging Spectrometer(AVIRIS)传感器在印第安纳州的西北部采集得到。图像大小为145×145，共有220个波段，覆盖0.4–2.5微米的光谱范围。去除水汽吸收的波段(104–108,150–163和220)后剩下200个波段，总共有16个类别的端元。

具体实施时，本发明技术方案可采用计算机软件技术实现自动运行流程。实施例执行步骤如下：

步骤1，对光谱图像的每个像素都进行领域平均操作；

令光谱图像为T∈R^K×P×B，其中K和P为图像在空间维度上的高度和宽度，B为光谱维度上的波段数，R为实数，则T中每个像素对应的光谱记为T_i,j∈R^B×1(i＝1,2…K；j＝1,2…P)，在以T_i,j为中心且大小为D×D(D为奇数，具体实施时本领域技术人员可预先设置，建议取值大于3小于图像长或者宽的1/3)的窗口Q_i,j内进行领域平均操作，计算平均结果Y_i,j方法如下：

$Y_{i,j}=\frac{1}{D^2}\underset{T_Q\∈Q_{i,j}}{\Σ}{T}_Q$

其中，考虑到边界，可设(D+1)/2≤i≤K-(D-1)/2，(D+1)/2≤j≤P-(D-1)/2，即Q_i,j为光谱图像T边界(D-1)/2范围外的像素，T_Q∈Q_i,j表示属于窗口Q_i,j内像素对应的光谱；特殊地，对于光谱图像T边界(D-1)/2范围内的像素，采用边界扩展技术后再进行领域平均操作。；最后得到领域平均操作后的光谱图像Y∈R^K×P×B；

边界扩展技术为现有技术，本发明不予赘述，可参见文献Samnotra M,Girdhar A.Ultrasound Image Enhancement Using Laplacian Kernel Set[J].Ultrasound,2015,4(12).

步骤2，建立光谱图像分类的数学模型；

令A＝[a₁,a₂,...,a_M]∈R^B×M表示M个端元(即M种待分类目标类别)的光谱，则光谱图像Y中各像素对应的光谱Y_1,1,Y_1,2,…Y_1,P,Y_2,1,…,Y_K,P可以由A的线性组合来表示，数学模型如下：

C＝[Y_1,1,Y_1,2,…Y_1,P,Y_2,1,…,Y_K,P]＝AX+N

其中，X∈R^M×KP是丰度系数矩阵，N∈R^B×KP表示误差矩阵，C表示光谱图像Y中从左到右从上到下的像素对应光谱按列排列成的新矩阵；

则光谱图像上各像素对应光谱的分类可以转化为如下最优化问题：

$\underset{X\≥0}{min}\left|\right|AX-C|{|}_{2,1}+\mathrm{\λ}|\left|X\right|{|}_{2,1}$

上式表示求令目标函数||AX-C||_2,1+λ||X||_2,1关于X的最小值；其中，min是最小化算子，λ＞0表示正则化参数，具体实施时本领域技术人员可预设取值，||·||_2,1表示计算l_2,1范数，具体定义如下：

$\left|\right|AX-C|{|}_{2,1}={\mathrm{\Σ}}_{b=1}^B\sqrt{{\mathrm{\Σ}}_{p=1}^{KP}{N}_{b,p}^2}$

$\left|\right|X|{|}_{2,1}={\mathrm{\Σ}}_{m=1}^M\sqrt{{\mathrm{\Σ}}_{p=1}^{KP}{X}_{m,p}^2}$

其中，N_b,p为误差矩阵N坐标(b,p)处数值，X_m,p为丰度系数矩阵X坐标(m,p)处数值。

步骤3，求解最优化模型；

采用交替方向乘子法(ADMM)求解上述最优化问题，得到令目标函数|AX-C||_2,1+λ||X||_2,1取最小值的解X；

交替方向乘子法为现有技术，本发明不予赘述，可参见文献S.Boyd,N.Parikh,E.Chu,B.Peleato,and J.Eckstein,“Distributed optimization and statistical learning via the alternating direction method of multipliers,”Foundations and Trends R in Machine Learning,vol.3,no.1,pp.1–122,2011.

步骤4，根据求解到的X对光谱图像进行分类；

光谱图像Y中各像素Y_i,j的类别Class(Y_i,j)由最小类重构误差得到：

$Class\left(Y_{i,j}\right)=\arg \underset{m=1,...,M}{min}\left|\right|Y_{i,j}-a_m{X}_{m,i\×P+j}\left|\right|,i=1,2...K;j=1,2...P$

上式表示如果1到M中的某个值m令||Y_i,j-a_mX_m,i×P+j||取到最小值，则Y_i,j就归属到端元库中的第m个类别；X_m,i×P+j表示丰度系数矩阵X坐标(m,i×P+j)处数值，即像素Y_i,j对应光谱在第m个类别物质的丰度系数。

实施例中，λ＝5*10^-6，D＝36，M＝16，B＝200，P＝K＝145。为便于理解本发明实施例技术方案的技术效果，选取了7种不同的方法做对比，即SVM[1]、NRS[2]、Gabor-NRS[3]、CRT[4]、OMP[5]、SOMP[5]和JRSRC[6]。采用总准确率OA、平均准确率AA来评价这些算法的分类效果。可以看到本发明在总准确率与平均准确率上都是最优的。

表1：对Indian Pines数据集采用不同的分类方法得到的OA(％)、AA(％)

[1]R.Archibald and G.Fann,“Feature selection and classification of hyperspectral images with support vector machines,”IEEE Geosci.Remote Sens.Lett.,vol.4,no.4,pp.674–677,2007.

[2]R.Archibald and G.Fann,“Feature selection and classification of hyperspectral images with support vector machines,”IEEE Geosci.Remote Sens.Lett.,vol.4,no.4,pp.674–677,2007.

[3]W.Li and Q.Du,“Gabor-filtering-based nearest regularized subspace for hyperspectral image classification,”IEEE J.Sel.Topics Appl.Earth Observ.Remote Sens.,vol.7,no.4,pp.1012–1022,2014.

[4]W.Li,Q.Du,and M.Xiong,“Kernel collaborative representation with tikhonov regularization for hyperspectral image classification,”IEEE Geosci.Remote Sens.Lett.,vol.12,no.1,pp.48–52,2015.

[5]Y.Chen,N.M.Nasrabadi,and T.D.Tran,“Hyperspectral image classification using dictionary-based sparse representation,”IEEE Trans.Geosci.Remote Sens.,vol.49,no.10,pp.3973–3985,2011.

[6]C.Li,Y.Ma,X.Mei,C.Liu,and J.Ma,“Hyperspectral image classification with robust sparse representation,”IEEE Geosci.Remote Sens.Lett.,vol.13,no.5,pp.641–645,2016.

具体实施时，本发明所提供方法可基于软件技术实现自动运行流程，也可采用模块化方式实现相应系统。

本发明提供一种基于领域与l_2,1范数的光谱图像分类系统，包括以下模块，

领域平均模块，用于令光谱图像为T∈R^K×P×B，其中K和P为图像在空间维度上的高度和宽度，B为光谱维度上的波段数，R为实数，T中每个像素对应的光谱记为T_i,j，i＝1,2…K，j＝1,2…P，对光谱图像的每个像素都进行领域平均操作，得到领域平均操作后的光谱图像Y∈R^K×P×B；模型构建模块，用于建立光谱图像分类的数学模型，得到相应最优化问题，实现如下，

设有M种待分类目标类别，令A＝[a₁,a₂,...,a_M]∈R^B×M表示M个端元的光谱，则光谱图像Y中各像素对应的光谱Y_1,1,Y_1,2,…Y_1,P,Y_2,1,…,Y_K,P的数学模型如下，

C＝[Y_1,1,Y_1,2,…Y_1,P,Y_2,1,…,Y_K,P]＝AX+N

其中，X∈R^M×KP是丰度系数矩阵，N∈R^B×KP表示误差矩阵，C表示光谱图像Y中从左到右从上到下的像素对应光谱按列排列成的矩阵；

则光谱图像上各像素对应光谱的分类转化为如下最优化问题，

$\underset{X\≥0}{min}\left|\right|AX-C|{|}_{2,1}+\mathrm{\λ}|\left|X\right|{|}_{2,1}$

上式表示求令目标函数||AX-C||_2,1+λ||X||_2,1关于X的最小值；其中，min是最小化算子，λ＞0表示正则化参数，||·||_2,1表示计算l_2,1范数，定义如下，

$\left|\right|AX-C|{|}_{2,1}={\mathrm{\Σ}}_{b=1}^B\sqrt{{\mathrm{\Σ}}_{p=1}^{KP}{N}_{b,p}^2}$

$\left|\right|X|{|}_{2,1}={\mathrm{\Σ}}_{m=1}^M\sqrt{{\mathrm{\Σ}}_{p=1}^{KP}{X}_{m,p}^2}$

其中，N_b,p为误差矩阵N坐标(b,p)处数值，X_m,p为丰度系数矩阵X坐标(m,p)处数值；

求解模块，用于求解最优化问题，得到令目标函数|AX-C||_2,1+λ||X||_2,1取最小值的解X；

分类模块，用于根据求解到的X对光谱图像Y进行分类，

光谱图像Y中各像素Y_i,j的类别Class(Y_i,j)由最小类重构误差得到如下，

$Class\left(Y_{i,j}\right)=\arg \underset{m=1,...,M}{min}\left|\right|Y_{i,j}-a_m{X}_{m,i\×P+j}\left|\right|,i=1,2...K;j=1,2...P$

上式表示如果1到M中的某个值m令||Y_i,j-a_mX_m,i×P+j||取到最小值，则Y_i,j就归属到端元库中的第m个类别；X_m,i×P+j表示丰度系数矩阵X坐标(m,i×P+j)处数值。

各模块具体实现可参见相应步骤，本发明不予赘述。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。