EN中文

Journal of Geo-information Science >

2023 , Vol. 25 >Issue 9: 1894 - 1907

DOI: https://doi.org/10.12082/dqxxkx.2023.220938

Construction of Sparse Non-Photosynthetic Vegetation Index Combining Sentinel-2 Data and Spectral Features

  • YUAN Hao , 1 ,
  • JI Cuicui , 1, * ,
  • YANG Xuemei 2 ,
  • CHEN Maolin 1 ,
  • PAN Jianping 1 ,
  • CAO Yiming 1
Expand
  • 1. Smart City Academy, Chongqing Jiaotong University, Chongqing 400074, China
  • 2. School of Tourism, Lanzhou University of Arts and Sciences, Lanzhou 730000, China
* JI Cuicui, Email:

Received date: 2022-12-01

  Revised date: 2023-03-26

  Online published: 2023-09-05

Supported by

National Key R&D Program of China(2020YFC1511602)

National Natural Science Foundation of China(32060373)

National Natural Science Foundation of China(41801394)

Natural Science Foundation of Gansu Province(20JR10RA465)

Fold

Abstract

Non-photosynthetic Vegetation (NPV) is a key component in the ecosystem of arid and semi-arid regions and plays a vital role in the carbon cycle of vegetations. Obtaining the quantitative information of NPV is of great significance for scientific assessment of land desertification and scientific prevention and control of desertification. This study takes Minqin County, Gansu Province as the research area, analyzes the NPV and other end-component spectral information collected by the hyperspectral instrument in the research area, and constructs the vegetation index for extracting sparse NPV coverage based on Sentinel-2 MSI image. The measured NPV coverage is used as the verification data, and the coefficient of determination (R2), root mean square error, and significant p-value are used for model accuracy evaluation. The results show that: ① Based on the spectral analysis using the original spectral method, the first derivative transformation method, the reciprocal logarithm method, and the continuous elimination method, it is found that the envelope value of the NPV spectrum obtained by the continuous elimination method is significantly different from that of Photosynthetic Vegetation (PV) and Bare Soil (BS); ② After calculating the envelope difference of the green light band and the spectral reflectance ratio of the short wave infrared band, the difference between the NPV and PV and BS spectral eigenvalues can be effectively detected, in order to extract the NPV information and obtain a new index representing the NPV coverage information; ③ Through precision verification, it is found that the newly constructed NPV index-Envelope Difference Vegetation Index (EDVI) has the best correlation with the measured NPV coverage on the ground (R2=0.71, RMSE=0.0170, p<0.01), followed by the four-band red edge soil adjusted NDTI (S_NDTI4RE) (R2 is 0.57 with weight γ=0). While other NPVIs models have poor accuracy in estimating sparse NPV coverage. Our results indicate that the EDVI index model constructed in this paper can effectively achieve the rapid extraction and monitoring of NPV at large scale in arid and semi-arid regions through spectral feature analysis.

Key words： non-photosynthetic vegetation; non-photosynthetic vegetation index; Sentinel-2 MSI; arid and semi-arid regions; spectral analysis; linear regression; feature space

Cite this article

YUAN Hao , JI Cuicui , YANG Xuemei , CHEN Maolin , PAN Jianping , CAO Yiming . Construction of Sparse Non-Photosynthetic Vegetation Index Combining Sentinel-2 Data and Spectral Features[J]. Journal of Geo-information Science, 2023 , 25(9) : 1894 -1907 . DOI: 10.12082/dqxxkx.2023.220938

1 引言

近年来，由于气候变化以及人类不合理的生产生活活动使得占地球陆地面积40%的干旱和半干旱区的土地荒漠化问题日趋严重[1]；如何及时、准确的监测干旱和半干旱区的植被分布状况，为科学评估土地荒漠化程度以及荒漠化科学防治提供科学依据是目前的一个重要议题[2-3]。植被按照能否进行光合作用可以分为光合植被（Photosynthetic Vegetation， PV）和非光合植被(Non-Photosynthetic Vegetation， NPV)，其中NPV主要包括植被的凋落物、枯叶以及非绿色枝干等。NPV作为干旱区生产力的重要组成部分，在维护干旱与半干旱区生态系统的结构、功能和动态平衡方面发挥了重要作用[1,4]。至今，光合植被的定量提取研究已经取得了较大的进展，而NPV的定量提取研究，特别是在干旱与半干旱地区的研究相对较少[5]，因此，快速、准确的非光合植被覆盖度（Fractional Coverage of Non-Photosynthetic Vegetation， fNPV）定量估算对调查干旱区植被状况具有重要意义[6]。
传统的fNPV获取方法大多为小范围的植被覆盖度调查，费时、费力、成本高[7-8]，不能满足干旱区大面积fNPV的估算需求。遥感技术的发展为大范围快速、准确地获取fNPV信息以及进行连续观测提供了有效的技术手段[9]。目前，基于遥感数据的fNPV定量估算主要有3种方法：① 回归模型法，建立fNPV与非光合植被指数或者光谱反射率之间的线性或非线性关系[10]，利用特定区域的光谱特征来检测混合像元中的成分；② 光谱混合分析法（Spectral Mixture Analysis，SMA ），假设每个混合像元由不同的“端元”构成，由混合像元光谱特征计算出各个端元的权值[11]；③ 数据驱动模型，使用机器学习方法挖掘空间信息[12]。其中，回归模型应用简单、实用性强，与植被的光谱特征联系紧密，被很多学者作为fNPV的遥感估算方法，并取得了一定的进展。
现有的NPV指数当中，基于高光谱数据构建的高光谱纤维素吸收指数(Cellulose Absorption Index, CAI) 利用NPV在短波红外波段特有的木质-纤维素吸收特征将NPV和裸土（Bare Soil, BS）进行区分，且其在玉米、大豆、作物茬的混合场景中被证明与作物茬覆盖度有较好的线性相关性（R2=0.88）[4]。但由于高光谱数据获取困难，使其无法用于fNPV的大面积估算。多光谱数据虽然缺乏详细的光谱特征，依靠海量数据、覆盖范围大等优势，广泛应用于fNPV的估算研究[3,13]。基于多光谱传感器构建的指数中，归一化差异指数5（Normalized Difference Index 5, NDI5）和NDI7指数利用Landsat近红外和中红外波段反演农田玉米残留物覆盖度[14]，但其仅能对覆盖度进行分类而无法精确估算，且当分类类别较多时精度下降较大。同样基于Landsat数据的归一化差异耕作指数（Normalized Difference Tillage Index, NDTI）、土壤耕作指数（Soil Tillage Index, STI）则利用短波红外波段实现了耕地场景作物残茬的覆盖估算，且取得了较好的效果[15]。归一化差异衰老植被指数（Normalized Difference Senescent Vegetation Index, NDSVI）将衰老植被组织中水分流失的特征用于识别牧草中的衰老植被，相关性系数达到了0.9[16]。基于MODIS数据B6和B7波段的比值指数（Shortwave Infrared Ratio, SWIR32）以及干枯燃料指数（Dead Fuel Index, DFI），前者与归一化植被指数（Normalized Difference Vegetation Index, NDVI）结合构建的三元线性混合模型较为准确地估算了澳大利亚草原地区的fNPV[17]，后者在模拟的混合场景fNPV的估算中表现出了较好的相关性(R2=0.94)[18]。基于Sentinel-2数据的改进的4波段红边NDTI（4-band Red Edge NDTI， NDTI4RE）、4波段红边S_NDTI（4-band Red Edge Soil Adjusted NDTI，S_NDTI4RE）、4波段红边STI（4-band red edge STI，STI4RE）[19]被证明在青藏高原地区较NDTI、S_NDTI、STI有更好的效果。在上述指数构建选用的传感器中， Landsat TM升空时间最早并得到了广泛的应用，具有较高的空间分辨率。与之相比，于1999年升空的MODIS则在时间分辨率上更具有优势，但受限于其空间分辨率，使其更适用于大尺度的动态监测，而非县域大小的NPV信息提取。于2015年升空的Sentinel-2 MSI则兼顾时间于空间分辨率，与前两者相比，具有最多的波段数和最高的分辨率，且在光学数据中，Sentinel-2 MSI数据是唯一一个在红边范围含有3个波段的数据，这使得Sentinel-2数据在获取PV和NPV指数方面具有明显优势[20]。因此，本文选用Sentinel-2多光谱数据进行实验。由于干旱与半干旱地区植被覆盖度低，土壤背景对NPV反射率的影响较大[21]。因此，发展一种对干旱与半干旱地区提取非光合植被信息有效、且对背景干扰不敏感的非光合植被指数至关重要。
本研究主要开展两方面的工作：① 对NPV光谱进行分析，了解稀疏NPV信息提取有效光谱特征；② 基于Sentinel-2 MSI影像发展适用于干旱与半干旱地区稀疏NPV信息提取的新的非光合植被指数模型，与常用非光合植被指数（Non-Photosynthetic Vegetation Index, NPVI）进行精度比较分析，评价新构建模型可靠性。

2 研究区概况及数据来源

2.1 研究区概况

研究区位于甘肃省民勤县西缘绿洲——荒漠过渡区域内(38°37′42.60″N, 102°55′11.25″E)。研究区全境大部分被腾格里沙漠和巴丹吉林沙漠覆盖，属于温带干旱气候区，具有典型的温带大陆性气候特点，区域内土壤干燥、大气湿度较低，气温较高，降雨量低，蒸发量大，植被以红柳、白刺、梭梭为主，呈点状零散分布（图1(c)），植被覆盖稀疏，种群少，PV与NPV混杂共存（图1(c)）。由于本文的主要研究内容为PV、NPV、BS混合场景下的NPV信息提取，需要选择一年中植被覆盖度较高且较为稳定的时间段获取实验数据，所以选择在8月开展前期的研究区野外调研及数据采集工作。
图1 研究区样地分布及其对应非光合植被覆盖度

Fig. 1 The distribution of sample plots in the study area and their fractional coverage of non-photosynthetic vegetation coverage

2.2 实验数据

2.2.1 遥感影像数据获取及预处理

遥感影像数据采用Google Earth Engine（GEE）平台的Sentinel-2 MSI地面反射率数据集，选取与地面验证数据相匹配的2019年8月的所有影像，并利用具有云掩码信息的QA60波段筛选云量较小的影像进行镶嵌。由于研究过程中使用的部分波段的空间分辨率不同，需要先将空间分辨率为20 m的RE1、RE2、RE3、RE4、SWIR1以及SWIR2六个波段重采样到10 m。Sentinel-2 MSI波段参数见表1。
表1 Sentinel-2 MSI波段参数

Table 1 Sentinel-2 MSI Band Parameters Table

波段编号 波段名 波长范围/nm 空间分辨率/m
B1 AEROSOLS 442.3~443.9 60
B2 BLUE 492.1~496.6 10
B3 GREEN 559.0~560.0 10
B4 RED 664.5~665.0 10
B5 RE1 703.8~703.9 20
B6 RE2 739.1~740.2 20
B7 RE3 779.7~782.5 20
B8 NIR 833.0~835.1 10
B8A RE4 864.0~864.8 20
B9 WATER VAPOR 943.2~945.0 60
B11 SWIR1 1 610.4~1 613.7 20
B12 SWIR2 2 185.7~2 202.4 20

2.2.2 野外地面数据获取

野外实地数据采集包括样地内植被（PV、NPV）以及裸土的端元光谱数据采集和植被覆盖度采集两部分内容，两部分工作同期进行，具体实施方案如下。
端元光谱数据采集采用ASD(Analytical Spectral Devices)Field Space Pro光谱仪采集各端元350~2 500 nm波谱范围的光谱数据。其中，采集了红柳、红砂等10种类型灌丛植被的PV和NPV组分光谱，刺蓬、五星蒿等10种类型灌草的PV端元光谱，粗沙、细沙、结皮等5种类型的BS端元光谱。端元光谱数据均为天气状况良好、阳光充足的时段采集的纯净端元冠层光谱， PV端元为纯绿色植被植株，在下面铺上全吸收黑布，去掉BS及枯落物干扰；NPV端元为完全枯死植株下层铺黑布，测量其纯净端元光谱，对于没有纯净PV和NPV的植株，将植被纯净叶子和枝干分别裁剪下来分开放在黑布上，分别测量纯净叶子和纯净枝干端元光谱；BS为完全暴露裸地在其上方10 cm测量端元光谱（图2(a)、 图2(b)）。所有光谱数据测量5次并取其平均值作为最终结果，并去掉受大气影响严重波段和水汽吸收波段，保留350~1 350 nm、1 450~1 750 nm、2 000~2 350 nm波谱范围的数据。
图2 地面数据采集

Fig. 2 Ground data collection

非光合植被覆盖度验证数据获取于2019年8月，在研究区内共选取了50块植被分布均匀、植被覆盖情况稳定、混合物种的自然植被群落作为样地，样地分布及其非光合植被覆盖如图1所示。为使采集的样地数据能够与遥感影像的单个像素相对应，采用横断面取样方法布设样地[22]：在样地内沿正南正北方向用卷尺拉出一条30 m的样线，并以15 m刻度处为中心，再分别拉两条30 m样线，三条线间成60°角。采用针刺法(间隔1 m)在整刻度线处分别获取样地PV和NPV的覆盖度信息，同时利用GPS记录样地中心位置（图2(c)）。

3 研究方法

3.1 端元光谱分析

本文分别采用原始光谱法、一阶导数变换法、倒数的对数法以及连续统去除法对采集的端元光谱数据进行处理和分析，从而掌握NPV、PV和BS的端元光谱特征，筛选出NPV与PV和BS光谱特征明显不同波段，为NPV指数模型构建做准备。
（1）原始光谱法。将预处理过后的光谱数据计算反射率并分为PV、NPV、BS 共3组（图3(a)—图3(c)），图中不同类型的PV、NPV、BS虽然在反射率数值上有所差异，但光谱曲线都具有相似的变化趋势， PV组光谱均具有明显的“两谷一峰”及“红边”特征， NPV组光谱反射率上升和下降的幅度则较为缓和，BS组光谱曲线较为集中且反射率呈平缓上升状态，可见各组中均值曲线均能够较好的体现出组内不同类别端元的光谱特征，此外，由于研究区内各物种呈混合分布的特征，单一物种的光谱对于PV、NPV、BS整体而言不具有代表性，所以，在后续的研究中将采用均值光谱用于数据的处理及分析。
图3 端元光谱反射率及斜率曲线

Fig. 3 End member spectral reflectance and slope curve

（2）一阶导数变换法（式（1））。通过将原始光谱数据求一阶导进而得到导数光谱的非线性数学分析方法。其优势在于能从重叠的吸收光谱中分离出各端元的吸收峰，并能对其定量，能够减弱大气效应、土壤成分、太阳高度角和地形等外界因素导致的反射率变化不敏感问题[23-24]。
R ' ( λ i ) = R λ ( i + 1 ) - R λ i Δ λ
式中： R ' ( λ i )表示 λ i波长处的一阶导数值； λ i表示 i点的波长。
将均值光谱数据重采样到Sentinel-2的10个中心波段上，并计算其曲线斜率。如图3(d)所示，在BLUE至RE4波段之间， PV与NPV一阶导数具有相似的变化趋势，但PV的值更大，且在GREEN至RED波段，PV的一阶导数值为负值也是PV区别于另外2个端元的主要区别之一，RE4至SWIR2波段，NPV与BS具有一定的相似性，NPV一阶导数值变化起伏更加明显。
（3）倒数的对数法（式（2））。倒数的对数法是将光谱数据先取倒数再取对数进而得到光谱线性变换的方法，应用此方法能够将光谱数据的差别进行扩大，进而扩大明显差别部分，有利于不同端元的分别[23]。
R ″ ( λ i ) = l o g 10 1 R λ i
式中： R ″ ( λ i )表示 λ i波长处的倒数的对数值； λ i表示 i点的波长。
将倒数的对数法处理过后的端元光谱数据以及原始端元光谱数据的各个波段分别进行比值和差值处理（图4）。总的来说，无论是原始光谱还是倒数的对数光谱PV绝大部分波段之间的差值和比值差异均要明显大于另外两端元，因此在差值色块中颜色显示较深的色块更多，在比值色块中颜色最深和最浅的色块都更多。BS整体波段间差异性较为均匀，色块图颜色过渡平缓，没有明显的颜色差异；NPV在大多数波段差异性与PV保持一致，但在原始光谱的SWIR1-BLUE差值以及倒数的对数光谱的GREEN/SWIR1比值与BS与NPV有较大的不同。
图4 端元原始光谱曲线和倒数的对数光谱曲线波段间差值和比值

Fig. 4 The difference (upper left) and ratio (lower right) between the bands of the original spectral curve of the end member and the logarithmic spectral curve of the reciprocal

（4）连续统去除法。连续统去除法常用于植被高光谱吸收特征研究，该方法优点是对反射率光谱曲线变化较小或较为接近的可见光范围波段处理极其有效，可最大程度增强植被光谱间差异，易于对不同植被光谱曲线的特征值进行对比，有利于光谱的识别分类[23,25]。本文将该方法用于处理多光谱数据，以连续的3个波段为一组，分别以BLUE和RED、GREEN和RE1、RED和RE2、RE1和RE3、RE2和NIR、RE3和RE4、NIR和SWIR1、RE4和SWIR2为起终点，计算GREEN、RED、RE1、RE2、RE3、NIR、RE4以及SWIR1共8个波段的包络线值（式（3）—式（5））。
C R i = R i / R H i
R H i = R s t a r t + k ( λ i - λ s t a r t )
k = R e n d - R s t a r t / λ e n d - λ s t a r t
式中： k为连接波峰或者波谷起终点直线（HULL）的斜率； R e n d及 R s t a r t为起终点的光谱反射率； λ e n d及 λ s t a r t为起终点的波长值； R H i为 i点对应的“HULL”值， C R i即为 i点的包络线值。包络线值的大小取决于端元光谱中中间波段相对于左右波段的凸出或凹陷程度，波峰越大，包络线值越大，波谷越深，则该点的包络线值越小，由图5可知，PV与NPV的包络线值随着波段的变化有着类似的变化趋势，且PV曲线的起伏程度更加剧烈，在GREEN波段， PV的包络线值明显大于另外2个端元，这是因为PV在绿光波段有一个明显的反射峰，这个特征是NPV与BS所不具备的；而在SWIR1波段， NPV则明显突出。包络线值方法的引入，可以放大不同端元光谱曲线间容易被忽视的微小波峰与波谷间的差异，是将NPV信息从混合光谱中分离出来的关键指标之一。
图5 PV、NPV、BS端元光谱对应Sentinel-2波段的包络线曲线

Fig. 5 The envelope curves of PV, NPV, and BS endmember spectra corresponding to Sentinel-2 bands

3.2 稀疏NPV指数模型构建

采用回归模型方法对PV-NPV-BS混合像元进行分解需要指数的支持，因此就要保证指数能够在混合端元中提取出NPV信息[26-27]，而NPV与PV、NPV与BS区分的难度均较大，为使构建的指数能够较好的将NPV进行分离，将指数分为2个部分，第一部分对NPV与PV进行区分，第二部分则是区分NPV与BS， 2个部分的乘积作为最终的非光合植被指数。经过上文对NPV端元光谱的特征分析，总结4种方法处理之后的光谱数据发现：在可见光波段， PV与NPV具有相近的光谱反射特征， BS反射率较高，但PV具有明显的“两谷一峰”特征， NPV在GREEN波段处（560 nm）的包络线值最大（图5），所以选取此处的包络线值用于扩大PV与NPV的差异，并以此作为指数的第一部分；此外，在短波红外波段， BS相较于PV和NPV曲线变化趋势较为平缓， PV与NPV变化趋势则几乎一致，因此将SWIR1、SWIR2波段取比值作为指数的第二部分。
研究中发现，指数要实现对低覆盖地区的NPV信息提取， NPV端元的指数值应当在BS和PV中处于优势地位，因为回归模型是依据指数值的大小来对fNPV进行定量，指数值越大， fNPV相对越高，指数值越小， fNPV相对越小。此外，研究区域的 fNPV很低， NPV对混合像元指数值的贡献很小（式（6））， NPV指数值优势不明显会导致NPV信息被背景掩埋，所以较大的指数值是将NPV与PV、BS区分开来的关键因素，并且有利于利用回归经验模型来对fNPV进行反演。基于上述条件以及对端元光谱数据的分析，构建了包络线差值植被指数（Envelope Difference Vegetation Index, EDVI）（式（7）—式（9））。
I = ∑ f i N i = f P V I P V + f N P V I N P V + f B S I B S
式中： I表示混合像元的指数值； I P V、 I N P V、 I B S，分别表示PV、NPV、BS的指数值； f P V、 f N P V、 f B S，分别表示像元内PV、NPV、BS的占比。
E D V I = a × b
a = b 11 / b 12
b = 1 / b 3 - 0.614 b 2 - 0.386 b 4 + 0.01
式中：b2、b3、b4、b11、b12分别为Sentinel-2 MSI影像数据BLUE、GREEN、RED、SWIR1以及SWIR2波段的光谱反射率。

3.3 其他常用NPVIs

在现有的常用NPVIs指数中，本文选取了基于TM、MODIS、Sentinel-2 MSI传感器的共8个NPVIs用于估算fNPV并进行精度评价。为保持在指数计算时的数据一致性，将基于MODIS和Landsat数据构建的指数采用等效的sentinel-2 MSI传感器波段进行计算。PV、NPV、BS光谱以及TM、MODIS、 Sentinel-2 MSI传感器波段对应关系见图6。
图6 PV、NPV、BS端元光谱曲线及Landsat TM、MODIS、Sentinel-2 MSI传感器波段分布

注:图中红色数字为对应传感器波段代码。

Fig. 6 Endmember spectral curves of PV, NPV, BS and band distribution of Landsat TM, MODIS, Sentinel-2 MSI sensors

（1）基于TM波段构建的NPV指数。常用的基于TM数据构建的非光合植被指数主要包括NDI5、NDI7、NDTI、NDSVI等归一化指数以及STI。如图6所示，Landsat TM的B1、B2、B3、B4、B5、B7波段分别与Sentinel-2 MSI的BLUE、GREEN、RED、NIR、SWIR1、SWIR2波段相对应，对波段进行等效替换后的公式见表2。
表2 常用多光谱NPVIs

Tab. 2 Commonly used multispectral NPVIs

NPVIs 公式 公式编号 参考文献
NDI5 NDI5=(b8-b11)/(b8+b11) （10） [14]
NDI7 NDI7=(b8-b12)/(b8+b12) （11） [14]
NDTI NDTI=(b11-b12)/(b11+b12) （12） [15]
NDSVI NDSVI=(b11-b4)/(b11+b4) （13） [16]
STI STI=b11/b12 （14） [15]
DFI DFI=100×(1-b12/b11)×b4/b8A （15） [18]
SWIR32 SWIR32=b12/b11 （16） [17]
S_NDTI4RE (γ=0) S_NDTI4RE (γ=0)=(b8-b7)×2/(b8+b7+1) （17） [19]

注：式中b4、b8、b8A、b11、b12分别为Sentinel-2 MSI影像数据RED、NIR、RE4、SWIR1以及SWIR2波段的光谱反射率。

（2）基于MODIS波段构建的NPV指数。近年来提出的DFI和SWIR32在fNPV估算方面取得了不错的效果。如图6所示，MODIS的B1、B2、B3、B4、B6、B7波段分别与Sentinel-2 MSI的RED、RE4、BLUE、GREEN、SWIR1、SWIR2波段相对应，对波段进行等效替换后的公式见表2。
（3）基于Sentinel波段构建的NPV指数。在NDTI、S_NDTI、STI指数的基础上加入Sentinel红边波段构建的4波段红边指数NDTI4RE、S_NDTI4RE、STI4RE，如式（18）—式（20）。
N D T I 4 R E = γ × b 11 - b 12 b 11 + b 12 + 1 - γ × b 8 - b 7 b 8 + b 7
S _ N D T I 4 R E = γ × b 11 - b 12 ∙ 2 b 11 + b 12 + 1 + 1 - γ × b 8 - b 7 ∙ 2 b 8 + b 7 + 1
S T I 4 R E = γ × b 11 b 12 + 1 - γ × b 8 b 7
式中：b7、b8、b11、b12分别为Sentinel-2 MSI影像数据RE3、NIR、SWIR1以及SWIR2波段的光谱反射率， γ为权重系数。

3.4 精度评估模型

本文构建指数模型精度检验，通过构建模型估算fNPV与地面实测值进行相关性分析，以决定系数(R2，式（21）)、均方根误差（Root Mean Square Error, RMSE，式（22））以及t检验显著性p值作为模型精度评价的依据。相关系数主要反映估算结果与地面实测端元丰度的关系密切程度，相关系数越高，关系越密切。RMSE主要反映模型估算结果与地面实测数据的接近性，值越小模型估算精度越高。 t检测先行计算统计量t值（式（24）），并在t分布表中查找对应的p值， p值反映的是模型估计值与地面实测端元丰度不存在线性关联的概率，通常阈值为5%，即p<0.05时，认为估计值与地面实测端元丰度存在显著的线性关联。
R 2 = ∑ i = 1 n X i - X - Y i - Y - 2 ∑ i = 1 n X i - X - 2 ∑ i = 1 n Y i - Y - 2
R M S E = ∑ i = 1 n X i - Y i 2 / n
s b ˆ = ∑ i = 1 n X i - X - 2 - b ˆ 2 ∑ i = 1 n Y i - Y - 2 n - 2
t = b ˆ s b ˆ
式中： X i、 Y i分别表示估算值和实测值； X -、 Y -分别为估算值和实测值的均值； n为验证样本数量； b ˆ为回归系数估计值。

4 结果与分析

对于3种4波段红边指数NDTI4RE、S_NDTI4RE、STI4RE，由于权重值γ的引入，采用线性回归方法筛选出与研究区内50个验证样地的实测fNPV相关性最好的指数，γ在[0,1]区间取值，步长为0.1， 3种指数的相关性系数（R）变化如图7所示。3个指数的相关性指数均在权重值γ为0时达到最大值，其中S_NDTI4RE指数相关性系数最大，达到了0.76，因此，本文仅选取S_NDTI4RE(γ=0)指数参与后续的比较与分析。
图7 NDTI4RE、S_NDTI4RE、STI4RE指数相关性系数与权重值γ关系

Fig. 7 Relationship between NDTI4RE, S_NDTI4RE, STI4RE index correlation index and weight value γ

4.1 基于影像的NPV指数特征空间

基于研究区的NDVI以及NPVIs指数值，绘制NPVIs-NDVI特征空间（图8），PV端元的NDVI值较高而NPV与BS端元的值处于一个较低的水平，NPV端元的NPVIs值较高而PV与BS的值较低（SWIR32反之），因此，像素点应在由NPVIs与NDVI组成的二维空间中应呈三角形分布。在图8的散点图中，点的颜色越紫代表散点分布越稀疏，颜色越红代表点的分布越密集，在所有的散点图中，绝大多数散点均密集分布在横坐标的左侧，说明研究区大多数像素的NDVI值偏低，这与干旱与半干旱地区植被稀疏，BS占据大部分面积的特征相匹配。其中，经EDVI（图8(a)）、DFI（图8(b)）、NDTI（图8(f)）处理的像素点的三角形分布特征较为明显。特别是在低NDVI区间（0~0.2），EDVI-NDVI的特征空间表现为明显的三角形，由于研究区植被覆盖度较低，本文主要关注在NDVI较低时的散点分布特征。
图8 NPVIs与NDVI之间特征空间

Fig. 8 Feature Space of NPVIs and NDVI

4.2 线性回归与精度评价

将EDVI与表2中的8种指数（DFI、 SWIR32、STI、 NDTI、 NDI7、 NDI5、 NDSVI、 S_NDTI4RE (γ=0)）结合50个建模样本的实测fNPV进行线性回归分析，计算其R2、RMSE以及显著性p值，评估9个模型的估算精度，并对其在干旱区反演fNPV的可行性进行评价（表3）。
表3 NPVIs回归模型精度

Tab. 3 NPVIs regression model precision

NPVIs RMSE R2 p
EDVI 0.017 0 0.71 <0.01
S_NDTI4RE ( γ = 0) 0.018 5 0.57 <0.01
DFI 0.017 4 0.30 <0.01
SWIR32 0.017 4 0.30 <0.01
NDTI 0.017 3 0.29 <0.01
STI 0.016 9 0.27 <0.01
NDI7 0.014 5 0.17 <0.01
NDI5 0.011 0 0.07 <0.05
NDSVI 0.008 2 0.03 <0.05
如表3所示，9个模型中，EDVI回归模型的指数值与地面实测fNPV的相关性最好， R2达到了0.71，RMSE=0.0170， p<0.01，通过了显著性F检验，达到了极显著性相关水平，且散点均匀分布在回归线两侧（图9(a)），与其余的模型性比，估算精度有了大幅度的提升。除EDVI外，S_NDTI4RE (γ=0)也表现出了较好的估算精度，对除EDVI外的剩余模型有着明显的优势，相关性也明显好于剩余模型，R2达到了0.57，但S_NDTI4RE(γ=0)对部分fNPV出现了较为明显的低估，这也是导致其相关性不如EDVI的重要原因。DFI、STI、NDTI、SWIR32、NDI7、NDI5和NDSVI则在精度评定中表现出较差的相关性，受BS背景效应和PV的影响较大，无法在PV、NPV和BS共3种成分混合情况下，对NPV进行有效分离，达不到用于干旱地区fNPV估算的要求。
图9 不同非光合植被指数指数值与实测值散点图

Fig. 9 The scatter diagram of different non-photosynthetic vegetation index values and measured values

5 讨论

5.1 NPV光谱分析及与最优指数的构建

为解决干旱与半干旱地区PV-NPV-BS混合场景下的NPV信息提取问题，本文发展了一种基于Sentinel-2 MSI 数据的非光合植被指数EDVI。EDVI是一种五波段多光谱指数，它分别考虑了可见光和短波红外波段中NPV与PV和BS的光谱差异，分别在多个NPV特征波段对PV和BS进行区分，采用包络线差值方法、多波段比值方法对光谱差异进行放大，实现稀疏NPV覆盖度精准估算。
本文采集了多种类型NPV、PV、BS的端元光谱数据，分析发现，同类物种的光谱差异主要集中在光谱反射率的大小上，曲线的起伏形态和变化趋势大体是相似的，所以本文在构建指数过程中，没有直接选择波段反射率参与计算，而是利用一阶导数、倒数的对数和连续统去除法对光谱数据进行先行处理，其中连续统去除法的加入显著提升了指数对于曲线间微小差异的提取能力[25]，在一定程度上减弱了不同物种的光谱差异对指数的影响。通过对研究区平均光谱的分析发现：NPV在短波红外波段与PV有着较高的相似性，差异主要集中在可见光和近红外波段，NPV在蓝光和红光波段具有优势而PV在红光到近红外波段的光谱曲线具有更高的斜率（图3（d）），这与Cao等[18]观察到的NPV光谱特征相似； PV在Sentinel-2数据的RE波段与可见光波段之间具有最大的光谱差异（图4），可有效将植被与土壤进行分离，这与Liu等[19]的结果一致，但NPV在RE波段与可见光波段之间同样具有较高的光谱差异，无法在PV背景中对NPV进行有效分离，故未将其加入指数的运算； BS在可见光波段的反射率更高且在短波红外波段与植被有明显的差异，同样适用于土壤与植被的分离。

5.2 不同NPVIs指数提取NPV特征信息差异

本文将新构建EDVI与常用NPVIs与实测fNPV通过线性回归进行分析，绘制散点图（图9），其中，EDVI的指数值与实测值散点图中，散点分布均匀，无明显的fNPV高估和低估问题，且无异常值存在，散点与拟合线较为接近，在精度评定中表现出较高的相关性。其余NPVIs模型与fNPV相关性从大到小依次为S_NDTI4RE (γ=0)、DFI、SWIR32、NDTI、STI、NDI7、NDI5、NDSVI，其中，S_NDTI4RE (γ=0)、DFI、SWIR32、NDTI、STI和NDI7达到了极显著相关水平。此外，同样基于PV、NPV和BS共3种成分混合情况下构建的S_NDTI4RE(γ=0)，也能够较为准确的估算PV、NPV和BS混合情况下干旱与半干旱地区的fNPV， Sentinel-2数据中特有的RE波段的加入对 fNPV的估算起到了促进作用，NDTI4RE、S_NDTI4RE、STI4RE的相关性均明显优于NDTI、S_NDTI、STI[19]。与其他指数相比，EDVI、S_NDTI4RE(γ=0)和DFI加入了更多的光谱波段参与运算，在一定程度上减少了背景效应的影响，提高了fNPV估算的准确度[18]。SWIR32和STI以短波红外波段的光谱差异为基础，能够对NPV和BS进行有效分析，但容易受到PV的影响，导致NPV的特征被掩盖[28-29]。NDTI、NDI7、NDI5和NDSVI分别选取了不同的波段进行归一化，缺乏对PV-NPV-BS混合光谱的综合考量，回归精度较差，此发现与Daughtry等 [30]比较NDTI、NDI7、NDI5和NDSVI估算作物残茬覆盖度的估算精度结果一致。

5.3 EDVI指数在不同植被景观类型应用可靠性有待研究

虽然新构建EDVI相比于现有的NPVIs在干旱地区以灌丛为主的植被景观区域估算fNPV更有优势，但针对不同区域和不同植被景观结构类型（如草地和林地）的fNPV反演模型构建上还有待于进一步研究。前期学者的研究中也主要针对单一植被景观结构群落开展实验，如王光镇等[31]基于NDVI-DFI构建的三元线性混合模型在锡林郭勒草原的 fNPV估算；郑国雄等[10]将指数与光谱混合模型结合，构建了中国北方沙地中草本fNPV的反演模型， Liu等[19]以西藏高山草甸为研究对象构建NPV指数模型； Ji[32]以民勤梭梭和白刺植被景观为研究对象，构建线性和非线性光谱混合模型估算NPV覆盖度。以上研究中，研究区域有限，植被结构和景观类型单一，本文存在同样问题，缺少在乔木林、草地、农作物等其他植被景观类型的可靠性验证。此外，由于本文采用平均端元光谱作为单一纯净端元光谱理论真值，忽略了端元光谱可变性，但光谱线性或非线性混合以及端元可变性造成误差依旧不可避免。在大尺度区域包括多种端元及植被类型景观情况下，端元类型越丰富，该问题会愈加突出，多端元光谱混合分析（Multiple Endmember Spectral Mixture Analysis, MESMA）[33]和自动蒙特卡洛解混法（Automated Monte Carlo Unmixing, AutoMCU）[34]为此问题提供了新的解决方案，将指数作为因子加入到解混算法中是否会有更好的效果有待于进一步的研究。

6 结论

本文选择甘肃民勤县稀疏植被覆盖区为研究区，非光合植被为研究对象，地面采集端元高光谱数据和Sentinel-2 MSI遥感数据为数据源，构建基于多光谱遥感影像提取稀疏NPV覆盖信息的指数模型，主要得出以下结论：
（1）为有效排除裸地光谱影响和PV光谱干扰，稀疏NPV指数构建应尽量以NPV的指数值与BS和PV在光谱特征上差异明显为准则，通过光谱分析发现，相较于原始光谱法、一阶导数变换法、倒数的对数法光谱分析，采用连续统去除法获取的NPV端元光谱包络线值与BS和PV端元光谱包络线值产生明显差异。
（2）在绿光以及短波红外波段，NPV端元光谱特征最突出，通过绿光波段的包络线差值和短波红外波段的光谱反射率比值计算后，可以有效拉开NPV与PV和BS光谱特征值差距，从而实现提取NPV信息，并以此为基础构建新NPV指数；
（3）通过地面量测非光合植被覆盖度验证得知，新构建NPV覆盖信息估算的EDVI指数模型能够有效应用于干旱与半干旱地区稀疏非光合植被覆盖度估算，与地面真实验证数据相关性较高，R2为0.71；与其他NPVIs比较发现，仅次于EDVI的为S_NDTI4RE(γ=0)，R2为0.57；再者为DFI和SWIR32，R2为0.30；其余指数相关性较差。

References

Publishing order | Descend order by publishing year | Descend order by cited within
[1]
李晓松. 干旱地区稀疏植被覆盖度高光谱遥感定量反演研究[D]. 北京: 中国林业科学研究院, 2008.

[Li X S. Quantitative inversion of hyperspectral remote sensing for sparse vegetation coverage in arid areas[D]. Beijing: Chinese Academy of Forestry Sciences, 2008.]

[2]
Guli J, Xi C, Bao A M. A comparison of methods for estimating fractional vegetation cover in arid regions[J]. Agricultural and Forest Meteorology, 2011, 151(2011):1698-1710. DOI:10.1016/j.agrformet.2011.07.004

DOI

[3]
姬翠翠, 骆义峡, 李晓松, 等. Sentinel-1和Sentinel-2协同反演稀疏非光合植被覆盖度方法研究[J]. 遥感学报, 2021.

[Ji C C, Luo Y X, Li X S, et al. Research on Synergistic Inversion of Sparse Non-Photosynthetic Vegetation Coverage by Sentinel-1 and Sentinel-2[J]. Journal of Remote Sensing, 2021.] DOI:10.11834/jrs.20211207

DOI

[4]
Daughtry C S T, Hunt E R, Mcmurtrey JE. Assessing crop residue cover using shortwave infrared reflectance[J]. Remote sensing of environment, 2004, 90(2004):126-134. DOI:10.1016/j.rse.2003.10.023

DOI

[5]
Okin G S. Relative spectral mixture analysis — a multitemporal index of total vegetation cover[J]. Remote Sensing of Environment, 2007, 106(4),467-479. DOI:10.1016/j.rse.2006.09.018

DOI

[6]
Arneth A. Uncertain future for vegetation cover[J]. Nature, 2015, 524(7563),44-45. DOI:10.1038/524044a

DOI

[7]
章文波, 符素华, 刘宝元. 目估法测量植被覆盖度的精度分析[J]. 北京师范大学学报(自然科学版),2001, 3(2001):402-408.

[Zhang W B, Fu S H, Liu B Y. Accuracy analysis of vegetation coverage by visual estimation method[J]. Journal of Beijing Normal University (Natural Science Edition), 2001, 3(2001):402-408.]

[8]
柴国奇, 王静璞, 王光镇. 基于MODIS数据的典型草原非光合植被覆盖度估算[J]. 国土资源遥感, 2019, 31(3):234-241.

[Chai G Q, Wang J P, Wang G Z, et al. Estimation of fractional cover of non-photosynthetic vegetation in typical steppe based on MODIS data[J]. REMOTE SENSING FOR LAND & RESOURCES, 2019, 31(3):234-241.] DOI:10.6046/gtzyyg.2019.03.29

DOI

[9]
Zheng G X, Bao A M, Li X S, et al. The potential of multispectral vegetation indices feature space for quantitatively estimating the photosynthetic, non-photosynthetic vegetation and bare soil fractions in northern China[J]. Photogrammetric Engineering & Remote Sensing, 2019, 85(1),65-76. DOI:10.14358/pers.85.1.65

DOI

[10]
姬翠翠. 干旱区稀疏光合与非光合植被多尺度光谱混合解析方法研究[D]. 武汉: 武汉大学. 2018.

[Ji C C. Research on muti-scale spectral mixture analysis method for sparse photosynthetic/non-photosynthetic vegetation in arid area[D]. Wuhan: Wuhan University. 2018.]

[11]
Guo Z C, Wang T, Liu S L, et al. Comparison of the backpropagation network and the random forest algorithm based on sampling distribution effects consideration for estimating non-photosynthetic vegetation cover[J]. International Journal of Applied Earth Observation and Geoinformation, 2021, 104:1569-8432. DOI:10.1016/j.jag.2021.102573

DOI

[12]
骆义峡, 姬翠翠, 李晓松. 植被指数模型估算干旱区稀疏光合/非光合植被覆盖度[J]. 遥感信息, 2022, 37(03):57-64.

[Luo Y X, Ji C C, Li X S, et al. Vegetation Index Model for Estimating Sparse Photosynthetic/Non- photosynthetic Vegetation Fractional Cover in Arid Zone[J]. Remote Sensing Information, 2022, 37(03):57-64.] DOI:10.3969/j.issn.1000-3177.2022.03.009

DOI

[13]
Mcnairn H and Protz R. Mapping corn residue cover on agricultural fields in Oxford county, Ontario, using thematic mapper[J]. Canadian Journal of Remote Sensing, 1993, 19(2):152-159. DOI:10.1080/07038992.1993.10874543

DOI

[14]
van Deventer A P, Ward A D, Gowda P H, et al. Using thematic mapper data to identify contrasting soil plains and tillage practices[J]. Photogrammetric Engineering and Remote Sensing, 1997, 63(1):87-93. DOI:0099-1112/97/6301-087$3.00/0

[15]
Qi J, Marsett R, Heilman P, et al. Ranges improves satellite - based information and land cover assessments in southwest United States[J]. EOS Transactions American Geophysical Union, 2002, 83(51):601-606. DOI:10.1029/2002eo000411

DOI

[16]
Guerschman J P, Hill M J, Renzullo L J, et al. Estimating fractional cover of photosynthetic vegetation, non-photosynthetic vegetation and bare soil in the Australian tropical savanna region upscaling the EO-1 Hyperion and MODIS sensors[J]. Remote Sensing of Environment, 2009, 113(5):928-945. DOI:10.1016/j.rse.2009.01.006

DOI

[17]
Cao X, Chen J, Matsushita B, et al. Developing a MODIS-based index to discriminate dead fuel from photosynthetic vegetation and soil background in the Asian steppe area[J]. International Journal of Remote Sensing, 2010, 31(6):1589-1604. DOI:10.1080/01431160903475274

DOI

[18]
Liu J L, Fan J R, Yang C, et al. Novel vegetation indices for estimating photosynthetic and non-photosynthetic fractional vegetation cover from Sentinel data[J]. International Journal of Applied Earth Observation and Geoinformation, 2022,109,102793. DOI:10.1016/j.jag.2022.102793

DOI

[19]
Ji C C, Li X S, Wei H D, et al. Comparison of different multispectral sensors for photosynthetic and non-photosynthetic vegetation-fraction retrieval[J]. Remote Sensing, 2020, 12(1):115-132. DOI:10.3390/rs12010115

DOI

[20]
Dong T, Liu J, Shang J, et al. Assessment of red-edge vegetation indices for crop leaf area index estimation[J]. Remote Sensing of Environment, 2019, 222:133-143. DOI:10.1016/j.rse.2018.12.032

DOI

[21]
Muir J, Schmidt M, Tindall D, et al. Field measurement of fractional ground cover[M]. Canberra: Australian Bureau of Agricultural and Resource Economics and Sciences. 2011:4-4.

[22]
王堃. 贵州威宁草海湿地优势植被群丛的光谱特征分析及其识别[D]. 贵州: 贵州师范大学, 2020.

[Wang K. Spectral characteristics analysis and identification of dominant vegetation clusters in Caohai wetland, Weining[D]. Guizhou: Guizhou. Guizhou Normal University, 2020.] DOI:10.27048/d.cnki.ggzsu.2019.000516

DOI

[23]
孙岩. 湿地植物高光谱特征分析与物种识别模型构建[D]. 北京: 清华大学, 2008.

[Sun Y. Analysis of hyperspectral characteristics of wetland plants and construction of species identification model[D]. Beijing: Tsinghua University, 2008.]

[24]
Alexander F H G. Three decades of hyperspectral remote sensing of the Earth: A personal view[J]. Remote Sensing of Environment, 2009, 113(1):S5-S16. DOI:10.1016/j.rse.2007.12.014

DOI

[25]
潘晓玲. 干旱区绿洲生态系统动态稳定性的初步研究[J]. 第四纪研究, 2001, 21(4):345-351.

[Pan X L. A preliminary study on the dynamic stability of oasis ecosystem in arid areas[J]. Quaternary Research, 2001, 21(4):345-351.]

[26]
张云霞, 李晓兵, 陈云浩. 草地植被盖度的多尺度遥感与实地测量方法综述[J]. 地球科学进展, 2003, 18(1):85-93.

DOI

[Zhang Y X, Li X B, Chen Y H. A review of multi-scale remote sensing and field measurement methods for grassland vegetation cover[J]. Advances in Earth Science, 2003, 18(1):85-93.]

[27]
Gelder B K, Kaleita A L, Cruse R M. Estimating mean field residue cover on midwestern soils using satellite imagery[J]. Agronomy Journal, 2009, 101(3):635-643. DOI:10.2134/agronj2007.0249

DOI

[28]
Jacques D C, Kergoat L, Hiernaux P, et al. Monitoring dry vegetation masses in semi-arid areas with MODIS SWIR bands[J]. Remote Sensing of Environment, 2014, 153(2014):40-49. DOI:10.1016/j.rse.2014.07.027

DOI

[29]
Daughtry C S T, Hunt E R, Doraiswamy P C, et al. Remote sensing the spatial distribution of crop residues[J]. Agronomy Journal, 2005, 97(3):864-871. DOI:10.2134/agronj2003.0291

DOI

[30]
王光镇, 王静璞, 韩柳, 等. 基于实测光谱模拟Landsat-8 OLI数据估算非光合植被覆盖度[J]. 地球信息科学学报, 2018, 20(11):1667-1678.

DOI

[Wang G Z, Wang J P, Hai L, et al. Estimating fractional cover of non-photosynthetic vegetation using field spectral to simulate Landsat-8 OLI[J]. Journal of Geo-information Science, 2018, 20(11):1667-1678.] DOI:10.12082/dqxxkx.2018.180196

DOI

[31]
Ji C C, Jia Y H, Gao Z H, et al. Nonlinear spectral mixture effects for photosynthetic/non-photosynthetic vegetation cover estimates of typical desert vegetation in western China[J]. PLOS ONE, 2017, 12(12):e0189292. DOI: 10.1371/journal.pone.0189292

DOI

[32]
Roberts D A, Gardner M, Church R, et al. Mapping chaparral in the santa monica mountains using multiple endmember spectral mixture models[J]. Remote Sensing of Environment, 1998, 65(3):267-279. DOI:10.1016/s0034-4257(98)00037-6

DOI

[33]
Asner G P, Lobell D B. A biogeophysical approach for automated SWIR unmixing of soils and vegetation[J]. Remote Sensing of Environment, 2000, 74(1):99-112. DOI:10.1016/s0034-4257(00)00126-7

DOI

Options
Outlines

/

〈 〉