CN111126203A

CN111126203A - 基于ndvi百分比匹配的浓密植被识别方法

Info

Publication number: CN111126203A
Application number: CN201911272623.3A
Authority: CN
Inventors: 孙林; 于会泳; 李睿博
Original assignee: Shandong University of Science and Technology
Current assignee: Shandong University of Science and Technology
Priority date: 2019-12-04
Filing date: 2019-12-04
Publication date: 2020-05-08

Abstract

本发明提供了一种基于NDVI百分比匹配的浓密植被识别方法，主要针对遥感卫星缺少短波红外波段，浓密植被识别困难。在假设浓密植被的百分比在时间上具有相对稳定的延续性，在同一区域相邻的几年里浓密植被所占的比重年际变化不大以及即使在观测角度、气溶胶等条件的影响下植被的NDVI值依然高于其它地表类型在相同条件下的NDVI值的两个前提下，首先构建MODIS地表反射率数据集，按最小合成法将长时间序列内同一日期的多景影像合成一景影像；然后对遥感卫星多光谱数据进行辐射定标、均值重采样、NDVI计算等处理；利用MODIS数据统计出遥感卫星多光谱数据对应时空的浓密植被像元百分比，最后根据植被指数NDVI阈值与浓密植被百分比确定遥感卫星数据中的浓密植被像元分布。

Description

基于NDVI百分比匹配的浓密植被识别方法

技术领域

本发明涉及一种浓密植被识别的方法，具体是利用NDVI百分比匹配的方法确定浓密植被像元，适用于缺少短波红外波段的遥感卫星影像浓密植被的确定。

背景技术

由于地表类型的复杂性、地物物候期的不同，地表参数处于不断变化过程中，因此，同一年不同季节以及不同年际之间植被覆盖状况是不同的，造成了浓密植被像元确定的复杂性。传统的浓密植被识别方法主要是针对有红外波段的影像，如MODIS影像和Landsat TM影像。浓密植被像元最初是根据可见光-近红外波段的植被指数来确定，其中归一化植被指数NDVI应用最为广泛。浓密植被像元在红光波段的反射率很低，而在近红外波段反射率较高，浓密植被区的NDVI值大于非浓密植被区，因此可以用归一化植被指数NDVI来提取图像中的浓密植被像元。NDVI为近红外波段的反射率值与红光波段的反射率值之差比上两者之和，即：NDVI＝(ρ_nir-ρ_r)/(ρ_nir+ρ_r)，其中ρ_nir为近红外波段的反射率值，ρ_r为红光波段的反射率值，NDVI取值范围为[-1，1]。经比值处理，能够部分消除与太阳高度角、地形、云、大气条件有关的辐照度条件变化等的影响。

Richeter等(2006)利用比值植被指数、红外波段和近红外波段反射率阈值识别TM影像的浓密植被像元。赵志强等(2015)采用归一化植被指数NDVI、归一化水体指数NDWI以及近红外波段表观反射率等掩膜处理非植被像元，利用浓密森林像元在可见光波段反射率低的特点，通过搜索红波段直方图的最小峰值自动识别山区浓密植被像元，并将该方法应用于环境减灾卫星(HJ-1)CCD影像浓密植被像元识别和气溶胶反演。

NDVI能够较好地反映地表植被分布状况，但是红波段信号受大气中气溶胶的散射和吸收的影响，使得NDVI本身随着气溶胶光学厚度的增加明显减小，当气溶胶光学厚度较大时，浓密植被的植被指数会降到很低，无法通过经验的阈值找到浓密植被像元，也很难用一个统一的阈值将浓密植被识别出来，因此直接利用NDVI方法确定浓密植被像元存在缺陷。

由于短波红外波段波长远远大于大多数气溶胶粒子直径，使得其对气溶胶散射不敏感，但其对地表差异非常敏感，因此，可以用短波红外波段反射率阈值确定浓密植被像元，用短波红外波段的归一化植被指数NDVI_SWIR取代NDVI，即：NDVI_SWIR＝(ρ_1.24-ρ_2.12)/(ρ_1.24+ρ_2.12)其中，ρ_1.24为波长1.24微米处的反射率值，ρ_2.12为波长2.12微米处的反射率值。由于短波红外波段受气溶胶影响相对较小，相对于归一化植被指数NDVI而言，短波红外波段的归一化植被指数NDVI_SWIR提高了浓密植被像元识别的精度。

Kaufman等(1997)提出利用短波红外波段识别浓密植被像元，并成功应用于Landsat TM影像(2.2μm波段)和MODIS影像(2.1μm波段)。然而由于Quickbird影像、IKONOS影像以及国产高分卫星影像如高分系列、HJ-1CCD影像、CBERS影像等仅设置了可见光-近红外波段，缺少短波红外波段，使得浓密植被像元的确定变得异常困难。

发明内容

针对遥感卫星缺少短波红外波段以至浓密植被像元识别困难的问题，本发明提出了一种基于NDVI百分比匹配的方法用于浓密植被的识别，利用现有的MODIS地表反射率产品协助确定遥感卫星多光谱数据的浓密植被像元分布。

本发明的方法有两点基本假设：一是假设浓密植被的百分比在时间上具有相对稳定的延续性，在同一区域相邻的几年里浓密植被所占的比重年际变化不大；二是假设即使在观测角度、气溶胶等条件的影响下植被的NDVI值依然高于其它地表类型在相同条件下的NDVI值。本发明通过遥感卫星数据与MODIS地表反射率产品配准，利用MODIS地表反射率产品数据统计出遥感卫星多光谱数据对应时空的浓密植被像元百分比，根据归一化植被指数NDVI与浓密植被百分比确定遥感卫星数据的浓密植被像元。

本发明采用的技术方案如下：

步骤一：构建地表反射率数据集

1)选取MOD09A1地表反射率数据产品。

2)利用MRT对下载的MOD09A1数据进行投影转换、拼接等预处理。

3)按最小值合成的方式将长时间序列内同一日期的多景影像合成为一景影像。

4)在地表反射率数据库构建时，选取与遥感卫星数据对应的波段及确定浓密植被像元的2.1μm波段，即MODIS蓝光波段、绿光波段、红光波段、近红外波段和2.1μm短波红外波段。

5)利用全国矢量对拼接的MOD09A1地表反射率产品进行裁剪获得全国陆地地表反射率数据集。

步骤二：数据预处理

1)遥感卫星多光谱数据辐射定标。

2)遥感卫星多光谱数据重采样、投影转换几何处理。

3)从MOD09A1地表反射率数据集中裁剪出与遥感卫星多光谱数据时相和空间都匹配的区域。

步骤三：计算归一化植被指数NDVI，并获得遥感卫星多光谱数据的浓密植被百分比。

1)NDVI计算公式如下：

NDVI＝(ρ_nir-ρ_r)/(ρ_nir+ρ_r)

其中，ρ_nir为近红外波段的反射率值，ρ_r为红光波段的反射率值。

2)利用MODIS地表反射率数据，统计出遥感卫星多光谱数据浓密植被的百分比P_ddv。

根据MODIS地表反射率数据确定浓密植被像元的条件为：

统计MODIS地表反射率数据中满足上述两者的像元占总像元个数的百分比P_ddv。

步骤四：根据区域浓密植被的百分比，确定遥感卫星多光谱数据浓密植被像元NDVI的阈值T_ndvi。

浓密植被像元NDVI阈值T_ndvi的确定公式为：

其中，n(NDVI＞T_ndvi)为NDVI值大于阈值T_ndvi的像元个数，N为像元总数。大于NDVI阈值T_ndvi的像元对应的遥感卫星多光谱影像像元为浓密植被像元。

本发明具有如下优点：

本发明针对直接利用NDVI方法确定浓密植被像元的不足以及遥感卫星数据缺少短波红外波段导致浓密植被像元确定困难的问题，提出了使用先验地表反射率产品协助确定遥感卫星多光谱数据浓密植被像元。MODIS地表反射率产品在反映真实地表反射率的变化方面具有较高的精度且2.1μm短波红外波段受大气影响小，可以使用阈值找到浓密植被像元。运用NDVI百分比匹配的方式确定遥感卫星多光谱数据中浓密植被像元的分布，回避了像元的直接匹配，减小异源数据像元配准引起的误差，提高了浓密植被识别的准确度。

附图说明

图1基于NDVI百分比匹配的浓密植被识别的技术流程图；

图2用于浓密植被像元识别的GF-1影像多光谱数据；

图3建立的第97～104天的地表反射率数据集的彩色合成影像；

图4建立的第121～128天的地表反射率数据集彩色合成影像；

图5建立的第145～152天的地表反射率数据集彩色合成影像；

图6建立的第297～304天的地表反射率数据集彩色合成影像；

图7 GF-1多光谱表观反射率数据重采样后的假彩色影像；

图8第121～128天MODIS地表反射率数据的浓密植被分布范围，用绿色矢量线表示；

图9 GF-1影像多光谱数据浓密植被像元的分布范围，用绿色矢量线表示；

具体实施方式

为了使本发明的目的、技术方案更加清楚明白，下面结合具体实施方式对本发明进行详细说明。

本发明选取GF-1多光谱数据作为浓密植被识别的数据源。当GF-1多光谱数据气溶胶光学厚度较大时，浓密植被的植被指数NDVI会降到很低，无法通过经验的阈值找到浓密植被像元。而且GF-1多光谱数据仅设置了可见光近红外波段，没有设置受气溶胶影响较小的短波红外波段，以至其浓密植被像元分布确定十分困难。本发明借助M0D09A1地表反射率产品识别GF-1多光谱影像的浓密植被像元。首先对GF-1多光谱数据进行辐射定标等处理，然后利用MODIS地表反射率产品数据统计出GF-1多光谱数据对应时空的浓密植被像元百分比，最后根据归一化植被指数NDVI与浓密植被百分比确定出GF-1多光谱数据中的浓密植被像元分布。具体步骤如下：

步骤一：GF-1多光谱数据准备

GF-1卫星是目前我国中高分辨率多光谱光学卫星中比较典型的代表。GF-1卫星搭载了两台2m分辨率全色/8m分辨率多光谱相机，四台16m分辨率多光谱宽幅相机。GF-1数据设有4个波段，分别为蓝波段(band1，0.45-0.52μm)、绿波段(band2，0.52-0.59μm)、红波段(band3，0.63-0.69μm)、近红外波段(band4，0.77-0.89μm)。图1是选取的中国山东东部地区的GF-1多光谱卫星影像，卫星观测时间为2013年5月5日，空间分辨率为16米，包括蓝波段、绿波段、红波段、近红外波段4个波段，UTM投影，基准面为WGS-84。

步骤二：构建全国地表反射率数据集

本发明利用MOD09A1产品构建地表反射率数据集。M0D09A1是8天合成的地表反射率产品，在生产过程中考虑了大气和气溶胶的散射和吸收的影响、土地覆盖类型变化、邻近象元效应、地表二向反射特性以及卷云的影响，与真实地表反射率的变化差距较小，空间分辨率为500米，具有7个波段，同一轨道号每年有46景影像。

1)从美国国家航空航天局(NASA)地球观测系统的数据和信息系统(EOSDIS)下载2008～2012五年间覆盖中国大部分陆地区域的所有8天合成的MOD09A1地表反射率数据。其在合成过程中，取临近的8天的数据，尽可能采用质量好的像元替代云和云阴影像元、高视角像元、高气溶胶像元，按最小反射率技术组成一幅影像，最大限度的降低了云的影响以及大气纠正的误差。

2)利用MRT对选取的每一景影像进行投影转换、拼接等预处理。处理后输出影像均采用正弦投影，投影参数设置为第一条标准纬线为25°，第二条标准纬线为47°，中央经线为105°，起始纬度、东偏和北偏均为0°，基准面为WGS-84。影像地面分辨率仍然保持为500m。

3)在地表反射率数据库构建时，选取与遥感卫星数据对应的波段及确定浓密植被像元的2.1μm波段，即MODIS蓝光波段、绿光波段、红光波段、近红外波段和2.1μm短波红外波段。

4)由于薄雾、碎云等因素影响浓密植被识别的精度，为了降低云覆盖的影响，对2008～2012五年内地表参数产品进行时间分析，假设浓密植被的百分比在时间上具有相对稳定的延续性，在同一区域相邻的五年里浓密植被所占的比重年际变化不大，按最小值合成的方式将五年间同一日期的5景影像合成为一景影像。

最小合成法描述如下：以2008～2012五年第121～128天的MOD09A1地表反射率数据为例，将五景影像的相应位置的反射率值按从小到大的顺序排列，每个位置选取最小值参与最后的数据合成，最终形成一景几乎没有云覆盖的地表反射率数据，如图3所示。数据集中包含MOD09A1的红光波段、蓝光波段、绿光波段、近红外波段以及2.1μm短波红外波段。

图2、图4、图5都是按上述方式合成的地表反射率数据，该影像基本全部由晴空像元组成，最终可以获得46景覆盖我国大部分陆地区域的地表反射率数据集。

步骤三：GF-1多光谱数据预处理

1)GF-1多光谱数据辐射定标

首先，建立数字量化值与对应视场中辐射亮度值之间的定量关系，利用式(1)将GF-1卫星多光谱数据的通道观测值DN转换为卫星载荷入瞳处等效表观辐亮度数据。

L＝Gain·DN+Bia (1)

式中，DN为卫星载荷观测值；Gain为定标斜率，Bias为定标截距。GF-1卫星各个载荷的绝对辐射定标系数Gain、Bias可以从GF-1影像元数据中获取。

然后，将表观辐亮度值转换为表观反射率值。

式中，ρ^*为表观反射率，L为式(1)所得的卫星载荷入瞳处等效表观辐亮度值，θ_S为太阳天顶角。D为日地距离订正因子，由式(3)计算得到，E₀为日地平均距离处波段太阳常数，可以由式(4)计算。

日地距离订正因子D的计算公式为：

x＝2π(d_n-1)/365 (4)

d_n为儒略日，如果遇到闰年，则用366代替式(4)中的365。

日地平均距离处波段太阳常数E₀的计算公式为：

式中，L(λ)为太阳光谱辐照度，f(λ)为卫星光谱响应函数。

2)GF-1多光谱数据表观反射率影像重采样

MOD09A1地表反射率数据集的空间分辨率为500米，而经辐射定标的GF-1多光谱数据表观反射率影像的空间分辨率仍为16米，因此需要将GF-1多光谱数据表观反射率影像转换到500米分辨率的尺度，由于采样前后影像的空间分辨率变化较大，本发明中采用均值重采样的方式对GF-1多光谱数据表观反射率影像进行重采样，即对原始表观反射率影像进行n×n像元的合成，将n×n个像元取平均值作为重采样后对应的像元值，其中n＝500/8。图6是GF-1多光谱表观反射率数据重采样后的假彩色影像，空间分辨率为500米；

3)GF-1多光谱数据表观反射率影像投影转换

将正弦投影下的MOD09A1地表反射率数据的投影方式转换为UTM投影，使其投影方式与GF-1多光谱数据表观反射率一致，基准面仍为WGS-84。

4)从MOD09A1地表反射率数据集中裁剪出与GF-1多光谱数据时相和空间都匹配的区域。选择的GF-1多光谱数据的观测时间为2013年5月5日，对应MOD09 A1地表反射率数据集中第121～128天的数据，如图3所示。

步骤四：计算归一化植被指数NDVI，获得GF-1多光谱数据的浓密植被百分比。

1)归一化植被指数NDVI的计算公式如下：

NDVI＝(ρ_nir-ρ_r)/(ρ_nir+ρ_r) (6)

其中，ρ_nir为近红外波段的反射率值，ρ_r为红光波段的反射率值，NDVI取值范围为[-1，1]。

利用式(6)计算步骤三中裁剪的MOD09A1地表反射率数据的NDVI值。

利用式(6)计算步骤三处理后的GF-1多光谱表观反射率影像的NDVI值。

使用MODIS地表反射率数据2.1μm波段确定浓密植被像元，通过对MODIS数据植被覆盖像元的大量统计，确定浓密植被像元的条件为：

统计MODIS地表反射率数据中满足上述两者的像元占总像元个数的百分比P_ddv。在本发明实验中，与GF-1时相和空间都匹配的MODIS区域的浓密植被像元所占的百分比为36.5％，其分布如图7中绿色矢量范围所示。

步骤五：根据区域浓密植被的百分比P_ddv，确定GF-1多光谱数据浓密植被像元NDVI的阈值。

在利用NDVI百分比匹配浓密植被时假设即使在观测角度、气溶胶等条件的影响下植被的NDVI值依然高于其它地表类型在相同条件下的NDVI值。GF-1多光谱数据浓密植被像元NDVI阈值T_ndvi的确定方式为：

其中，n(NDVI＞T_ndvi)为NDVI大于T_ndvi的像元个数，N为像元总数，大于该NDVI阈值T_ndvi的像元对应的GF-1表观反射率影像数据像元为浓密植被像元。根据本发明方法确定的GF-1多光谱数据浓密植被的阈值为：0.665，其分布范围用图8中的绿色矢量表示。

本发明的方法有效解决了缺少短波红外波段的遥感卫星数据浓密植被识别的问题，利用NDVI百分比匹配的方式进行先验地表反射率产品与遥感卫星多光谱数据配准，回避了像元的直接匹配，减小像元配准引起的误差，提高了浓密植被识别的准确度。

Claims

1.一种基于NDVI百分比匹配的浓密植被识别方法，包括以下步骤：

第一步：构建地表反射率数据集

1-1)选取MOD09A1地表反射率数据产品。

1-2)利用MRT对下载的MOD09A1数据进行投影转换、拼接等预处理。

1-3)按最小值合成的方式将长时间序列内同一日期的多景影像合成为一景影像。

1-4)在地表反射率数据库构建时，选取与遥感卫星数据对应的波段及确定浓密植被像元的2.1μm波段，即MODIS蓝光波段、绿光波段、红光波段、近红外波段和2.1μm短波红外波段。

1-5)利用全国矢量对拼接的MOD09A1地表反射率产品进行裁剪获得全国陆地地表反射率数据集。

第二步：数据预处理

2-1)遥感卫星多光谱数据辐射定标。

2-2)遥感卫星多光谱数据重采样、投影转换几何处理。

2-3)从MOD09A1地表反射率数据集中裁剪出与遥感卫星多光谱数据时空相匹配的区域。

第三步：计算归一化植被指数NDVI，并获得遥感卫星多光谱数据的浓密植被百分比。

第四步：根据区域浓密植被的百分比，确定遥感卫星多光谱数据浓密植被像元NDVI的阈值。

2.在第三步中，归一化植被指数NDVI以及浓密植被的百分比P_ddv按如下公式计算：

3-1)NDVI计算公式如下：NDVI＝(ρ_nir-ρ_r)/(ρ_nir+ρ_r)

3-2)根据MODIS地表反射率数据确定浓密植被像元的条件为：

3.在第四步中，遥感卫星多光谱数据浓密植被像元NDVI阈值T_ndvi的确定方式为：

其中，n(NDVI＞T_ndvi)为NDVI大于T_ndvi的像元个数，N为像元总数。大于NDVI阈值T_ndvi的像元对应的原始影像像元为浓密植被像元。