高光谱图像处理

高光谱传感器在电磁光谱的可见光、近红外、中红外和短波红外部分产生数百个窄光谱带的图像。这意味着在每个点上都会获得整个频谱。这些丰富的信息在许多应用中都很有用，比如地质学和农业。

本节介绍OTB在高光谱图像处理中的应用。

拆分

由于高光谱传感器的低空间分辨率、微观物质混合和多次散射，所测得的光谱是场景中实际材料光谱的混合物。因此，假设像素是几种材质的混合物，称为端元。在本节中，假设混合是线性的，即

R=A.S+N

在哪里，如果 l 是高光谱图像的波段数， n 像素数和 k 端子成员的数量：

  • R is the matrix of observed pixel, of size l.n
  • A is the endmember matrix, of size l.k
  • S is the abundance matrix, of size k.n
  • N is the noise matrix, of size l.n

分解问题是估计矩阵的问题 A 和 S 从… R 。下面是由几种地质材料组成的场景中的高光谱分解示例。该场景是AVIRIS传感器获取的Cuprite数据集的摘录。摘录是可用的 here 并且整个数据集可以从 AVIRIS 美国宇航局网站。

../_images/cuprite_rgb.png

Fig. 3 考虑到铜矿图像，这里显示了波段16、100和180。

由于输入图像的端元数目是未知的，高光谱分解的第一步是估计这个数目。这里使用了HFC虚拟维度算法。

otbcli_EndmemberNumberEstimation  -in inputImage.tif
                                  -algo vd
                                  -algo.vd.far 1e-5

该算法利用协方差矩阵的特征值与相关矩阵的特征值之差的Neyman-Pearson统计检验来估计端元数。如果一个分量的特征值之间的差为零，这意味着除了该特定分量的噪声之外，没有端元对相关特征值有贡献，因为噪声能量由协方差特征值表示。该算法的输出为：

Output parameters value:
number: 19

高光谱的下一步是估计端元。这可以使用顶点分量分析算法(VCA)来完成。

otbcli_VertexComponentAnalysis  -in inputImage.tif
                                -ne 19
                                -outendm endmembers.tif

应用程序的输出是一个包含以下内容的图像 k 像素，每个像素表示一个端元。

该算法基于对数据中端元的研究。这意味着每个末端成员必须至少关联一个纯像素。这个想法背后的基本原理是，高光谱混合数据包含在维度的单纯形中 k ，其中端元是单形的顶点。该算法具有算法复杂度低、无噪声情况下端元估计无偏等优点，在高光谱分解中得到了广泛的应用。然而，如果不尊重纯像素假设，则会对端元产生估计误差。

分解的最后一步是使用优化算法估计丰度矩阵。

otbcli_HyperspectralUnmixing  -in inputImage.tif
                              -ie endmembers.tif
                              -out unmixedImage.tif
                              -ua ucls

这里使用了一种无约束最小二乘算法。生成的丰度图如下所示：

../_images/hyperspectralUnmixing_rgb.png

Fig. 4 生成的未混合图像，这里显示的是前三个波段。

分类

高光谱图像处理中的一项常见任务是对数据立方体进行分类。例如，可能需要将像素与参考端元进行匹配。光谱角度映射器(SAM)算法 [1] 这是通过计算每个像素和参考之间的光谱测量来实现的。

像素的光谱角度 x 使用参考像素 r 由以下因素定义：

sam[x, r] = cos^{-1}(\frac{<x,r>}{\|x\| * \|r\| } )

哪里 <x,r> 表示x和r之间的标量积。这也称为 x 和 r 。在SAM分类中，计算具有一组参考像素的每个像素的光谱角度，与该像素相关联的类别是具有最低光谱角度的参考像素的索引。

otbcli_SpectralAngleClassification -in inputImage.tif
                                   -ie endmembers.tif
                                   -out classification.tif
                                   -mode sam

分类结果如下：

../_images/classification.png

Fig. 5 由此产生的未混合的机密图像。

在该应用中，基于光谱信息发散的另一种算法是可用的 [1] ，如果我们通过以下公式定义像素x的概率质量函数p：

x = [x_1, x_2 ,..., x_L ]^T \\ p = [p_1, p_2 ,..., p_L ]^T \\ p_i = \frac{x_i}{\sum_{j=1}^L x_j}

像素之间的光谱信息发散 x 和参考资料 r 由以下因素定义：

sid[x, r] = \sum_{j=1}^{L} p_j *log(\frac{p_j}{q_j}) + \sum_{j=1}^{L} q_j * log(\frac{q_j}{p_j})

其中p和q分别是x和r的概率质量函数。与SAM算法一样，与像素相关联的类别是具有最低光谱信息发散度的参考像素的索引。

请注意， classification recipe 也可以应用于高光谱数据。

异常检测

异常是指在场景中不会发现的元素。这种不同寻常的元素很可能与它的环境不同，它的存在是少数人的场景。通常，田野中的岩石和森林中的木屋都是异常现象，可以使用高光谱图像进行检测。这里将使用一种异常检测算法--局部Rx检测器来检测城市环境中的小目标。输入图像已采集完毕 Pavia 大学的ROSIS传感器，有103个光谱波段。

../_images/pavia.png

Fig. 6 帕维亚大学的高光谱图像，这里显示了50、30和10波段。

第一个可选步骤是降低输入图像的维度，以便在保持与异常相关的信息的同时减少光谱数据的大小，这里使用主成分分析算法。

otbcli_DimensionalityReduction  -in inputImage.tif
                                -out reducedData.tif
                                -method pca
                                -nbcomp 10

由于局部Rx需要计算每个像素上的相关矩阵并求逆，因此将降维作为预处理步骤将显著降低算法的计算代价。

局部Rx检测使用滑动窗口来计算每个像素的异常分数。该窗口由两个RADIUS子窗口组成 ir 和 er 同 ir < er ，并以感兴趣的像素为中心。通过比较中心像素和属于环形的像素来计算局部Rx分数。

otbcli_LocalRxDetection  -in reducedData.tif
                         -out RxScore.tif
                         -ir 1
                         -er 5

这里的异常被认为是小的，因此应该选择小的内半径。 ir=1 。另外，由于环境是城市的，由于传感器的空间分辨率较低，背景的统计量随距离变化很快，所以外半径不宜太大，在这里 er=5 已经被选中了。这些参数实际上取决于输入图像和感兴趣的对象。

然后，可以对得到的分数图像应用阈值，以产生异常图。

otbcli_BandMath  -il RxScore.tif
                 -out anomalyMap.tif
                 -exp "im1b1>100"

阈值的值取决于异常检测器的灵敏度。

image_1 image_2
左：计算的处方分数，右：检测到的异常(红色)

[1] Du, Yingzi & Chang, Chein-I & Ren, Hsuan & Chang, Chein-Chi & Jensen, James & D'Amico, Francis. (2004). New Hyperspectral Discrimination Measure for Spectral Characterization. Optical Engineering - OPT ENG. 43.