一种基于空谱特征联合约束的高光谱异常检测方法
技术领域
本发明属于高光谱图像处理领域，特别涉及一种基于空谱特征联合约束的高光谱异常检测方法。
背景技术
在高光谱图像处理领域，未能获取目标光谱先验的目标检测称为异常检测。鉴于获取目标先验具有一定的难度和局限性，所以无需先验的无监督异常检测算法具有更为广泛的应用场景和更鲜明的实用价值。高光谱图像中的“异常”，是被定义为具有与绝大部分背景像元具有显著光谱特征差异的小部分其它像元。例如，由高光谱相机拍摄获得的大面积空旷的停机坪上，一两架单独的飞机可被定义为该幅高光谱图像中的异常目标，而其他部分例如跑道等则可被视为背景。
徐超，詹天明(《基于低秩全变差正则化的高光谱异常检测方法》，计算机科学与探索，2020，14(12)：2140—2149)通过对建立低秩分解模型，并在光谱维度施加全变分正则化约束增强异常检测结果。但该方法采用核范数近似矩阵秩，具有放缩误差，本发明中采用逼近程度更高的S1/2范数近似矩阵秩，提高了精度；该论文中仅对光谱维度进行了变分约束，忽视了对高光谱图像中空间域信息的利用，本发明同时对光谱维度和空间维度进行变分约束，有效利用了高光谱图像中的空间信息；同时，本发明在建模过程中，考虑了图像中的噪声并在求解过程中对噪声予以约束，进一步增强了异常检测的效果。
发明内容
本发明解决的技术问题是：为降低噪声对高光谱异常检测的影响，克服现有技术未考虑空间域信息的缺陷，同时本发明提出一种基于空谱特征联合约束的高光谱异常检测方法，在传统低秩模型的基础上，首先使用矩阵的S1/2范数替代核范数逼近矩阵秩，提高近似精度；为降低噪声对高光谱异常检测的影响，进一步添加去噪约束项，提高异常检测的精度，实现了对光谱信息的充分挖掘与高度利用；本发明进一步采用全变分正则化约束项，分别引入2个空间维度和1个光谱维度上的变分算子，实现图像平滑的同时提高了高光谱图像对空间信息的利用率。
本发明的技术方案是：一种基于空谱特征联合约束的高光谱异常检测方法，步骤1：建立基于S1/2范数与去噪的光谱信息低秩模型，包括以下子步骤：
步骤1.1：设任意一幅原始的高光谱图像表示为其中H和W分别表示高光谱图像中每一个波段的高度和宽度，d表示高光谱图像的波段数；
步骤1.2：将Ω按波段维度进行张量的矩阵化，得到二维矩阵其中d为波段数，n表示像元数目，数量上有n＝H×W；
步骤1.3：对矩阵化后的高光谱D进行建模得到
D＝Z+S+G
其中，表示背景矩阵，表示异常矩阵，表示噪声矩阵；
步骤1.4：建立S1/2低秩表示模型：
s.t.D＝Z+S+G
其中，最小化表征的是背景矩阵Z的低秩性，而表示的是异常矩阵的稀疏性；
步骤1.5：对步骤1.4中的S1/2低秩表示模型施加去噪约束项，
步骤1.6：建立优化目标函数：
步骤2：对背景矩阵Z进行变分约束，得到基于全变分正则化的空间信息约束模型：
SSTV(Z)＝||DhZDs||1+||DvZDs||1
其中，算子是空间域变分运算符，是光谱域变分运算符；矩阵的1范数是矩阵中元素的绝对值之和。
步骤3：根据步骤1得到的光谱信息低秩模型和步骤2得到的空间信息约束模型，建立最终高光谱异常检测优化目标函数如下所示：
其中，λ1,λ2,λ3是用于平衡各项重要性的系数；
步骤4：根据步骤3得到的高光谱异常检测优化目标函数，建立拉格朗日函数，基于交替方向乘子法，在优化其中一个变量时，固定其余所有变量进行求解，得到优化后的背景矩阵Z和异常矩阵S后，计算异常矩阵S每一行的l2范数，大于给定阈值时，可被判定为异常。
本发明进一步的技术方案是：针对空间域变分运算符，实质上是对每一个波段(Di,i＝1,2,...,d)上的H×W即n个像元进行变分操作。在沿光谱维度展开的二维矩阵数据模型下计算时，空间域变分运算符Dh和Dv均对二维高光谱矩阵的每一列进行差分，区别在于通过建立不同的运算符矩阵，Dh运算符差分的是每一个波段上水平相邻的像元，Dv运算符差分的是每一个波段上垂直相邻的像元。针对光谱域变分运算符，实质上是对同一个像元的不同波段进行变分操作，在沿光谱维度拉伸展开的二维矩阵数据模型下计算时，光谱域变分算子Ds是对二维高光谱矩阵的行进行差分。
本发明进一步的技术方案是：系数取值为[0.0001,0.001,0.01,0.1,1,10,100,1000]。
本发明进一步的技术方案是：采用ROC曲线和AUC值对本方法进行定性和定量分析。
发明效果
本发明的技术效果在于：采用本发明的有益效果主要包括：
(1)传统的低秩表示模型中一致采用矩阵的核范数，即矩阵奇异值之和近似矩阵的秩，而本发明提出以矩阵的S1/2范数近似矩阵的秩，即奇异值开根号之和，相对于核范数可以以更高的精度逼近矩阵的秩；
(2)传统的高光谱异常检测算法并未考虑噪声的影响，仅对背景和异常进行建模，本发明针对这一问题，在目标函数中对噪声建模，考虑了噪声对高光谱图像异常检测的影响，提高了异常检测的精度；
(3)已有的高光谱异常检测算法中绝大部分仅考虑了光谱特征，忽视了对高光谱图像空间特征的利用和对空间域信息的约束，本发明针对这一特性，对空间维度和光谱维度施加全变分正则化约束，实现了对高光谱图像空间特征的利用。
附图说明
图1是方法流程图：
图2分别是测试数据集的真值图和本发明方法以及其它几个对比方法的检测效果图对比示意图，其中(a)为真值图，(b)为RX方法图，(c)为LRX方法图，(d) 为CRD方法图，(e)为LRASR方法图，(f)为LSMAD方法图，(g)GTVLRR方法图， (h)本发明方法图
图3是本发明在测试数据集下与其它几个对比方法的检测结果ROC曲线图，该曲线是以阈值截断的方式绘制在不同阈值下模型的效果，上方曲线模型优于被其包住的下方曲线。
图4为集中在图像中间区域建筑群的异常目标
具体实施方式
在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”、“顺时针”、“逆时针”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
参加图1-图3，本发明的技术方案为：
(1)基于Schatten 1/2quasi-norm(S1/2范数)与去噪的光谱信息低秩模型
对任意一个m行n列的矩阵其秩函数、S1/2范数和核范数的定义分别如下所示：
其中，σi(i＝1,...,min(m,n))是矩阵的奇异值。
设任意一幅原始的高光谱图像表示为其中H和W分别表示高光谱图像中每一个波段的高度和宽度，d表示高光谱图像的波段数。本发明是在二阶矩阵空间上对高光谱图像进行建模处理，因此将原图像按波段维度进行张量的矩阵化，得到二维矩阵其中d为波段数，n表示像元数目，数量上n等于H和W的乘积，即n＝H×W。对拉伸后的高光谱数据矩阵D进行建模可得：
D＝Z+S+G
其中，表示背景矩阵，表示异常矩阵，表示噪声矩阵。低秩表示模型则可以分别采用S1/2范数和矩阵的L2,1范数对低秩性和稀疏性进行约束表示。故可以建立如下所示的S1/2低秩表示模型：
s.t.D＝Z+S+G
其中，最小化表征的是背景矩阵Z的低秩性，而即表示异常矩阵S每一行的l2范数之和，表征了异常矩阵S的稀疏性。上述低秩表示模型实现了对高光谱数据矩阵整体结构的刻画。
同时，为进一步增强异常检测的效果，对上述模型施加去噪约束项，即：
因此建立优化目标函数如下：
上述目标函数实现了基于光谱信息对高光谱的建模和约束。
(2)基于全变分正则化的空间信息约束模型
上述模型充分利用了高光谱的光谱信息，但忽略了对高光谱图像空间信息的利用。本发明针对背景数据平滑的先验对矩阵Z进行变分约束，分别左乘空间域变分运算符、右乘光谱域变分运算符，即：
SSTV(Z)＝||DhZDs||1+||DvZDs||1
其中，算子是空间域变分运算符，是光谱域变分运算符。具体说明如下：
针对空间域变分运算符，实质上是对每一个波段(Di,i＝1,2,...,d)上的H×W即n个像元进行变分操作。在沿光谱维度展开的二维矩阵数据模型下计算时，空间域变分运算符Dh和Dv均对二维高光谱矩阵的每一列进行差分，区别在于通过建立不同的运算符矩阵，Dh运算符差分的是每一个波段上水平相邻的像元，Dv运算符差分的是每一个波段上垂直相邻的像元。
针对光谱域变分运算符，实质上是对同一个像元的不同波段进行变分操作，在沿光谱维度拉伸展开的二维矩阵数据模型下计算时，光谱域变分算子Ds是对二维高光谱矩阵的行进行差分。
(3)基于空谱特征联合约束的高光谱异常检测模型
综合考虑上述空间域信息与光谱域信息，建立最终高光谱异常检测优化目标函数如下所示：
其中，λ1,λ2,λ3是用于平衡各项重要性的参数即系数，在实验中通过人为设定的方式确定，在一定范围内调整以获得最优效果。本发明中的调节范围为 [0.0001,0.001,0.01,0.1,1,10,100,1000]。
(4)模型优化与求解
根据目标函数和约束条件，建立拉格朗日函数，基于交替方向乘子法，在优化其中一个变量时，固定其余所有变量，具体求解过程参见下方具体实例描述。得到优化后的背景矩阵Z和异常矩阵S后，计算每一行的l2范数，大于给定阈值的可被判定为异常。
针对高光谱异常检测方法，本发明采用经典的受试者操作特性曲线(ReceiverOperating Characteristics,ROC曲线)和曲线下面积值(Area Under Curve,AUC)对检测性能进行定性和定量分析。ROC曲线是由一系列的目标检测率和虚警率点对绘制而成，可以直接反映两种检测概率之间的这种关系，每当给定一个判别阈值，可以准确计算出对应的一组检测率和虚警率数值，分别设置不同的阈值得到一系列的点对。为进行准确的定量分析，通过对ROC曲线下的面积进行积分，可以得到AUC值，衡量方法的优劣，直观地评估异常检测方法的性能。
下面以某幅沿海岸城市背景高光谱图像为例，结合附图对本发明的技术方案进行详细说明。在该数据中，背景是沿海岸的城市图像，被认定为异常目标的是集中在图像中间区域的建筑群，如下图所示。本发明算法可以在无光谱先验信息的条件下检测出该幅城市图像中的建筑群的位置，可用于进一步城市规划建设。 H＝W＝100,d＝207，即每一个波段上都是100×100大小的图像，共207个波段，经张量的矩阵化展开后，得到张量矩阵化后的结果为其中n＝10000,d＝207，即每个波段上包括10000个像元，共207个波段。最终模型如下：
即：
步骤1：交替方向乘子法求解模型。为求解目标函数，引入辅助变量P,Q,V1，模型变形为：
s.t.P＝DhZDs,Q＝DvZDs,Z＝V1
构建如下的增广拉格朗日函数：
其中，B1,B2,B3为引入的拉格朗日乘子。
根据优化原则，分别对函数中的变量进行如下的依次优化：
1)优化变量V1：相关文献证明得出，上述优化方程结果为：
V1 *＝Hλ(Z+B3)＝Udiag(Hλ(σ))VT
其中，Hλ(σ)＝(hλ(σ1),hλ(σ2),...,hλ(σr))T，U和V是对矩阵Z+B3进行奇异值分解所得的左奇异矩阵和右奇异矩阵，即：
Z+B3＝Udiag(σ)VT
其中，λ＝1/ν。
2)优化变量P：
3)优化变量Q：相关文献证明得出，上述形式的优化方程可以采用软阈值算子求解，即：
4)优化变量Z：
根据相关文献以及矩阵克罗内克积的性质：即若已知矩阵A＝BCD，有小写字母分别代表大写字母所对应的矩阵的向量形式，即矩阵列的垂直堆叠。在上式优化矩阵Z的过程中，可令将优化矩阵Z的函数重写为：
经求导后，可得：
其中，可用LSQR方法仿真求解。
5)优化变量S：相关文献证明得出，形如的优化结果为：
步骤2：最终得出稀疏异常矩阵S，计算每一行的l2范数，大于设定阈值η的对应列的所对应的像元即可认定为异常像元。根据不同的阈值绘制ROC曲线，可判定方法性能。
表1以量化的方式和具体数值的形式证明本发明模型的优越性。数值结果是AUC值，指代ROC曲线的面积，数值越大模型效果越好
表1 AUC值
