[发明专利]基于背景字典学习和结构稀疏表示的高光谱异常目标检测方法

摘要

本发明公开了一种基于背景字典学习和结构稀疏表示的高光谱异常目标检测方法，用于解决现有高光谱异常目标检测方法目标检测效率低的技术问题。技术方案是在基于局部RX算法选择初始背景像元后，利用主成分分析字典学习法学习得到鲁棒性的背景字典。在稀疏向量求解和图像重建过程中，引入重加权拉普拉斯先验，提高稀疏向量求解精度。最后，根据原始图像与重建图像之间的误差实现异常目标的精确提取。在真实的高光谱卫星图像AVIRIS和仿真的高光谱数据集上的试验结果表明，本发明获得的检测结果相对于背景技术在恒虚警率的前提下检测率提高了8％～15％。

主权项

一种基于背景字典学习和结构稀疏表示的高光谱异常目标检测方法，其特征在于包括以下步骤：步骤一、基于主成分分析的鲁棒性背景字典学习；(1)采用双窗局部RX算法获取背景像素集合；根据输入高光谱图像分辨率选择窗口尺寸都为奇数的矩形外窗Outer和内窗Inner，Outer＝n_outer×n_outer,Inner＝n_inner×n_inner,n_outer>n_innner,n_outer、n_inner分别表示外窗和内窗的尺寸，都为奇数；以每一个输入像素x_i∈X,i＝1,…,n_p为中心，根据外窗和内窗尺寸获取局部背景区域n_local＝Outer‑Inner，n_local为n_b×N_local大小的局部背景矩阵，然后计算n_local的均值和协方差矩阵分别得到μ_local＝[μ₁,…,μ_N]^T和∑_local；局部RX计算公式如下：D_LocalRX(x_i)＝(x_i‑μ_local(i))^T(∑_local(i))^‑1(x_i‑μ_local(i))    (1)式中，μ_local(i)表示第i个输入像元所对应的局部背景的均值，(∑_local(i))^‑1表示第i输入像元所对应的局部背景的协方差矩阵的逆矩阵，利用公式(1)遍历所有高光谱像元后选取阈值进行背景分割，小于阈值的对应像元为背景像元，得到背景像元矩阵Z_local，大小为n_b×N_Z，N_Z表示局部背景像元的个数；(2)基于主成分分析的鲁棒性背景字典学习；将背景像元矩阵Z_local转置后求其协方差矩阵，定义为CovZ_local，大小为n_b×n_b；然后求协方差矩阵CovZ_local的特征值V和特征向量，将特征值按照从大到小排列后，相应的特征向量也依次排列，得到的新的特征向量矩阵即为最终学习到的背景字典D；步骤二、建立基于重加权拉普拉斯稀疏先验的结构稀疏表示模型；在稀疏表示框架下，高光谱图像像元能够被学习到的字典D稀疏表示，即X＝DY+N    (2)式中，D为已学习到的背景字典，Y＝[y₁,…,y_n]为稀疏表示矩阵，N为稀疏表示误差和图像噪声；假设N服从的矩阵正态分布，则对应的似然函数为$p\left(X\right|Y,\mathrm{\&lambda;})=\frac{\exp \{-\frac{1}{2}||DY-X|{|}_{{\mathrm{\&Sigma;}}_n}^2\}}{{\left(2\mathrm{\&pi;}\right)}^{n_b{n}_p/2}{\left|{\mathrm{\&Sigma;}}_n\right|}^{n_p/2}},{\mathrm{\&Sigma;}}_n=diag\left(\mathrm{\&lambda;}\right),---\left(3\right)$其中，Σ_n＝diag(λ)是以λ元素为对角线元素的对角矩阵，用于表示误差和噪声的强度；表示Q矩阵的加权迹范数；为了表示稀疏向量内部的结构稀疏性，将重加权拉普拉斯稀疏先验引入Y；首先假设Y服从如下分布$p\left(Y\right|\mathrm{\&gamma;})=\frac{\exp \{-\frac{1}{2}|\left|Y\right|{|}_{{\mathrm{\&Sigma;}}_y}^2\}}{{\left(2\mathrm{\&pi;}\right)}^{n_b{n}_p/2}|{\mathrm{\&Sigma;}}_y{|}^{n_p/2}},{\mathrm{\&Sigma;}}_y=diag\left(\mathrm{\&gamma;}\right),---\left(4\right)$其中，Σ_y＝diag(γ)表示以γ的元素为对角线元素的对角矩阵，控制Y中每一行的稀疏度，γ_i＝0表示Y的第i行为0；假设则其中任一列y_i服从高斯分布；假设超参数γ服从以下的伽马分布，$p\left(\mathrm{\&gamma;}\right|\mathrm{\&kappa;})=\underset{i=1}{\overset{n_b}{\&Pi;}}Gamma(1,\frac{2}{{\mathrm{\&kappa;}}_i})=\underset{i=1}{\overset{n_b}{\&Pi;}}\frac{{\mathrm{\&kappa;}}_i}{2}\exp (-\frac{{\mathrm{\&kappa;}}_i{\mathrm{\&gamma;}}_i}{2})---\left(5\right)$以上两级先验，等价于重加权拉普拉斯分布，因为对于y_i有$p\left(y_i\right)=\&Integral;p\left(y_i\right|\mathrm{\&gamma;})p\left(\mathrm{\&gamma;}\right|\mathrm{\&kappa;})d\mathrm{\&gamma;}=\frac{1}{2}\underset{j=1}{\overset{n_b}{\&Pi;}}{\mathrm{\&kappa;}}_j^{1/2}\exp (-|\left|{\mathrm{Ky}}_i\right|{|}_1)---\left(6\right)$其中，为保证求解便利，采用级联先验，而且λ，γ和κ均为待估计参数；由于λ，γ和κ未知，根据经验贝叶斯框架，先基于输入的高光谱数据X利用MAP估计未知参数λ，γ和κ，如下$\{{\mathrm{\&lambda;}}_{opt},{\mathrm{\&gamma;}}_{opt},{\mathrm{\&kappa;}}_{opt}\}=\arg \underset{\mathrm{\&lambda;},\mathrm{\&gamma;},\mathrm{\&kappa;}}{max}p(\mathrm{\&lambda;},\mathrm{\&gamma;},\mathrm{\&kappa;}|X)\&Proportional;\arg \underset{\mathrm{\&lambda;},\mathrm{\&gamma;},\mathrm{\&kappa;}}{max}\&Integral;p\left(X\right|Y,\mathrm{\&lambda;})p\left(Y\right|\mathrm{\&gamma;})p\left(\mathrm{\&gamma;}\right|\mathrm{\&kappa;})dY---\left(7\right)$其中，λ_opt,γ_opt,κ_opt分别表示最优的λ，γ，κ；通过积分，并引入‑2log运算，容易得知式子(7)等价于最小化如下的式子其中，tr(·)表示迹范数，Σ_by＝Σ_n+DΣ_yD^T，为代价函数；通过变形式子(8)，得到稀疏信号Y的非分离稀疏约束模型；首先，对式子(8)的第一份部分进行变形$tr\left(n_{n_p}^{-1}{X}^T{\mathrm{\&Sigma;}}_{by}^{-1}X\right)=tr\&lsqb;\left(\frac{X^T}{\sqrt{n_p}}\right){\mathrm{\&Sigma;}}_{by}^{-1}\left(\frac{X}{\sqrt{n_p}}\right)\&rsqb;=\underset{Y}{min}\left|\right|DY-\frac{X}{\sqrt{n_p}}|{|}_{{\mathrm{\&Sigma;}}_n}^2+|\left|Y\right|{|}_{{\mathrm{\&Sigma;}}_y}^2.---\left(9\right)$然后，将式子(9)带入到式子(8)中，得接着，引入新的代价方程如下显然，而且能够证明，最小化式子(8)再对稀疏信号Y进行MAP估计，与直接最小化式子(11)得到的λ，γ和κ相同，关于Y的解仅相差一个常量因此，式子(11)看作是关于稀疏信号Y的正则化回归模型，其中为稀疏信号的非分离稀疏约束；该约束不能拆分成对于Y中每一行的独立约束，因此该约束能同时约束稀疏信号中非零元素，潜在地考量非零元素之间的相关性；此外，Σ_by中包含了表征噪声强度的λ，因此，得到的稀疏约束随着估计的噪声强度自适应的变化，具有噪声鲁棒性；最终得到了如下的基于重加权拉普拉斯稀疏先验的结构稀疏表示优化模型：$\underset{Y,\mathrm{\&lambda;}\&GreaterEqual;0,\mathrm{\&gamma;}\&GreaterEqual;0,\mathrm{\&kappa;}}{\min }\left(\left|\right|DY-\frac{X}{\sqrt{n_p}}|{|}_{{\&Sigma;}_n}^2+|\left|Y\right|{|}_{{\&Sigma;}_y}^2+\log \left|{\&Sigma;}_{by}\right|+\underset{i=1}{\overset{n_b}{\&Sigma;}}\frac{({\mathrm{\&kappa;}}_i{\mathrm{\&gamma;}}_i-2{\mathrm{log\&kappa;}}_i)}{n_p}\right)---\left(12\right)$步骤三、模型求解并重建高光谱图像；已知待检测高光谱数据X，采用坐标下降法求解式子(12)，每次迭代中仅优化一个变量而固定剩余的变量；λ⁰，γ⁰，κ⁰分别表示初始的数值，t记录迭代次数，η表示更新阈值，λ^t，γ^t，κ^t，Y^t分别表示迭代到第t次时的数值，λ^t+1，γ^t+1，κ^t+1，Y^t+1分别表示迭代到第t+1次时的数值；具体步骤如下：①初始化，λ⁰，γ⁰，κ⁰均初始化为对应长度的全1向量，计数变量t＝0；②更新中间变量Σ_n＝diag(λ^t)，Σ_y＝diag(γ^t)，Σ_by＝Σ_n+DΣ_yD^T；③固定λ^t，γ^t和κ^t，根据式子(11)得到关于Y的优化形式，如下求解得到Y的更新规则如下，$Y^{t+1}={\mathrm{\&Sigma;}}_y{D}^T{\mathrm{\&Sigma;}}_{by}^{-1}\frac{X}{\sqrt{n_p}}---\left(14\right)$④固定Y^t+1，λ^t和κ^t，得到关于γ的优化形式，如下求解得到如下的更新形式：${\mathrm{\&gamma;}}_i^{t+1}=\frac{\sqrt{4{\mathrm{\&kappa;}}_i{n}_p{z}_i+n_p^2}-n_p}{2{\mathrm{\&kappa;}}_i}---\left(16\right)$其中，为γ^t+1的第i个元素，代表V^T+Y^t+1(Y^t+1)^T的对角线元素组成的向量，z_i为z的第i个元素；⑤固定Y^t+1，γ^t+1和κ^t，得到关于λ的优化形式，如下求解得到如下的更新形式：${\mathrm{\&lambda;}}^{t+1}=\sqrt{diag\left({\mathrm{QQ}}^T\right)./\mathrm{\&alpha;}}---\left(18\right)$其中，根号运算表示向量每一个元素开方后组成的向量，./运算代表两个向量对应元素相除后组成的向量，代表对角线元素组成的向量；⑥固定Y^t+1，γ^t+1和λ^t+1，得到关于κ的优化形式，如下:求解得到如下的更新形式：${\mathrm{\&kappa;}}^{t+1}=\frac{2}{\mathrm{\&gamma;}+d}---\left(20\right)$上式中的加法和除法运算均作用在向量的每一个元素上，得到一个新向量，d＝10‑⁶的引入是为了确保γ中出现0时，式子(20)依然有意义；⑦计算稀疏信号Y更新前后的差异，如下$\mathrm{\&eta;}=\frac{\left|\right|Y^{t+1}.*\sqrt{n_p}|{|}_F}{\left|\right|Y^t.*\sqrt{n_p}|{|}_F}---\left(21\right)$其中，表示对Y^t+1内的每一个元素乘以||·||_F表示弗罗贝尼乌斯范数(Frobenius norm),如果计数器t>5或者更新差异η<10^‑3，则退出循环；否则t+1，循环执行步骤②至⑦；⑧假设上述循环结束得到的最优稀疏信号为Y_rec，则待重建的高光谱图像X_rec通过如下方式得到：$X_{rec}=D(Y_{rec}.*\sqrt{n_p})---\left(22\right)$步骤四、基于重建误差提取异常目标；根据原图像与重建图像之间的重建误差提取异常目标，重建误差计算公式如下：r(X)＝||X‑X_rec||₂    (23)其中，X、X_rec分别表示原图像和重建后图像，都为n_b×n_p二维矩阵，||·||₂表示向量的2范数；得到的误差项r(X)为1行、n_p列大小的行向量，先将其按照最小值和最大值归一化到[0,1]区间，再重新排列为n_row行、n_col列的灰度图像；根据灰度图选取分割阈值δ，图像中大于δ的像素标记为1表示目标，小于等于δ的像素标记0表示背景；最终得到只有0和1的二值结果图，完成异常目标检测。