多光谱图像可以为森林定位提供数据信息，已经用于生态预测、环境政治和城市发展等相关研究^[1-2]。由于外界环境的影响和传感器的局限性，原始多光谱图像中普遍存在大量的混合像素，使得多光谱图像的分辨率变得模糊。图像分辨率的模糊性给包括森林在内的土地覆盖类型的空间分布带来了诸多挑战^[3-5]。虽然光谱解混技术可以估计出像素属于每个地物类别的比例，但是仍然不能准确获得地物类别的分布情况^[6-9]。为了解决这个问题，提出了亚像素定位(subpixel mapping, SPM)技术^[10-12]。它处理丰度图像(原始图像的光谱解混结果)得到亚像素级的精细分辨率下的地物类别分类图像。

近年来，关于SPM的研究成果越来越多。文献[13-14]提出了基于像素交换的SPM算法，在混合像素内交换2个需要交换的亚像素的位置，然后迭代逼近SPM结果。文献[15-16]提出了基于Hopefield神经网络的SPM算法，该方法基于能量最小化原理，应用Hopefield神经网络获得代表亚像素的神经元输出高分辨率丰度图像，进而获得SPM结果。除此之外，可以利用BP神经网络、指示协克里格、空间-光谱相关性、时空信息和一些超分辨率改善SPM效果^[17-23]。特别地，可以采用不同类型的空间吸引模型(space attraction model, SAM)对空间相关性进行量化，得到SPM结果，由于考虑到两尺度空间相关性，混合空间吸引模型(mixed space attraction model, MSAM)^[24-25]可以得到更精确的定位结果。

SPM技术已成功用于获取森林区域定位，命名为亚像素森林定位(subpixel forest mapping, SPFM)^[26]。特别地，由于物理意义简单，不需要任何先验信息，基于MSAM的SPFM(SPFM-MSAM)是一种常用的获得森林定位结果的方法。然而，SPFM-MSAM中的森林光谱信息往往得不到充分利用，影响了最终定位结果的精度。为了改进SPFM-MSAM，本文提出了一种面向空-谱信息的多光谱图像亚像素森林定位方法(subpixel forest mapping for multispectral image based on space-spectrum information, SPFM-SSI)。SPFM-SSI利用MSAM生成的空间项和归一化差分植被指数(normalized difference vegetation index, NDVI)^[26]生成的光谱项共同提高定位精度。

1 多光谱图像亚像素森林定位方法(SPFM-SSI) 1.1 空间项

在提出的SPFM-SSI中，根据空间相关性原理，将空间项T^SPA定义为最大优化问题。MSAM用于计算中心亚像素p_j(j=1，2, …, MS²; M是混合像素的数目; MS²是亚像素数目)和相邻的混合像素或相邻的亚像素之间的空间相关性。空间项T^SPA的数学模型定义为：

$ \begin{gathered} T^{S P A}=\max \sum\limits_{j=1}^{M S^{2}} o_{j} \times\left[\sum\limits_{C=1}^{N} w_{j C} \times F\left(P_{C}\right)+\sum\limits_{c=1}^{N^{*}} w_{j c}\right]+\\ \left(1-o_{j}\right) \times\left[\sum\limits_{C=1}^{N} w_{j C} \times\left(1-F\left(P_{C}\right)\right)+\sum\limits_{c=1}^{N^{*}} w_{j c}\right] \end{gathered} $

(1)

式中：F(P_C)表示通过光谱解混获得的第C个与中心亚像素p_j相邻的混合像素P_C属于森林类别的概率值；N是相邻混合像素的数目，N^*是相邻亚像素的数目，在本文中，它们都被认为是最多8邻域；o_j是每个亚像素的二进制标签(1表示森林标签，0表示背景标签)。w_jC用于计算中心亚像素p_j和相邻混合像素P_C之间的相关性，w_jc用于计算中心亚像素p_j和相邻亚像素p_c之间的相关性，如图 1所示为这2种尺度的空间相关性示意图。

w_jC和w_jc被定义为：

$ w_{j C} =\exp \left(-d\left(p_{j}, P_{C}\right)^{2} / r_{1}\right) $

(2)

$ w_{j c} =\exp \left(-d\left(p_{j}, p_{c}\right)^{2} / r_{2}\right) $

(3)

式中：r₁和r₂为非线性参数；d(p_j, P_C)是中心亚像素p_j与相邻混合像素P_C之间的欧几里德距离；d(p_j, p_c)是中心亚像素p_j与相邻亚像素p_c之间的欧几里德距离。根据空间相关性原理，在最大化T^SPA的前提下，可以得到较好的亚像素分布。由于MSAM考虑了2种尺度的空间相关性，产生的空间项T^SPA具有更加丰富的空间信息。

1.2 光谱项

为了更充分地利用多光谱图像中森林的光谱信息，本文采用NDVI，提出了一种新型的具有光谱信息的光谱项T^SPE。将森林光谱反射率从红色波段(red band, RB)到近红外波段(near-infrared band, NB)的急剧增加的情况考虑进NDVI中。因此通过最小化观测到的NDVI值(NDVI^OBE)和模拟的NDVI值(NDVI^SIM)之间的光谱差导出T^SPE。

NDVI^OBE被定义为：

$ \mathrm{NDVI}^{\mathrm{OBE}}=\sum\limits_{1}^{M} \frac{r_{\mathrm{NB}}-r_{\mathrm{RB}}}{r_{\mathrm{NB}}+r_{\mathrm{RB}}} $

(4)

式中：r_NB和r_RB是2个波段中每个混合像素的观测光谱反射率; M是混合像素的数量。

通过计算NB中各混合像素的模拟光谱反射率o_NB和RB中各混合像素的模拟光谱反射率o_RB，得到NDVI^SIM：

$ \mathrm{NDBI}^{\mathrm{SIM}}=\sum\limits_{1}^{M} \frac{o_{\mathrm{NB}}-o_{\mathrm{RB}}}{o_{\mathrm{NB}}+o_{\mathrm{RB}}} $

(5)

每个混合像素的模拟光谱反射率被认为是其内所有亚像素光谱值的线性混合。o_NB和o_RB的数学模型为：

$ o_{\mathrm{NB}}=\left(o_{\mathrm{NB}}^{\mathrm{F}} \times h_{\mathrm{NB}}^{\mathrm{F}}\right)+\left[o_{\mathrm{NB}}^{\mathrm{BC}} \times\left(1-h_{\mathrm{NB}}^{\mathrm{F}}\right)\right] $

(6)

$ o_{\mathrm{RB}}=\left(o_{\mathrm{RB}}^{\mathrm{F}} \times h_{\mathrm{RB}}^{\mathrm{F}}\right)+\left[o_{\mathrm{RB}}^{\mathrm{BC}} \times\left(1-h_{\mathrm{RB}}^{\mathrm{F}}\right)\right] $

(7)

式中：o_NB^F和o_RB^F分别是NB和RB中森林的光谱反射率；o_NB^BG和o_RB^BG是2个波段中的背景反射率；h_NB^F和h_RB^F是2个波段中森林在每个混合像素中的比例，它们是通过森林亚像素数除以总亚像素数得到的。则背景所占的比例由1-h_NB^F和1-h_RB^F给出。

因此，光谱项T^SPE被表述为：

$ T^{\mathrm{SPE}}=\min \left(\mathrm{NVBI}^{\mathrm{OBE}}-\mathrm{NVBI}^{\mathrm{SIM}}\right)^{2} $

(8)

1.3 优化函数

SPFM-SSI的目标是通过一个权重参数δ(0≤δ < 1)最小化包括空间项T^SPA和光谱项T^SPE的优化函数E：

$ \min E=\delta T^{\mathrm{SPE}}-(1-\delta) T^{\mathrm{SPA}} $

(9)

为了得到更好的映射结果，本文采用了模拟退火算法优化函数E。首先，每个亚像素被随机分配森林标签或背景标签。然后迭代地改变亚像素的标签，直到得到E的最小值。在每次迭代中，森林标签都会更改为背景标签，反之亦然。如果E减少，则接受此更改，否则拒绝此更改。如果需要更改的标签少于总数的0.1%，则停止更改，完成模拟退火算法处理，得到最终的亚像素森林定位结果。

本文提出的SPFM-SSI通过空间项T^SPA和光谱项T^SPE充分利用了空-谱信息，因此使得最终的森林定位结果精度得到提高。

2 实验与分析 2.1 实验设置

选取2幅来自巴西亚马逊雨林的Landsat 8 OLI多光谱图像作为实验数据。图像2的森林分布比图像1复杂。2幅800×800像素图像的测试区域如图 2所示。为了定量地估计SPFM-SSI的性能，利用亚像素定位最常用的实验设置，模拟的粗糙多光谱图像是通过对原始的精细多光谱图像通过比例尺度S进行降采样获得的。在这种情况下，可以知道亚像素级的土地覆盖类型，因此可以评估图像配准误差对本文所提出方法的影响。如图 3所示，通过S=10降采样产生模糊图像。采用基于最小二乘支持向量机的光谱解混方法得到丰度图像^[4]。通过多次参数测试，对于2幅图像选择了参数δ=0.6和δ=0.5。

本文比较了3种SPFM方法：基于亚像素吸引模型的SPFM-SPSAM、SPFM-MSAM和提出的SPFM-SSI。利用总体精度(overall accuracy, OA (%))和kappa系数(Kappa)对3种SPFM方法进行评价。使用Matlab 2018a软件平台进行实验。

2.2 结果分析

如图 4(a)和5(a)所示，2个测试图像的参考图像是对图 2(a)、(b)精细多光谱图像通过基于最小二乘支持向量机分类算法获得的。从图 4和图 5所示的数据1和数据2的4个SPFM结果可以看出，利用SPFM方法可以获得森林的空间分布。但是由于SPFM-SPSAM和SPFM-MSAM没有充分利用原始图像的光谱信息，使得结果中边界处也存在许多毛刺，平滑区存在许多断孔。由于考虑了更丰富的空-谱信息，可以观察到SPFM-SSI产生了更加理想的结果，并且更加接近参考图像。

表 1显示的3种方法的评价指标，可以观察到提出的SPFM-SSI各项评价指标均高于其他SPFM方法。例如对于图像1得到的结果中的OA，SPFM-SPSAM产生的OA值为95.42%，SPFM-MSAM产生的OA值为95.77%。与SPFM-MSAM相比，SPFM-SSI使OA值增加约1.0%。根据OA的定义，图像1包含800×800个像素，提高了约1.0%意味着将会增加6 400个像素被正确定位，因此所提出的方法可以明显改森林定位的精度。与图像1得到的结果相似，针对图像2中的各项评价指标。SPFM-SSI方法仍然可以得到最佳的森林定位结果。

此外，为了测试比例尺度S对SPFM的影响，对图像1中进行另外2个比例尺度S=8和S=20的测试。图 6所示为3种比例尺度S下的3种SPFM方法的OA和Kappa。与表 1中的结果类似，SPFM-SSI仍然可以得到最高的OA和Kappa值。

最后，研究所选参数δ对SPFM-SSI的影响。以0.1为间隔，从0到0.9，10个组合对2个测试图像进行测试。总体精度OA与参数δ的关系如图 7所示。由于参数δ是平衡空间项T^SPA和光谱项T^SPE对E的影响，因此选择合适的δ可以改善最终的森林定位结果。针对2个测试图像，当分别δ=0.6和δ=0.5时，可以得到最高的总体精度OA，即此时为更合适的参数δ值。

3 结论

1) 利用混合空间引力模型充分提取多光谱图像多尺度空间信息，同时使用归一化差分植被指数充分利用多光谱图像森林光谱信息；

2) 通过充分地利用空-谱信息，SPFM-SSI提高了亚像素森林定位效果。

3) 对2幅Landsat 8 OLI卫星图像进行了测试，结果表明提出的SPFM-SSI较现有的SPFM方法得到了获得的亚像素森林定位精度更高。

本文通过多次试验选择了合适的参数δ。因此，一种自适应的选择参数δ方法是值得研究的。此外，SPFM-SSI在大面积真实图像中的性能值也得进一步研究。