CN107609573B - 基于低秩分解和空谱约束的高光谱图像时变特征提取方法 - Google Patents

Abstract

本发明涉及一种基于低秩分解和空谱约束的高光谱图像时变特征提取方法，主要包括以下步骤：计算在不同时间同一地点采集的两组高光谱遥感图像的差异图像；利用差异图像的内在数据结构特性建立带空谱约束的低秩矩阵分解模型；通过交替迭代的方式求解该模型的各分量，进而提取时域变化特征。本发明充分利用数据的内在结构，提出一种带空谱约束的低秩矩阵分解模型及其求解算法，有效提取时域变化特征、去除多种形式的噪声，增强真实变化、抑制噪声所导致的虚假变化，从而提高时变高光谱遥感图像的变化检测精度。

Description

基于低秩分解和空谱约束的高光谱图像时变特征提取方法

技术领域

本发明涉及一种基于低秩分解和空谱约束的高光谱图像时变特征提取方法。

背景技术

遥感技术是本世纪六十年代发展起来的新兴综合技术，与空间、电子光学、计算机、地理学等科学技术紧密相关，是研究地球资源环境的有力技术手段。高光谱遥感是将成像技术与光谱技术相结合的多维信息获取技术。高光谱成像仪在电磁波谱的数十至数百个狭窄且连续的波段上同时探测目标的二维几何空间与一维光谱信息，为地物信息的提取和分析提供了极其丰富的数据，从而被广泛应用于地质科学、水文科学、精准农业以及军事领域等。

近年来，高光谱遥感技术被用于对同一地区进行长时间、多时段的观测，随之产生了时变高光谱图像。变化检测是识别多幅图像中变化的技术，已广泛用于灰度图像、彩色图像、多光谱遥感图像的处理，并且可用于时变高光谱图像。现有的时变高光谱图像变化检测方法主要由非时变高光谱图像的处理方法(如目标探测、像元解混)或其他图像(如多光谱图像)的变化检测方法改进而来。这些方法不乏效果优良者，对时变高光谱变化检测的研究做出了重大贡献，例如S.Liu，L.Bruzzone，F.Bovolo，and P.Du，“Hierarchicalunsupervised change detection in multitemporal hyperspectral images，”IEEETrans.Geosci.Remote Sens.，vol.53，no.1，pp.244-260，Jan.2015.提出了“变化端元”的概念，通过层次化的象元解混方法识别出多种主次变化。然而现有研究成果仍存在一定不足：其一，忽视了时变高光谱图像中普遍存在且形式多样化的噪声，包括配准误差、光谱偏差、硬件限制等因素所导致的异常点和热噪声，从而无法有效抑制噪声所产生“虚假变化”；其二，只考虑了高光谱图像的光谱信息，没有充分利用其空间特征或空谱特性。

发明内容

本发明的目的是提供一种带空谱约束的低秩矩阵分解模型及其优化求解算法，从而有效地提取时变高光谱图像的时域变化特征，抑制异常点和热噪声，提高变化检测精度。

为了达到上述目的，本发明的技术方案是提供了一种基于低秩分解和空谱约束的高光谱图像时变特征提取方法，其特征在于，包括以下步骤：

步骤1、计算在不同时间同一地点采集的两组高光谱遥感图像T₁∈R^M×B和T₂∈R^M×B的差异图像Y∈R^M×B，其中Y＝T₁-T₂；M为样本数，即像素个数，令I和J分别表示图像的行数和列数，则有M＝I×J；B为样本维度，即波段数；

步骤2、建立带空谱约束的低秩矩阵分解模型：

Y＝L+S+N (1)

式(1)中，L∈R^M×B为低秩数据，且具有空谱约束所刻画的空谱特性，用于表示时域变化特征(以下简称“时变特征”)；S∈R^M×B为稀疏矩阵，用于描述异常点，即分布稀疏、强度较高的噪声；N∈R^M×B为高斯噪声，一般表示服从高斯分布、强度较低的热噪声。引入变量X∈R^M×B，得到式(1)所述模型的优化求解函数：

式(2)中，τ为约束系数；ε为阈值；|·||_SS为本发明提出的空谱约束项：

式(3)中，x_ij∈R^B×1为ΘX∈R^I×J×B中当前被处理的像素，Θ是将J×B矩阵重新排列成I×J×B三阶张量的算子；

为以x_ij为中心的W×W空间邻域中除x_ij之外的像素；W和

均为预先设定的常数，调节W和

的值可以控制时变特征的平滑程度，防止过平滑；令W为正奇数，

在ΘX中的空间坐标为(i_n，j_n)，i，i_n∈{1，2，…，I}，j，j_n∈{1，2，…，J}，且(i_n，j_n)≠(i，j)，则有

步骤3、采取交替迭代的方式对式(1)和式(2)的L，S，X，N进行求解，包括以下步骤：

步骤3.1、通过增广拉格朗日乘子法将式(2)转化为：

式(4)中，Λ₁和Λ₂是拉格朗日乘子；μ是惩罚因子。由式(4)得到：

式(5)、(6)、(7)中，t表示迭代次数，且t＝0，1，…，t_max，t_max为迭代次数上限；

步骤3.2、整理式(5)、(6)、(7)、(1)依次得到L，S，X，N的解。

优选地，所述步骤3.2包括：

步骤3.2.1、结合式(5)和双边投影法，得：

式(8)中：

其中：

A₁∈R^B×r和A₂∈R^M×r为随机生成的双边投影矩阵；Q₁∈R^M×M和R₁∈R^M×r是Z₁的QR分解矩阵，Q₂∈R^{B ×B}和R₂∈R^B×r是Z₂的QR分解矩阵；q为双边投影的加速因子，r为L的秩的限定阈值，q和r均为是预先设定的常数；

步骤3.2.2、由式(6)及软阈值限定法，得：

式(9)中，

为软阈值限定算子，且

步骤3.2.3、对于式(7)，令

X＝[x₁，x₂，…，x_m，…，x_M]^T(其中m＝1，2，…，M，m＝(j-1)×I+i，i＝1，2，…，I，j＝1，2，…，J)，且令J(X)对X的偏导为零，得X^(t+1)中各向量

的闭式解如下：

式(10)中，

为全1向量；

表示X^(t)中以

为中心的W×W空间邻域内的像素；

对于

和x_ij的空间位置相同，

和

的空间位置相同，故式(10)中空间邻域的定义、W和ω_w的定义以及w与

的空间坐标的对应关系均同式(3)；

步骤3.2.4、由式(1)及式(8)至式(10)，得：

N^(t+1)＝Y-L^(t+1)-S^(t+1) (11)

步骤3.2.5、根据式(8)至式(11)，交替更新L，X，S，N，迭代求解其最优解；当迭代误差Error_i达到既定阈值ε_i，ε_i＞0，i＝1，2，或迭代次数t达到t_max时，迭代终止，所得L或X即为时域变化特征。

优选地，步骤3.2.5中，迭代误差Error_i定义如下：

Error₂＝||L^(t+1)-X^(t+1)||_∞ (13)

综上，将本发明所实现的算法记为LRSD_SS，考虑到式(4)中拉格朗日乘子以及式(2)到式(13)中各系数、因子的迭代更新，其具体步骤如下：

a)输入T₁∈R^M×B和T₂∈R^M×B，τ，λ，r，q，ε₁，ε₂，t_max，μ的初始值μ₀，μ的更新步长ρ，以及μ的阈值μ_max；

b)计算Y＝T₁-T₂，并初始化：令迭代次数t＝0，L⁽⁰⁾＝Y，X⁽⁰⁾＝0，S⁽⁰⁾＝0，Λ₁ ⁽⁰⁾＝0，Λ₂ ⁽⁰⁾＝0，μ＝μ₀；

c)由式(8)、(9)、(10)、(11)分别计算L^(t+1)，X^(t+1)，S^(t+1)，N^(t+1)；

d)令

e)令μ←min(ρμ，μ_max)；

f)由式(12)、(13)分别计算Error₁，Error₂；

g)令t←t+1；

h)若Error₁＞ε1，Error₂＞ε₂，或t＜t_max，则返回c)；否则转至i)；

i)令L＝L^(t+1)，S＝S^(t+1)，N＝N^(t+1)，X＝X^(t+1)，并输出。

本发明充分利用高光谱图像的时变特征与噪声在固有数据结构上的差异以及时变特征的局部空谱特性，提出了一种带空谱约束的低秩矩阵分解模型及其优化求解算法，从而有效地提取时域变化特征，抑制异常点和热噪声，提高变化检测精度。在设计空谱约束项和求解算法时，本发明还考虑到计算复杂度和收敛性等问题，以提高算法LRSD_SS的实用性。

本发明具有如下特点：

1)本发明充分利用时变高光谱图像的内在数据结构特性来设计形如式(2)的时变特征表达模型：分别对低秩分解模型中的低秩分量、稀疏分量、高斯分量和差异图像中的时变特征、异常点、热噪声建立联系，并且通过形如式(3)的空谱约束来利用时变特征的局部空谱特性，同步实现时变特征提取和鲁棒性降噪。

2)所述空谱约束反应的是邻域的空谱距离，而非某波段像素的空间距离或某两个像素间的光谱距离，因此能够有力刻画时变特征的局部空谱特性；

3)所述空谱约束项包含邻域尺寸W和权值

这两个参数，可通过预先调节W和

的值来控制时变特征的局部平滑程度，防止过平滑。

4)所述空谱约束项采用l2范数来刻画邻域的空谱距离，使得X具有闭式解；该约束项还能同步利用高光谱图像的空间和光谱信息，而非逐波段进行空间平滑，节约了本发明所提出算法的运行时间；

5)在求解低秩分量L时采用了双边投影法，其计算复杂度为r²(M+3B+4r)+(4q+4)MBr；

6)得益于X的闭式解以及拉格朗日乘子法和双边投影法自身良好的收敛性，本发明的求解算法LRSD_SS具备良好的收敛性。

由于采用了上述技术方案，本发明和现有技术相比，具有以下的优点和积极效果：

1)针对时变高光谱图像的内在数据结构进行分析和建模，并以此来设计时变特征提取方法：因为(a)高光谱图像为高维数据，而拍摄区域内地表变化类别数量有限，每种变化对应的像素一般具有较强的类内相关性，故时变特征一般存在于某个低维数据空间中，具有低秩性；(b)时变高光谱图像的噪声由配准误差、光谱偏差、硬件限制等多重因素导致，主要分为异常点和热噪声，前者具有稀疏性且强度较高，后者一般遵循高斯分布且强度较低，噪声可表现为“虚假变化”，干扰真实变化的检测，所以建立式(2)所述的表达模型，采用能够同时分离低秩数据、稀疏数据和热噪声的低秩矩阵分解方法进行特征提取和鲁棒性降噪；

2)本发明通过式(3)所述的空谱约束项进一步挖掘数据的内在结构特性，例如某些特殊噪声具有类似于时变特征的内在数据结构(例如坏死点兼具低秩性和稀疏性)，而时变特征则具有一般噪声不具备的局部平滑性，从而提升时变特征提取效果；特别地，该空谱约束项反应的是邻域的空谱距离，而非某波段像素的空间距离或某两个像素间的光谱距离，因此能够有力刻画时变特征的局部空谱特性；该空谱约束还允许以预设并调节参数的方式来防止时变特征的过平滑；

3)本发明不仅考虑到时变特征提取的有效性，还兼顾了求解算法的实用性：(a)空谱约束项通过l₂范数来刻画邻域空谱距离，使得X具有闭式解，并且能够同步处理空间和光谱信息，而非逐波段进行空间平滑，保障了算法的效率；(b)采用收敛性良好的拉格朗日乘子法和双边投影法以及具有闭式解的空谱约束，从而使求解算法具有良好的收敛性。

附图说明

图1为本发明的流程图；

图2(a)为仿真时变高光谱图像的时域变化地物真实图，其中白色代表变化区域，黑色代表非变化区域；

图2(b)和图2(c)分别为仿真时变高光谱图像中T₁∈R^M×B和T₂∈R^M×B的伪彩色图，其R、G、B通道分别为相应高光谱图像的第10、20、40个波段；

图3为时变特征的幅度值分布的箱线图，其中每个“箱”的上、中、下边缘分别对应第三个四分位数、第二个四分位数(中位数)、第一个四分位数，从上、下边缘延伸的“须”表示适度离群值；

图4(a)及图4(b)为本发明的求解算法LRSD_SS的收敛曲线；

图5(a)为真实时变高光谱图像的时域变化地物真实图，其中白色代表变化区域，黑色代表非变化区域；

图5(b)和图5(c)分别为真实时变高光谱图像中T₁和T₂的伪彩色图，其R、G、B通道分别为相应高光谱图像的第20、40、60个波段；

具体实施方式

下面结合具体实施例，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围，此外应理解，在阅读本发明讲授的内容之后，本领域技术人员可以对本发明做各种改动和修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

本发明的实施方式涉及一种基于低秩分解和空谱约束的高光谱图像时变特征提取方法，该算法包括以下步骤：计算在不同时间同一地点采集的时变高光谱遥感图像的差异图像；根据差异图像的内在数据结构特性建立带空谱约束的低秩矩阵分解模型；通过交替迭代求解该模型的各分量，进而提取时域变化特征。将所得时变特征L的幅度向量A∈R^M×1(其中A(m)＝||L(m，：)||_F，m＝1，2，…，M)代入K-means或其他能够实现无监督二分类的分类器中，可以将A(m)标记为“变化类”或“非变化类”，从而在图像空间域上划分出变化区域和非变化区域，实现时变高光谱图像的变化检测。

仿真数据实验

仿真时变高光谱图像由真实的非时变高光谱图像“Pavia大学数据”构造而来。该图像由http：//www.ehu.es/ccwintco/uploads/e/ee/PaviaU.mat提供，其采集时间为2001年，地点为意大利的Pavia大学，成像仪为ROSIS，具有610×340个像素(即I＝610，J＝340)，去掉噪声污染、水吸收波段后，余下103个波段(即B＝103)。仿真时变高光谱图像的构造方法如下：1)将原Pavia大学数据作为在时刻1采样的高光谱图像T₁₀∈R^M×B(M＝I×J＝207400)；2)在T₁₀中选择6个区域，将两两区域内的像素互换，模拟地物的时域变化，得到时刻2采样的高光谱图像T₂₀∈R^M×B以及变化区域的地物真实图(如图2(a)所示)；3)分别对T₁₀和T₂₀添加如表1所述的异常点和热噪声，得到仿真时变高光谱图像T₁∈R^M×B和T₂∈R^M×B，其伪彩色图分别如图2(b)和图2(c)所示。

表1仿真热噪声和异常点(均为归一化的随机噪声)

采用K-means作为时变特征提取后的分类器，实现变化检测。本发明算法LRSD_SS的参数设置如下：τ＝0.01，μ_max＝10⁶，ρ＝1.05，

q＝3，W＝3，

ε₁＝ε₂＝10^-6，t_max＝30，其中i＝1，2，…I，j＝1，2，…J。所有算法运行的硬件环境为Intel(R)Xeon(R)X5667CPU 3.00GHz(双核)24GB内存，软件平台为Windows 7及MATLAB R2013b。

实验1算法效果的验证

算法效果的评价主要包括以下三种方式：(a)考察时变特征的幅度值分布特性；(b)比较分类标识图与地物真实图；(c)计算总体分类精度(Overall Accuracy，OA)，平均分类精度(Average Accuracy，AA)以及Kappa系数(κ)。

分别对图2(a)所示的变化类(Change)和非变化类(Non-Change)，考察无噪差异图像Y₀(Y₀＝T₁₀-T₂₀)、差异图像Y(Y＝T₁-T₂)、本发明算法LRSD_SS所得时变特征L、由W.He，H.Zhang，L.Zhang，and H.Shen，“Total-variation-regularized low-rank matrixfactorization for hyperspectral image restoration，”IEEE Trans.Geosci.RemoteSens.，vol.54，no.1，pp.176-188，Jan.2016.提出的带全变分空间约束的低秩矩阵分解算法LRSD_TV所得L以及由T.Zhou and D.Tao，“GoDec：randomized low-rank&sparse matrixdecomposition in noisy case，”Proc.28th ICML，Jan.2011，pp.33-40.提出的经典的低秩矩阵分解算法LRSD所得L的幅度值分布特性，如图3所示，本发明LRSD_SS的时变特征幅值具有较为紧密的类内分布和较为分散的类间分布，易于被分类器正确识别。如表2所示，当分类器取为经典的无监督聚类算法K-means时，本发明提出的LRSD_SS算法能够得到较高的分类精度，其变化检测效果优于表中经典的(SAM和PCA)或先进的时变特征提取算法(ASCD，LRSD和LRSD_TV)。

表2时变特征提取算法+K-means的变化检测结果(仿真时变高光谱图像)

实验2算法效率的验证

令LRSD_SS和LRSD_TV均采用双边投影法求解低秩分量。如表3所示，在求解L和S上，这两种带约束的低秩分解算法所耗时间相近；而在求解X上，因为本发明LRSD_SS采用的空谱约束使得X具有闭式解，所以该算法的耗时远少于LRSD_TV(其空间约束需要逐波段计算，且无闭式解，故耗时较多)。由此可知，本发明的计算复杂度不高，运行时间可以接受，其效率高于先进的低秩矩阵分解算法LRSD_TV。

表3单次迭代耗时(单位：秒；实验数据：仿真时变高光谱图像)

实验3算法收敛性的验证

如图4(a)及图4(b)所示，本发明提出的算法LRSD_SS经过20次左右的迭代即可收敛，具有一定的实用价值。

真实数据实验

本实验所用真实时变高光谱图像由c.Wu，B.Du，and L.Zhang，“A subspace-basedchange detection method for hyperspectral images，”IEEEJ.Sel.Topics.Appl.Earth Observ.Remote Sens.，vol.6，no.2，pp.815-830，Apr.2013.提供，其成像仪为Hyperion，T₁∈R^M×B和T₂∈R^M×B(其中M＝I×J＝63000，I＝450，J＝140，B＝155)的采集时间分别为2006年5月3日和2007年4月23日，采集地点为中国江苏盐城某农业用地。变化区域的地物真实如图5(a)所示，T₁和T₂的伪彩色图分别如图5(b)和图5(c)所示。除了实验数据，本实验其余设置同仿真数据实验。

如表4所示，对于图5所示的真实时变高光谱图像，当分类器取为经典的无监督聚类算法K-means时，本发明提出的LRSD_SS算法能够得到较高的分类精度，其变化检测效果优于表中经典的(SAM和PCA)或先进的时变特征提取算法(ASCD，LRSD和LRSD_TV)。

表4时变特征提取算法+K-means的变化检测结果(真实时变高光谱图像)

Claims (3)

1.一种基于低秩分解和空谱约束的高光谱图像时变特征提取方法，其特征在于，包括以下步骤：

步骤1、计算在不同时间同一地点采集的两组高光谱遥感图像T₁∈R^M×B和T₂∈R^M×B的差异图像Y∈R^M×B，其中Y＝T₁-T₂；M为样本数，令I和J分别表示图像的行数和列数，则有M＝I×J；B为样本维度；

步骤2、建立带空谱约束的低秩矩阵分解模型：

Y＝L+S+N (1)

式(1)中，L∈R^M×B为低秩数据，且具有空谱约束所刻画的空谱特性，用于表示时域变化特征；S∈R^M×B为稀疏矩阵，用于描述异常点；N∈R^M×B为高斯噪声；引入变量X∈R^M×B，得到式(1)所述模型的优化求解函数：

式(2)中，τ为约束系数；

ε为阈值；||·||_SS为空谱约束项：

式(3)中，x_ij∈R^B×1为ΘX∈R^I×J×B中当前被处理的像素，Θ是将J×B矩阵重新排列成I×J×B三阶张量的算子；

为以x_ij为中心的W×W空间邻域中除x_ij之外的像素；W和

均为预先设定的常数，调节W和

的值可以控制时变特征的平滑程度，防止过平滑；令W为正奇数，

在ΘX中的空间坐标为(i_n,j_n)，i,i_n∈{1,2,…,I}，j,j_n∈{1,2,…,J}，且(i_n,j_n)≠(i,j)，则有

步骤3、采取交替迭代的方式对式(1)和式(2)的L，S，X，N进行求解，包括以下步骤：

步骤3.1、通过增广拉格朗日乘子法将式(2)转化为：

式(4)中，Λ₁和Λ₂是拉格朗日乘子；μ是惩罚因子，由式(4)得到：

式(5)、(6)、(7)中，t表示迭代次数，且t＝0,1,…,t_max，t_max为迭代次数上限；

步骤3.2、整理式(5)、(6)、(7)、(1)依次得到L，S，X，N的解，包括：

步骤3.2.1、结合式(5)和双边投影法，得：

式(8)中：

其中：

A₁∈R^B×r和A₂∈R^M×r为随机生成的双边投影矩阵；Q₁∈R^M×M和R₁∈R^M×r是Z₁的QR分解矩阵，Q₂∈R^B×B和R₂∈R^B×r是Z₂的QR分解矩阵；q为双边投影的加速因子，r为L的秩的限定阈值，q和r均为是预先设定的常数；

步骤3.2.2、由式(6)及软阈值限定法，得：

式(9)中，

为软阈值限定算子，且

步骤3.2.3、对于式(7)，令

X＝[x₁,x₂,…,x_m,…,x_M]^T，其中m＝1,2,…,M，m＝(j-1)×I+i，i＝1,2,…,I，j＝1,2,…,J，且令J(X)对X的偏导为零，得X^(t+1)中各向量

的闭式解如下：

式(10)中，

为全1向量；

表示X^(t)中以

为中心的W×W空间邻域内的像素；

对于

和x_ij的空间位置相同，

和

的空间位置相同，故式(10)中空间邻域的定义、W和ω_w的定义以及w与

的空间坐标的对应关系均同式(3)；

步骤3.2.4、由式(1)及式(8)至式(10)，得：

N^(t+1)＝Y-L^(t+1)-S^(t+1) (11)

步骤3.2.5、根据式(8)至式(11)，交替更新L，X，S，N，迭代求解其最优解；当迭代误差Error_i达到既定阈值ε_i，ε_i>0,i＝1,2，或迭代次数t达到t_max时，迭代终止，所得L或X即为时域变化特征，迭代误差Error_i定义如下：

Error₂＝||L^(t+1)-X^(t+1)||_∞ (13)。

2.如权利要求1所述的一种基于低秩分解和空谱约束的高光谱图像时变特征提取方法，其特征在于，通过预先调节并设置W和

的值能够控制时域变化特征的局部平滑程度。

3.如权利要求1所述的一种基于低秩分解和空谱约束的高光谱图像时变特征提取方法，其特征在于，所实现的算法LRSD_SS包含式(4)中拉格朗日乘子，以及式(2)到式(13)中所有系数、因子的迭代更新，具体步骤如下：

a)输入T₁∈R^M×B和T₂∈R^M×B，τ，λ，r，q，ε₁，ε₂，t_max，μ的初始值μ₀，μ的更新步长ρ，以及μ的阈值μ_max；

b)计算Y＝T₁-T₂，并初始化：令迭代次数t＝0，L⁽⁰⁾＝Y，X⁽⁰⁾＝0，S⁽⁰⁾＝0，Λ₁ ⁽⁰⁾＝0，Λ₂ ⁽⁰⁾＝0，μ＝μ₀；

c)由式(8)、(9)、(10)、(11)分别计算L^(t+1)，X^(t+1)，S^(t+1)，N^(t+1)；

d)令

e)令μ←min(ρμ,μ_max)；

f)由式(12)、(13)分别计算Error₁，Error₂；

g)令t←t+1；

h)若Error₁>ε₁，Error₂>ε₂，或t<t_max，则返回c)；否则转至i)；

i)令L＝L^(t+1)，S＝S^(t+1)，N＝N^(t+1)，X＝X^(t+1)，并输出。

Images