CN104036509B

CN104036509B - Method for unmixing hyperspectral mixed pixel based on compressed sensing

Info

Publication number: CN104036509B
Application number: CN201410276372.7A
Authority: CN
Inventors: 付宁; 徐红伟; 殷聪如; 乔立岩
Original assignee: Harbin Institute of Technology Shenzhen
Current assignee: Harbin Institute of Technology Shenzhen
Priority date: 2014-06-19
Filing date: 2014-06-19
Publication date: 2017-01-11
Anticipated expiration: 2034-06-19
Also published as: CN104036509A

Abstract

The invention relates to a method for unmixing a hyperspectral mixed pixel based on compressed sensing, and relates to the compressed sensing field and the hyperspectral remote sensing field. The method is used for solving the problem of low mixed pixel unmixing speed during hyperspectral image data acquisition by use of the traditional Nyquist sampling theorem. The method comprises the steps of firstly, inputting a measurement matrix Phi and a compressed measurement matrix Y and establishing a spectral mixing model Y=PhiXT=Phi(AS)T by use of the compressed sensing theory, secondly, performing iterative processing on the estimated value S^ of an end member abundance matrix S and the estimated value A^ of an end member spectrum matrix A, and if the difference of absolute values of every corresponding element in the estimated values A^ of the end member spectrum matrix A obtained by two adjacent times of iterative processing is smaller than 0.1, stopping iteration and outputting the end member abundance matrix S^ and completing the unmixing of the hyperspectral mixed pixel, otherwise, continuing the iterative processing. The method is mainly applied to unmixing the hyperspectral mixed pixels.

Description

A Method of Hyperspectral Mixed Pixel Decomposition Based on Compressed Sensing

技术领域technical field

本发明涉及压缩感知领域与高光谱遥感领域。The invention relates to the fields of compressed sensing and hyperspectral remote sensing.

背景技术Background technique

高光谱成像技术是遥感领域发展起来的一种新型的对地观测技术，典型的硬件设备为成像光谱仪。成像光谱仪通过分光技术，将电磁波信号分解为许多微小的、相邻的波段，对应波段上的能量被不同的传感器接收。因此与传统的遥感成像技术相比，高光谱成像技术具有图谱合一、光谱波段多、光谱分辨率高等特点，且在地表物质识别和分类、特征提取等方面有较大的优势。高光谱图像的空间分辨率较低，因此不可避免的会产生混合像元，混合像元的存在成为了高光谱遥感应用的进一步发展的障碍。如果对高光谱图像的混合像元分解得到各端元所对应的丰度，即对高光谱遥感数据进行“盲”分解将有非常重要的意义。Hyperspectral imaging technology is a new type of earth observation technology developed in the field of remote sensing. The typical hardware device is an imaging spectrometer. The imaging spectrometer decomposes the electromagnetic wave signal into many tiny, adjacent bands through spectroscopic technology, and the energy on the corresponding bands is received by different sensors. Therefore, compared with traditional remote sensing imaging technology, hyperspectral imaging technology has the characteristics of integrated map, multiple spectral bands, and high spectral resolution, and has great advantages in surface material identification, classification, and feature extraction. The spatial resolution of hyperspectral images is low, so it is inevitable to produce mixed pixels, and the existence of mixed pixels has become an obstacle to the further development of hyperspectral remote sensing applications. If the mixed pixels of the hyperspectral image are decomposed to obtain the corresponding abundance of each end member, that is, the "blind" decomposition of the hyperspectral remote sensing data will be of great significance.

盲源分离技术是指在源信号和混合方式均未知的情况下，根据源信号的统计特性，仅由观测到的混合信号恢复出源信号的过程。“盲”有两种解释：一种是除观测信号外，其他的信息均未知，称为“全盲”；另一种是信号或者混合系统的某些先验信息已知，称为“半盲”。针对混合像元分解的问题，盲源分离技术同时获得端元光谱和端元丰度，即实现混合像元的盲分解。在处理混合像元分解问题之前，需要建立光谱混合模型。光谱混合方式主要包括线性混合和非线性混合，线性光谱混合模型是目前应用最广泛的混合像元分解模型，其表达式为：Blind source separation technology refers to the process of recovering the source signal only from the observed mixed signal according to the statistical characteristics of the source signal when the source signal and the mixing method are unknown. There are two explanations for "blindness": one is that other information is unknown except for the observed signal, which is called "full blindness"; the other is that some prior information of the signal or hybrid system is known, which is called "semi-blindness". ". Aiming at the problem of mixed pixel decomposition, blind source separation technology obtains endmember spectrum and endmember abundance at the same time, that is, realizes blind decomposition of mixed pixels. Before dealing with the problem of mixed pixel decomposition, it is necessary to establish a spectral mixture model. The spectral mixing methods mainly include linear mixing and nonlinear mixing. The linear spectral mixing model is currently the most widely used mixed pixel decomposition model, and its expression is:

$X x = = {Σ Σ}_{i i = = 11}^{P P} {A A}_{i i} {s the s}_{i i} = = AS AS - - - - - - ((11 - - 11))$

式中，X＝[x₁,...x_l,...x_L]^T∈R^L×1为一个像元L个波段的混合像元光谱，是已知观测量；A∈R^L×P为端元光谱矩阵，其中每一列对应一个端元的光谱向量；P为端元数目，S∈R^P ^×1为该混合像元的丰度向量。丰度向量需要满足两个约束条件：非负性约束、和为1约束。即In the formula, X=[x ₁ ,...x _l ,...x _L ] ^T ∈ R ^L×1 is a mixed pixel spectrum of L bands in one pixel, which is a known observation; A∈R ^{L ×P} is the endmember spectral matrix, where each column corresponds to a spectral vector of an endmember; P is the number of endmembers, and S∈R ^P ^×1 is the abundance vector of the mixed pixel. The abundance vector needs to satisfy two constraints: the non-negativity constraint and the sum being 1 constraint. which is

s_i≥0,i∈{1,2,...,P} (2-1)s _i ≥ 0, i∈{1,2,...,P} (2-1)

${Σ Σ}_{i i = = 11}^{P P} {s the s}_{i i} = = 11,, i i &Element; &Element; {{1,2 1,2,, . . . . . .,, P P}} - - - - - - ((33 - - 11))$

若高光谱图像具有N个像元，则公式(1-1)可以进一步扩展为公式(4-1)。If the hyperspectral image has N pixels, formula (1-1) can be further extended to formula (4-1).

X＝AS (4-1)X=AS (4-1)

公式(4-1)中X∈R^L×N,其列向量表示每个像元在L个波段的混合光谱向量，S∈R^P×N构成了端元丰度矩阵。端元丰度矩阵满足以下约束条件：In formula (4-1), X∈R ^L×N , its column vector represents the mixed spectrum vector of each pixel in L bands, and S∈R ^P×N constitutes the endmember abundance matrix. The endmember abundance matrix satisfies the following constraints:

s_i,j≥0,i∈{1,2,...,P},j∈{1,2,...,N} (5-1)s _i,j ≥0,i∈{1,2,...,P},j∈{1,2,...,N} (5-1)

${Σ Σ}_{i i = = 11}^{P P} {s the s}_{i i,, j j} = = 11,, j j &Element; &Element; {{1,2 1,2,, . . . . . .,, N N}} - - - - - - ((66 - - 11))$

基于盲源分离技术的高光谱图像混合像元分解问题即在端元矩阵和丰度矩阵均未知的情况下，利用混合像元光谱的观测值来恢复端元光谱和端元丰度。The mixed pixel decomposition problem of hyperspectral images based on blind source separation technology is to use the observed values of the mixed pixel spectrum to recover the endmember spectrum and endmember abundance when both the endmember matrix and the abundance matrix are unknown.

高光谱图像的波段数多的特点使高光谱图像的数据量巨大，若采用传统的奈奎斯特采样定理对数据采集，将会对硬件的采样速度、传输速度以及存储能力产生巨大的压力。压缩感知理论的提出解决了这一难题，该理论指出当信号满足稀疏性或者可压缩时，能以远低于奈奎斯特采样率的速度对信号进行全局观测，然后采用压缩感知重构算法恢复源信号。高光谱图像在小波域是稀疏的，因此可以将压缩感知理论应用到对高光谱图像的采样中。具体的采样过程可以表示为下式：The large number of bands in hyperspectral images makes the data volume of hyperspectral images huge. If the traditional Nyquist sampling theorem is used to collect data, it will put a huge pressure on the sampling speed, transmission speed and storage capacity of the hardware. The compressed sensing theory solves this problem. The theory points out that when the signal satisfies sparsity or compressibility, the signal can be globally observed at a speed much lower than the Nyquist sampling rate, and then the compressed sensing reconstruction algorithm can be used. Restore the source signal. Hyperspectral images are sparse in the wavelet domain, so compressive sensing theory can be applied to the sampling of hyperspectral images. The specific sampling process can be expressed as the following formula:

Y＝ΦX^T＝Φ(AS)^T Y＝ΦX ^T ＝Φ(AS) ^T

在上述公式中，Φ∈R^M×N为观测矩阵,X∈R^L×N为混合像元光谱矩阵，Y∈R^M×L为压缩观测矩阵。In the above formula, Φ∈RM ^×N is the observation matrix, X∈R ^L×N is the mixed pixel spectrum matrix, and Y∈RM ^×L is the compressed observation matrix.

发明内容Contents of the invention

本发明是为了解决现有采用传统的奈奎斯特采样定理对高光谱图像数据采集时，混合像元分解速度慢的问题，本发明提供了一种基于压缩感知的高光谱混合像元分解的方法。The present invention aims to solve the problem of slow decomposition speed of mixed pixels when traditional Nyquist sampling theorem is used to collect hyperspectral image data. The present invention provides a hyperspectral mixed pixel decomposition method based on compressed sensing method.

一种基于压缩感知的高光谱混合像元分解的方法，该方法为，A method for decomposing hyperspectral mixed pixels based on compressed sensing, the method is,

步骤一，输入观测矩阵Φ和压缩观测矩阵Y，利用压缩感知理论，建立光谱混合模型：Step 1: Input the observation matrix Φ and the compressed observation matrix Y, and use the compressed sensing theory to establish a spectral mixing model:

Y＝ΦX^T＝Φ(AS)^T (1)Y＝ΦX ^T ＝Φ(AS) ^T (1)

Φ∈R^M×N为M×N的观测矩阵，R为实数，Φ∈R ^M×N is the observation matrix of M×N, R is a real number,

X∈R^L×N为L×N的混合像元光谱矩阵，X∈R ^L×N is L×N mixed pixel spectral matrix,

Y∈R^M×L为M×L的压缩观测矩阵，Y∈R ^M×L is the compressed observation matrix of M×L,

S∈R^P×N为P×N的端元丰度矩阵，S∈R ^P×N is the endmember abundance matrix of P×N,

A∈R^L×P为L×P的端元光谱矩阵，A∈R ^L×P is the endmember spectral matrix of L×P,

步骤二，初始化，随机选取一个端元光谱矩阵A作为端元光谱矩阵A的估计值且为L×P的矩阵，Step 2, initialization, randomly select an endmember spectral matrix A as the estimated value of the endmember spectral matrix A and is a matrix of L×P,

令端元丰度矩阵S的估计值：Let the estimated value of the endmember abundance matrix S be:

$\overset{^^}{S S} = = {[\begin{matrix} 00 & . . . . . . & 00 \\ . . . . . . \\ 00 & 00 \end{matrix}]}_{N N \times \times P P} - - - - - - ((22))$

其中，为N×P的矩阵，in, is an N×P matrix,

步骤三，令迭代次数变量t的初始值为1；Step 3, let the initial value of the iteration number variable t be 1;

步骤四，对端元丰度矩阵S的估计值和端元光谱矩阵A的估计值进行迭代处理；Step 4, the estimated value of the endmember abundance matrix S and an estimate of the endmember spectral matrix A perform iterative processing;

步骤五，在公式(1)两边同时乘以的伪逆，则公式(1)变形为Step five, multiply both sides of formula (1) by The pseudo-inverse of the formula (1) is transformed into

Y₁＝ΦS^T (3)，Y ₁ =ΦS ^T (3),

步骤六，设小波基为Θ，Σ为端元丰度的稀疏系数，则公式(3)变形为Step 6, set the wavelet base as Θ, and Σ as the sparse coefficient of endmember abundance, then formula (3) is transformed into

Y₁＝Φ(ΘΣ)^T (4)，Y ₁ = Φ(ΘΣ) ^T (4),

采用BP算法求解下式：The BP algorithm is used to solve the following equation:

argmin||Σ||₁ s.t. Y₁＝Φ(ΘΣ)^T (5)，argmin||Σ|| ₁ st Y ₁ = Φ(ΘΣ) ^T (5),

获得端元丰度矩阵S的估计值Obtain an estimate of the endmember abundance matrix S

$\overset{^^}{S S} = = - - - - - - ΘΣ ΘΣ ((66)),,$

步骤七，在概率密度函数为P(S_k)∝α exp(α|S_k|)的条件下，归一化端元光谱矩阵A的估计值的列向量，更新端元光谱矩阵A的估计值，Step 7, under the condition that the probability density function is P(S _k )∝α exp(α|S _k |), the estimated value of the normalized endmember spectral matrix A A column vector of , updating the estimated value of the endmember spectral matrix A ,

$\overset{^^}{A A} = = {A A}_{11} + + λ λ {{- - {A A}_{11} ((B B {\overset{^^}{S S}}^{T T} + + I I))}} - - - - - - ((77)),,$

其中，A₁为上一次迭代得到的端元光谱矩阵的估计值，λ为迭代步长，I为单位矩阵，B为向量集合，且B＝{B₁，B₂，B₃，......B_k}，为第k个源信号的先验概率分布，k为整数，α为大于零的实数，Among them, A ₁ is the estimated value of the endmember spectral matrix obtained in the last iteration, λ is the iteration step size, I is the identity matrix, B is the vector set, and B={B ₁ , B ₂ , B ₃ ,... ... _Bk }, is the prior probability distribution of the kth source signal, k is an integer, and α is a real number greater than zero,

步骤八，如果相邻两次获得的端元光谱矩阵A的估计值中对应的每个元素的绝对值之差小于0.1时，则停止迭代，执行步骤九，否则，令t＝t+1，且返回步骤四；Step 8, if the estimated value of the endmember spectral matrix A obtained twice adjacent When the difference of the absolute value of each corresponding element in is less than 0.1, then stop the iteration and execute step nine, otherwise, set t=t+1, and return to step four;

步骤九，输出端元丰度矩阵完成对高光谱混合像元的分解。Step 9, output the endmember abundance matrix Complete the decomposition of hyperspectral mixed pixels.

本发明带来的有益效果是，对端元丰度矩阵S的评价指标采用平均信噪比SNR_avg实现，The beneficial effect brought by the present invention is that the evaluation index of the endmember abundance matrix S is realized by using the average signal-to-noise ratio _SNRavg ,

${SNR SNR}_{avg avg} = = \frac{11}{P P} {Σ Σ}_{k k = = 11}^{P P} 2020 \times \times {log log}_{1010} ((\frac{| | | | {S S}_{k k} {| | | |}_{22}}{{| | | | {S S}_{k k} - - {\overset{^^}{S S}}_{k k} | | | |}_{22}})) - - - - - - ((88))$

实验中采用的高光谱图像的端元数目P＝3，波段数L＝40，高光谱的像元数N＝32×32。在压缩感知采样时，压缩率为0.8，对混合像元分解的端元丰度矩阵的平均信噪比为17.78dB，可以较好的实现分离，使得混合像元分解速度提高了40％以上。The hyperspectral image used in the experiment has the number of end members P=3, the number of bands L=40, and the number of hyperspectral pixels N=32×32. In compressed sensing sampling, the compression ratio is 0.8, and the average signal-to-noise ratio of the end-member abundance matrix decomposed into mixed pixels is 17.78dB, which can achieve better separation and increase the decomposition speed of mixed pixels by more than 40%.

具体实施方式detailed description

具体实施方式一：本实施方式所述的一种基于压缩感知的高光谱混合像元分解的方法，该方法为，Specific embodiment 1: A method for decomposing hyperspectral mixed pixels based on compressed sensing described in this embodiment, the method is:

Y＝ΦX^T＝Φ(AS)^T (1)Y＝ΦX ^T ＝Φ(AS) ^T (1)

其中，为N×P的矩阵，in, is an N×P matrix,

Y₁＝ΦS^T (3)，Y ₁ =ΦS ^T (3),

Y₁＝Φ(ΘΣ)^T (4)，Y ₁ = Φ(ΘΣ) ^T (4),

$\overset{^^}{S S} = = - - - - - - ΘΣ ΘΣ ((66)),,$

步骤七，在概率密度函数为P(S_k)∝α exp(α|S_k|)的条件下，归一化端元光谱矩阵A的估计值的列向量，更新端元光谱矩阵A的估计值 Step 7, under the condition that the probability density function is P(S _k )∝α exp(α|S _k |), the estimated value of the normalized endmember spectral matrix A A column vector of , updating the estimated value of the endmember spectral matrix A

Claims

1. A method for decomposing a hyperspectral mixed pixel based on compressed sensing is characterized in that,

inputting an observation matrix phi and a compressed observation matrix Y, and establishing a spectrum mixed model by using a compressed sensing theory:

Y＝ΦX^T＝Φ(AS)^T(1)

Φ∈R^M×Nan observation matrix of M × N, R is a real number,

X∈R^L×Na mixed pixel spectral matrix of L × N,

Y∈R^M×La compressed observation matrix of M × L,

S∈R^P×Nis an end-member abundance matrix of P × N,

A∈R^L×Pan end-member spectral matrix of L × P,

initializing, and randomly selecting an end member spectrum matrix A as an estimated value of the end member spectrum matrix AAnd isIs a matrix of L × P and,

let the estimate of the end-member abundance matrix S:

\hat{S} = {[\begin{matrix} 0 & ... & 0 \\ ... \\ 0 & 0 \end{matrix}]}_{P \times N} - - - (2)

wherein,is a matrix of P × N and,

step three, setting an initial value of an iteration variable t as 1;

step four, estimating the abundance matrix S of the end-to-end elementEstimate of sum-end spectral matrix ACarrying out iterative processing;

step five, multiplying the two sides of the formula (1) simultaneouslyThe pseudo-inverse of (1) is transformed into

Y₁＝ΦS^T(3)，

Step six, if the wavelet basis is theta and sigma is a sparse coefficient of the abundance of the end members, the formula (3) is deformed into

Y₁＝Φ(ΘΣ)^T(4)，

Solving the following formula by using a BP algorithm:

arg min||Σ||₁s.t. Y₁＝Φ(ΘΣ)^T(5)，

obtaining an estimate of an end-member abundance matrix S

\hat{S} = Θ Σ - - - (6),

Step seven, the probability density function is P (S)_k)∝αexp(α|S_k|) normalized end-member spectral matrix AUpdating the estimated value of the end-member spectral matrix A

\hat{A} = A_{1} + λ {- A_{1} (B {\hat{S}}^{T} + I)} - - - (7),

Wherein A is₁For the estimated value of the end member spectrum matrix obtained in the last iteration, λ is the iteration step length, I is the identity matrix, B is the vector set, and B ═ B₁，B₂，B₃，......B_k}， Is the prior probability distribution of the kth source signal, k being an integer, α being a real number greater than zero, S_kRepresents the k-th end-member abundance matrix S,represents an estimate of the k-th end-member abundance matrix S,

step eight, if the estimated value of the end member spectrum matrix A obtained by two adjacent timesWhen the difference between the absolute values of each corresponding element in (b) is less than 0.1, thenStopping iteration, executing the step nine, otherwise, making t equal to t +1, and returning to the step four;

ninth, output end member abundance matrixAnd (4) decomposing the hyperspectral mixed pixel.