CN107389036B - 一种结合无人机影像的大空间尺度植被覆盖度计算方法 - Google Patents

Abstract

本发明公开了一种结合无人机影像的大空间尺度植被覆盖度计算方法，包括：对遥感影像大气校正和几何纠正，计算NDVI，并根据预定阈值获取有效区域；对无人机图片拼接并获取正射影像，结合几何纠正后的卫星数据配准空间位置，在无人机影像上选择典型样区，利用非监督分类解译典型样区内各地物占比；随机选取一部分样区，利用样区内各地物占比和对应卫星遥感波段反射率，结合最小二乘法求解各地物端元的反射率；利用光谱分解模型和各地物端元的反射率，求解有效影像区域内各像元的植被覆盖度；利用剩余一部分样区数据校正植被覆盖度计算结果。本发明核心在基于无人机的端元反射率获取方法和植被覆盖校正模型，可有效提高大空间尺度植被覆盖的计算精度。

Description

一种结合无人机影像的大空间尺度植被覆盖度计算方法

技术领域

本发明属于遥感技术应用领域，特别涉及一种结合无人机影像的大空间尺度植被覆盖度计算方法。

背景技术

植被覆盖度为植被在地面的垂直投影面积与区域总面积的比值，是刻画地表植被覆盖状况的重要参数，是模拟大空间尺度范围内大气-植被-土壤系统水分循环、能量交换和生物化学循环等地表过程的重要基础性数据，在农业、林业、资源环境管理、土地利用、水文生态过程、水土流失调查、灾害风险监测、干旱监测等领域都有广泛的应用。遥感技术因其大范围的观测能力，已成为获取区域甚至全球植被覆盖度的重要手段。经过近几十年的发展，利用遥感技术的植被覆盖度估算方法已经取得了长足的进步，并形成了一批成熟的区域性、全球性产品，但在实际应用过程中，还存在以下几个问题：

(1)无人机等低空遥感估算精度高，但空间范围有限

随着近年来无人机等低空遥感技术的普及，其高灵活性、云下飞行、影像分辨率高和成本低等优势，使快速获取景观小尺度上的真实植被覆盖度有了可能，例如发明专利“一种基于无人机的景观尺度植被覆盖度的计算方法及系统”(申请号：201610913357.8)。但是，因续航能力、数据量等方面的问题，无人机只能获取景观小尺度范围内的植被覆盖度。

(2)以混合像元分解方法为代表的卫星遥感植被覆盖度估算方法具备一定的理论意义和应用基础，但端元的选择存在诸多问题

目前，遥感测算植被覆盖度的方法主要有回归模型法、训练模型法和混合像元分解法三类。其中，回归模型方法主要是通过多元统计建立植被覆盖度与波段反射率或者植被指数的经验关系，该方法虽然简单易行，但是仅适用于特定的区域和特定的植被类型，并且需要大量的高质量地面观测以保证精度。训练模型法是通过一些训练样本对学习模型进行训练，进而建立起计算植被覆盖度的模型，该方法的关键点和难点都在于训练样本的选择，在复杂的地表环境下，样本的准确性和代表性难以保证。混合像元分解法依据传感器收集的地面反射光谱信息是植被光谱与其它光谱的综合信息的原理，获取地表植被的贡献率，因其具备一定的物理意义、应用方便且精度较高，成为目前应用最为广泛的方法，如二分法。

混合像元分解方法包括端元(代表某一类地物的纯像元)光谱获取和丰度(混合像元内各类地物贡献率)计算两个步骤。其中，端元光谱的获取是最关键步骤，也是面临问题较多的步骤。常用的提取方法有两种，一种是通过地面观测或者是已有的波谱库，通常称为“参考端元”，另一种是在遥感影像中直接选择，称为“影像端元”。参考端元虽然理论上比较精确，但由于影像上地物光谱受到大气、地形和传感器等的影响，与实测或者地物波谱库中的地物波谱曲线存在很大差别，即使经过辐射纠正等预处理，也难以很好的匹配，同时基于地面实测的光谱需要投入较多的人力及物力。影像端元虽然有经济、简单和方便的优势，且可以保证端元与影像具有相同的度量尺度，但是在获取大空间尺度所需的中分辨率影像(如TM/ETM+,环境星HJ)上，因地物的复杂性可能并不存在“纯净”的端元，因此，不论采用什么样的方法提取影像端元，获取的光谱信息均不纯粹，导致后续丰度计算出现较大误差。

(3)受到植被类型的复杂性、不同生长状况和影像成像质量的影响，遥感估算的植被覆盖度存在偏差

遥感估算植被覆盖度的核心理论是不同地物类型的光谱反射率差异，不同植被类型，以及植被的不同生长期条件下，均可以导致光谱反射率的变化，造成植被的光谱反射率与其它地物类型难以区分，甚至出现“同物异谱”、“同谱异物”现象，造成植被覆盖度估算值存在偏差。同时，在遥感影像成像过程中，由于受大气散射和吸收的影响，改变了传感器接收到的地表反射能量，极大地影响着遥感光谱信息的提取精度，虽然现阶段已经开发出大量的大气校正模型，但是这些模型通常严重依赖地面观测和探空数据，在这些数据无法保障的情况下，校正的精度往往不高，尤其是在大空间尺度范围内，造成遥感估算的植被覆盖度无法满足应用级的需求。

鉴于遥感估算植被覆盖度在实际应用中存在的上述问题，研究一种可充分利用低空遥感和卫星遥感的优势，提高大空间尺度植被覆盖度估算精度的方法，不仅具有重要的理论意义，也同时具有重要的实际应用价值。

发明内容

本发明针对混合像元分解方法中存在的端元选择问题，以及植被类型复杂性和影像质量等引起的误差问题，在充分利用无人机技术优势的基础上，提供了一种大空间尺度植被覆盖度计算方法。该方法基于无人机影像和最小二乘法计算端元光谱反射率，利用光谱分解方法，计算大空间尺度范围内各像元的植被覆盖度，并结合无人机对可能存在的偏差提供校正，具有适应性强、精度高、易于操作等优点，可大量减少野外实地观测，因而具备较强的推广性。

本发明的目的通过以下的技术方案实现：一种结合无人机影像的大空间尺度植被覆盖度计算方法，包括步骤：

(1)对遥感影像进行预处理，计算NDVI，根据预定阈值获取有效影像区域；

(2)对无人机获得的图像进行拼接并获取目标区域正射投影影像，将目标区域正射投影影像、预处理后的遥感影像进行空间位置配准；在无人机影像上选择典型样区，解译各典型样区内各地物类型的丰度；

(3)随机选取一部分样区，结合样区内各地物类型的丰度和对应卫星遥感波段反射率，利用最小二乘法求解各地物端元的反射率；

(4)根据各地物端元的反射率，利用光谱分解方法，逐像元计算遥感影像中有效影像区域的植被覆盖度；

(5)评估及校正步骤：用步骤(3)随机选取后剩余的样区检验及评估步骤(4)获取的同区域的植被覆盖度；若评估不合格，则建立二者的线性拟合关系，利用该线性拟合关系校正步骤(4)获取的有效影像区域的植被覆盖度，进而得到最终应用级的植被覆盖度。

优选的，所述步骤(1)中，对遥感影像进行预处理，包括步骤：

(1-1)进行大气纠正，将影像的DN值转化为真实地表反射率；

(1-2)对大气纠正后的影像进行几何校正，获取空间精确定位的影像。

优选的，所述步骤(1)中，计算NDVI的公式如下：

NDVI(＝(ρ_NIR-ρ_R)/(ρ_NIR+ρ_R))；

其中，ρ_NIR、ρ_R分别为纠正后的近红外波段、红光波段的反射率。

优选的，所述步骤(2)中，对无人机获得的图像经过如下处理：

(2-1)对无人机获得的图像进行拼接，以获取目标区域正射投影影像，并获取正射影像的空间分辨率；

(2-2)基于同名点追踪的配准模型，对目标区域正射投影影像和预处理后的遥感影像进行空间位置配准；

(2-3)在无人机影像上筛选出若干典型样区，要求样区的长和宽均为遥感影像空间分辨率的整数倍，样区的数量为N，N≥2n，n为地物类型的数量；

(2-4)对筛选出的样区进行监督分类，获取各地物类型的丰度(也称为所占的比率)。

优选的，所述步骤(3)中，随机选取若干个样区，依据地表反射率是像元内的各端元地表反射率的线性组合的原理，得到公式：

R_j＝∑F_ir_ij+e_j

其中，R_j是某一样区遥感影像中第j波段的反射率，1≤j≤n；F_i为i类地物类型在某地块上的丰度，1≤i≤n；r_ij是第i类地物类型端元在j波段的反射率；e为随机误差。

利用最小二乘法，当满足|E|²＝e₁ ²+e₂ ²+...+e_n ²取值最小时，得到一组地物类型端元的反射率：

矩阵A即为需要求解的各地物类型端元在不同波段上的反射率。

更进一步的，在选取遥感影像的波段时，优先选取植被在波段内存在明显变化特征的波段，例如，植被存在一个低反射峰、高反射峰和吸收谷的绿光、近红外和红光波段。

优选的，由于地表的复杂性，以及不同遥感影像成像时间和条件差异，通过步骤(1)-步骤(4)计算的植被覆盖度，可能存在一些系统系的偏差，步骤(5)所述的线性拟合，主要目的是找到两者的线性拟合曲线，以便根据拟合的曲线校正偏差。

本发明通过一部分无人机影像，利用最小二乘法求解了与遥感影像相匹配的各端元的反射率，基于混合像元光谱分解的思路，计算了大空间尺度范围内的植被覆盖度，并利用剩余一部分无人机影像对光谱分解的植被覆盖度进行校正。与现有的其它研究技术相比，本发明具有如下优点和有益效果：

(1)本发明创新性的提出了一种高分辨率的无人机影像和中低分辨率遥感影像联动的大空间范围植被覆盖度的方法，具有理论基础扎实、可操作性强和普适性高的特点，既提高了无人机影像的应用范围，又提高了中低分辨率遥感影像。

(2)本发明提出的端元反射率计算方法，一方面计算过程简单、方便，不需要野外地面实测地物光谱，降低了野外调查所需的人力、物力成本，并且获取的端元反射率与遥感影像数据具有相匹配的度量尺度，从而可以提高估算精度；另一方面，不需要直接从图像中选择端元，从而提升了可应用的范围。

(3)本发明提出的基于无人机的混合像元分解模型校正方法，可有效的降低因植被类型的复杂性、不同生长状况和遥感影像成像质量等因素而造成的计算误差。

附图说明

图1是本发明方法的流程图；

图2是基于Landsat8遥感影像求解的植被覆盖度与无人机实测植被覆盖度的散点图；

图3是经过校正后的Landsat8遥感影像植被覆盖度与无人机实测植被覆盖度的散点图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。本领域技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本实例以Landsat8OLI为卫星遥感数据源，以我国南方某岩溶区为试验区，结合岩溶区地物特征，将地面植被类型分为土壤、植被和岩石三类。

步骤一，获取Landsat8影像，并对影像进行大气纠正，以消除大气散射引起的辐射误差影响，将影像的DN值转化为真实地表反射率。

对大气纠正后的影像进行几何校正，获取空间精确定位的影像。

步骤二，利用纠正后的近红外和红光波段的反射率(分别为ρ_NIR和ρ_R)，结合公式(1)计算NDVI，依据预先设定的阈值确定有效影像区域，剔除诸如水体等不需要计算植被覆盖度的区域；

NDVI＝(ρ_NIR-ρ_R)/(ρ_NIR+ρ_R) (1)

步骤三，使用专业软件拼接无人机拍摄的图片以获取目标区域正射投影影像，在ArcGIS等制图软件中空间配准无人机正射影像和Landsat8影像，在有效影像区域随机选取样区，样区的大小为30m*30m，且位置与Landsat8影像像元对应，样区数量为36(满足条件N≥2n)，利用监督分类获取样区内土壤、植被和岩石所占的比率，获取样区中Landsat8中绿、红和近红三个波段的反射率。

步骤四，随机选取N/2个(18个)样区，依据所选样区内三类地物所占比率和对应的Landsat8影像中绿、红和近红三个波段的反射率，利用最小二乘法(见公式(2)))，求解各个地物类型的端元反射率，

R_j＝∑F_ir_ij+e_j (2)

式中，R_j是Landsat8中绿、红和近红三个波段中第j波段的反射，1≤j≤3；F_i为i类地物类型端元在某地块上所占的比率(1≤i≤3)，由步骤三所述的监督分类方法获取；r_ij是土壤、植被和岩石中第i类地物类型端元在j波段的反射率；ej为第j波段的随机误差。e为随机误差。

利用最小二乘法，当满足|E|²＝e₁ ²+e₂ ²+e₃ ²取值最小时，可求得一组地物类型的端元反射率A：

步骤五，利用最小二乘法(见公式(2)，并设定e_j为0)，结合端元反射率A和Landsat8影像中三个波段的反射率，逐像元求解Landsat8影像中植被所占比率，求解过程类似于步骤四。具体地，在公式(2)中，R_j是已知的卫星遥感影像，e_j＝0，r_ij是端元的反射率，在这种状况下，如果地物类型不多于卫星影像的波段数，那么F_i是可以通过解多项式求解出来。上述求解多项式的过程就是光谱分解模型。

步骤六，利用步骤四剩余的N/2个样区的无人机解译结果，检验步骤五计算结果的精度，如果精度不够，则建立二者的线性拟合关系，根据两者线性拟合关系，校正步骤五获取的有效影像区域的植被覆盖度，进而得到最终应用级的植被覆盖度。

参见图2，为基于Landsat8遥感影像求解的植被覆盖度(估算值)与无人机实测植被覆盖度(真实值)的散点图，两者确定性系数(R²)为0.78，平均误差为0.076(估算值-真实值)，均方根误差为0.12。图中虚线表示参考曲线y＝x，理论上来说，散点均匀分布在参考曲线两侧说明估算效果较好，但是从图中看，散点多在该线的下方，说明估算值系统性偏大，需要进行校正。

为此，利用步骤五获得的Landsat8影像的植被覆盖度(估算值)，与无人机实测植被覆盖度(真实值)建立拟合关系，其关系为：

Y＝0.9303X-0.0331 (3)

式中，X代表估算值，Y代表真实值。在进行校正时，获得估算值之后，将估算值代入上述公式(3)，即可得到Landsat8遥感影像中有效影像区域内校正后的植被覆盖度Y’，即最终应用级的植被覆盖度。

参见图3，为经过校正后的Landsat8遥感影像植被覆盖度与无人机实测植被覆盖度的散点图，这里是将无人机实测植被覆盖度作为真实值Y’，以此数据为基准对校正后的数据进行验证。真实值Y和校正后的估算值Y’之间的关系是：

Y’＝Y-4e^-5；

二者线性拟合关系的斜率为1，截距接近于0，两者的平均误差无限接近于0，均方根误差为0.09。因此，校正后得到的估算值已经可近似认为是真实值，精度很高。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (5)

1.一种结合无人机影像的大空间尺度植被覆盖度计算方法，其特征在于，包括步骤：

(1)对遥感影像进行预处理，计算NDVI，根据预定阈值获取有效影像区域；

(2)对无人机获得的图像进行拼接并获取目标区域正射投影影像，将目标区域正射投影影像、预处理后的遥感影像进行空间位置配准；在无人机影像上选择典型样区，解译各典型样区内各地物类型的丰度；

(3)随机选取若干个样区，依据地表反射率是像元内的各端元地表反射率和丰度的线性组合的原理，得到公式：

R_j＝∑F_ir_ij+e_j

其中，R_j是某一样区遥感影像中第j波段的反射率，1≤j≤n，n为地物类型的数量；F_i为i类地物类型在某地块上的丰度，1≤i≤n；r_ij是第i类地物类型端元在j波段的反射率；e为随机误差；

利用最小二乘法，当满足|E|²＝e₁ ²+e₂ ²+...+e_n ²取值最小时，得到一组地物类型端元的反射率：

矩阵A即为需要求解的各地物类型端元在不同波段上的反射率；

(4)根据各地物端元的反射率，利用光谱分解方法，逐像元计算遥感影像中有效影像区域的植被覆盖度；

(5)评估及校正步骤：用步骤(3)随机选取后剩余的样区检验和评估步骤(4)获取的同区域的植被覆盖度；若评估不合格，则建立二者的线性拟合关系，利用该线性拟合关系校正步骤(4)获取的有效影像区域的植被覆盖度，进而得到最终应用级的植被覆盖度。

2.根据权利要求1所述的结合无人机影像的大空间尺度植被覆盖度计算方法，其特征在于，所述步骤(1)中，对遥感影像进行预处理，包括步骤：

(1-1)进行大气纠正，将影像的DN值转化为真实地表反射率；

(1-2)对大气纠正后的影像进行几何校正，获取空间精确定位的影像。

3.根据权利要求1所述的结合无人机影像的大空间尺度植被覆盖度计算方法，其特征在于，所述步骤(1)中，计算NDVI的公式如下：

NDVI(＝(ρ_NIR-ρ_R)/(ρ_NIR+ρ_R))；

其中，ρ_NIR、ρ_R分别为纠正后的近红外波段、红光波段的反射率。

4.根据权利要求1所述的结合无人机影像的大空间尺度植被覆盖度计算方法，其特征在于，所述步骤(2)中，对无人机获得的图像经过如下处理：

(2-1)对无人机获得的图像进行拼接，以获取目标区域正射投影影像，并获取正射影像的空间分辨率；

(2-2)基于同名点追踪的配准模型，对目标区域正射投影影像和预处理后的遥感影像进行空间位置配准；

(2-3)在无人机影像上筛选出若干典型样区，要求样区的长和宽均为遥感影像空间分辨率的整数倍，样区的数量为N，N≥2n，n为地物类型的数量；

(2-4)对筛选出的样区进行监督分类，获取各地物类型的丰度。

5.根据权利要求1所述的结合无人机影像的大空间尺度植被覆盖度计算方法，其特征在于，在选取遥感影像的波段时，选取植被存在一个低反射峰、高反射峰和吸收谷的绿光、近红外和红光波段。