CN111402182B - 一种基于土地覆盖信息的中分影像合成方法 - Google Patents

Abstract

本发明一种基于土地覆盖信息的中分影像合成方法。该方法首先对目标区的时序影像进行预处理，将原始影像各波段的像元量化值转换为大气顶层反射率；对影像进行正射校正，保证影像间的配准误差小于1个像元；计算每景影像的归一化植被指数；其次对目标区域划分格网，然后用格网裁切影像得到分块数据；接着对每个分块数据进行NDVI最大值影像合成得到对应分块影像A；然后对每个分块数据进行蓝波段反射率最小值影像合成得到对应分块影像B；最后基于土地覆盖类型合成影像A和B得到分块最终影像C，拼接分块影像，得到覆盖目标区的最终合成影像D。

Description

一种基于土地覆盖信息的中分影像合成方法

技术领域

本发明属于遥感图像处理领域，具体与时间序列中分光学遥感卫星数据的影像合成方法相关。

背景技术

任意一天平均有64％的地球表面被云覆盖，有超过一半的对地观测的光学遥感卫星影像都或多或少受到云的影响。降低云及其阴影对光学影像数据分析的影响是遥感影像处理领域的一个重要问题。提供无云的光学遥感数据对土地覆盖、农业监测等应用都非常重要。影像合成是从一系列纠正好的时间序列影像中按照一定的标准挑选出“最好”的“观测”，合成一幅在该时间段能够代表该区域、在空间上连续且去除云和大气影响的数据集。

AVHRR、MODIS和Landsat数据由于有热红外波段，而云的温度相对地表明显偏低，所以这些数据的影像合成方法大多采用了热红外波段。本文提出了一种有效的中分辨率时间序列影像合成新方法，该方法基于已有的NDVI最大值和蓝波段反射率最小值的影像合成方法，并引入土地覆盖先验数据辅助影像合成，生成无云影像，不需要热红外波段。该方法解决了以往影像合成方法中水体区域效果不好的问题，同时最大程度保留了陆地区域的有效信息。

发明内容

本发明解决的技术问题是：克服现有技术的不足，本发明提出了一种基于土地覆盖信息的中分影像合成方法，针对无热红外波段的中分时序数据合成无云影像，解决了水体区域合成效果不佳的问题，最终合成影像上，陆地和水体区域的云和阴影均得到了有效去除，效果较好。

本发明的技术方案是：一种基于土地覆盖信息的中分影像合成方法，包括以下步骤：

步骤1、时间序列影像预处理；

对固定时间段(比如30天)内时间序列影像做3步预处理工作，包括大气层顶反射率计算、正射校正和NDVI计算，具体如下：

(1)大气层顶反射率计算

使用辐射定标系数把影像各个波段的DN值转为辐亮度，具体公式如下：

L＝Gain*DN+Bias (1)

其中L是卫星载荷波段入瞳处等效辐射亮度(W·m^-2·sr^-1·μm^-1)，Gain和Bias分别为定标系数增益、偏移量。

计算各波段的大气顶层反射率B，具体公式如下：

其中B为大气顶层反射率(无量纲)，i为第i波段，π为常量，L为辐亮度，D为日地距离(天文单位)，ESUN为大气层顶的平均太阳光谱辐照度(W·m^-2·sr^-1·μm^-1)，θ为太阳的天顶角。

(2)正射校正

对时间序列影像进行正射校正，消除原始影像的几何形变，并保证时间序列影像之间的几何偏差要小于一个像元。

(3)NDVI计算

单景影像的NDVI计算需要红波段和近红外波段的反射率，具体公式如下：

NDVI＝(NIR–R)/(NIR+R) (3)

其中R和NIR分别为红波段和近红外波段的大气层顶反射率。

步骤2、对目标区域划分格网，然后用格网裁切覆盖目标区域的时间序列影像得到分块数据；

为了降低计算压力，减少重复运算，采用分块的方式进行影像合成。预先划定方形的矢量格网覆盖目标区域，每块大小可根据计算能力确定，例如每块3000*3000个像元(900万个像元)；用矢量格网去裁切经过预处理的单景影像和NDVI影像，得到每个分块的时间序列影像。

步骤3、对分块数据进行NDVI最大值影像合成得到对应分块合成影像A；

对于每个像元所在位置，有若干个时相的数据对应，选取其中NDVI最大值所对应的时相的影像像元作为合成后影像像元值，得到合成影像A。

步骤4、对分块数据进行蓝波段反射率最小值影像合成得到对应分块合成影像B；

类似于步骤3，但像元的选取标准为蓝波段大气层顶反射率最小值，进而得到合成影像B。

步骤5、基于土地覆盖类型合成影像A和B，得到分块最终影像C，拼接分块影像，得到最终的合成影像。

基于目标区的土地覆盖先验数据将目标区分为植被和非植被区，植被区包括农田、森林、草地、灌木等，非植被区包括例如水体、不透水地表、裸地、雪/冰等。对于每个像元所在位置，如果土地覆盖先验数据为植被区，选取影像A对应的像元作为合成后影像的像元值；反之，如果土地覆盖先验数据为非植被区，选取影像B对应的像元作为合成后影像的像元值；影像A和B合成得到分块最终影像C；拼接分块影像得到最终的合成影像。

本发明相比现有技术的优点：

(1)建立了一种基于土地覆盖信息的中分影像合成方法，该方法对于无热红外波段的中分时序数据合成无云影像非常有效。

(2)本文提出的方法引入先验的土地覆盖信息，结合了NDVI最大值合成和蓝波段反射率最小值合成算法各自的优势，该方法简单可靠。

附图说明

图1为分块(编号7)时间序列影像和不同方法的合成影像结果，其中图1a为分块7成像时间2014年7月22日影像，红波段灰度显示；图1b为分块7成像时间2014年7月27日影像，红波段灰度显示；图1c为分块7成像时间2014年7月30日影像，红波段灰度显示；图1d为分块7成像时间2014年7月31日影像，红波段灰度显示；图1e为分块7成像时间2014年7月31日影像(不同相机)，红波段灰度显示；图1f为土地覆盖图，植被区黑色，非植被区白色；图1g为NDVI最大值合成结果影像，红波段灰度显示；图1h为蓝波段最小值合成结果影像，红波段灰度显示；图1m为最终合成结果影像，红波段灰度显示。

图2为覆盖湖南省岳阳市行政辖区的时间序列单景正射影像叠合图和最终的合成结果，其中图2a为覆盖岳阳市的时间序列单景正影像叠合图，红波段灰度显示，岳阳市行政边界显示为白色曲线，分块的边界线显示为白色直线，每个分块的编号显色为白色的数字；图2b为覆盖岳阳市的最终合成结果影像，红波段灰度显示，岳阳市行政边界显示为白色曲线。

具体实施方式

一种基于土地覆盖信息的中分影像合成方法，以覆盖湖南省岳阳市行政辖区的2014年7月的10景GF-1WFV影像(见图2a)合成为例，具体工作步骤如下：

步骤1、时间序列影像预处理；

对GF-1WFV时间序列影像做3步预处理工作，包括大气层顶反射率计算、正射校正和NDVI计算，具体如下：

(1)大气层顶反射率计算

GF-1WFV影像的辐射定标系数可在中国资源卫星应用中心的官方网站下载。使用辐射定标系数把影像各个波段的DN值转为辐亮度，具体公式如下：

L＝Gain*DN+Bias (1)

其中L是卫星载荷波段入瞳处等效辐射亮度(W·m^-2·sr^-1·μm^-1)，Gain和Bias分别为定标系数增益、偏移量。

计算GF-1WFV的4个波段的大气顶层反射率B，具体公式如下：

其中B为大气顶层反射率(无量纲)，i为第i波段，π为常量，L为辐亮度，D为日地距离(天文单位)，ESUN为大气层顶的平均太阳光谱辐照度(W·m^-2·sr^-1·μm^-1)，θ为太阳的天顶角。

(2)正射校正

对10景GF-1WFV时间序列影像进行正射校正，消除原始影像的几何形变，并保证时间序列影像之间的几何偏差要小于一个像元。

(3)NDVI计算

单景影像的NDVI计算需要红波段和近红外波段的反射率，具体公式如下：

NDVI＝(NIR–R)/(NIR+R) (3)

其中R和NIR分别为GF-1WFV第3波段红波段和第4波段近红外波段的大气层顶反射率。

步骤2、对目标区域划分格网，然后用格网裁切覆盖目标区域的时间序列影像得到分块数据；

为了降低计算压力，减少重复运算，采用分块的方式进行影像合成；预先划定方形的矢量格网覆盖岳阳市行政辖区，每块3000*3000个像元(900万个像元)，共12块，见图2a；用矢量格网去裁切经过预处理的单景影像和NDVI影像，得到每个分块的时间序列影像和对应的NDVI影像，例如编号7的分块就有5个时间序列影像，见图1a、图1b、图1c、图1d和图1e。

步骤3、对分块数据进行NDVI最大值影像合成得到对应分块合成影像A；

每个像元所在位置，有若干个时相的数据与之对应，选取其中NDVI最大值所对应的时相的影像像元作为合成后影像像元值，得到合成影像A，见图1g。

步骤4、对分块数据进行蓝波段反射率最小值影像合成得到对应分块合成影像B；

类似于步骤3，但像元的选取标准为蓝波段大气层顶反射率最小值，得到合成影像B，见图1h。

步骤5、基于土地覆盖类型合成影像A和B得到分块最终影像C，拼接分块影像，得到最终的合成影像。

基于覆盖岳阳市的土地覆盖先验数据将每一个分块数据所覆盖区域分为植被和非植被区。图1f为编号为7的分块的土地覆盖展示图，植被区显示为黑色，非植被区显示为白色。对于每个像元所在位置，如果土地覆盖先验数据为植被区，选取影像A对应的像元作为合成后影像的像元值；反之，如果土地覆盖先验数据为非植被区，选取影像B对应的像元作为合成后影像的像元值；影像A和B合成得到分块最终影像C，见图1m；拼接分块影像得到最终的合成影像,见图2b。

Claims (6)

1.一种基于土地覆盖信息的中分影像合成方法，其特征在于包括以下步骤：

步骤1、时间序列影像预处理；

步骤2、对目标区域划分格网，然后用格网裁切覆盖目标区域的时间序列影像得到分块数据；

步骤3、对分块数据进行NDVI最大值影像合成得到对应分块合成影像A；

步骤4、对分块数据进行蓝波段反射率最小值影像合成得到对应分块合成影像B；

步骤5、基于土地覆盖类型合成影像A和B，得到分块最终影像C，拼接分块影像，得到最终的合成影像D，具体如下：

基于目标区的土地覆盖先验数据将目标区分为植被区和非植被区，所述植被区包括农田、森林、草地、灌木，非植被区包括水体、不透水地表、裸地、雪、冰；对于每个像元所在位置，如果土地覆盖先验数据为植被区，选取影像A对应的像元作为合成后影像的像元值；反之，如果土地覆盖先验数据为非植被区，选取影像B对应的像元作为合成后影像的像元值；影像A和B合成得到分块最终影像C；拼接分块影像得到最终的合成影像。

2.根据权利要求1所述的一种基于土地覆盖信息的中分影像合成方法，其特征在于：所述步骤1的具体过程为：对固定时间段内时间序列影像做三步预处理工作，包括大气层顶反射率计算、正射校正和NDVI计算。

3.根据权利要求2所述的一种基于土地覆盖信息的中分影像合成方法，其特征在于：所述大气层顶反射率计算具体为：

使用辐射定标系数把影像各个波段的DN值转为辐亮度，具体公式如下：

L＝Gain*DN+Bias

其中L是卫星载荷波段入瞳处等效辐射亮度，Gain和Bias分别为定标系数增益、偏移量；

计算各波段的大气顶层反射率B，具体公式如下：

其中B为大气顶层反射率，i为第i波段，π为常量，L为辐亮度，D为日地距离，ESUN为大气层顶的平均太阳光谱辐照度，θ为太阳的天顶角。

4.根据权利要求2所述的一种基于土地覆盖信息的中分影像合成方法，其特征在于：所述正射校正具体为：对时间序列影像进行正射校正，消除原始影像的几何形变，并保证时间序列影像之间的几何偏差要小于一个像元。

5.根据权利要求2所述的一种基于土地覆盖信息的中分影像合成方法，其特征在于：所述NDVI计算具体为：

单景影像的NDVI计算需要红波段和近红外波段的反射率，具体公式如下：

NDVI＝(NIR–R)/(NIR+R)

其中R和NIR分别为红波段和近红外波段的大气层顶反射率。

6.根据权利要求1所述的一种基于土地覆盖信息的中分影像合成方法，其特征在于：所述步骤2的具体过程为：采用分块的方式进行影像合成；预先划定方形的矢量格网覆盖目标区域，每块大小根据计算能力确定。

Images