一种中分辨率遥感影像中塑料大棚提取及盖度反演的方法与流程

本发明涉及一种塑料大棚提取的新方法，特别是针对利用中分辨率遥感影像进行塑料大棚空间覆盖的提取以及影像像元内塑料大棚盖度的反演，属于图像处理领域。

背景技术：

自19世纪50年代，日本及欧洲国家利用塑料膜进行植被温床栽培的实验取得成功以来，塑料大棚(plasticgreenhouse，pg)农业种植获得了空前发展，并迅速遍及世界各地。截止到2016年，全世界pg农业的种植面积已经达到了3.019×10⁶公顷。参见briassoulis,d.,dougka,g.,dimakogianni,d.,vayas,i.,2016.analysisofthecollapseofagreenhousewithvaultedroof.biosystemsengineering.151,495-509,(《拱形大棚倒塌分析》，briassoulis,d等，生物系统工程杂志，2016年)。作为传统农业向精细化农业的一种转变方式，塑料大棚材料具有良好的保温保湿作用，从而可以有效协助农作物避免不利的生长条件，提高农作物的生产量。然而，大量的塑料制品的使用给区域环境带来了极大的压力，例如产生大量的塑料垃圾、破坏农业种植区的生物多样性、造成不良的视觉影响等。

此外，由于塑料材质特殊的物理特征(光学透射性、气密性以及高反射性)，pg农业的大范围覆盖能够改变能量与水汽循环，从而影响区域气候。其对气候的影响机理主要表现在三个方面。1)增强地表的反射率：pg农业中广泛使用的白色透明塑料膜具有较高的反射率，因此可以大大增加地表反射，减少地表吸收的能量。2)减少pg内部辐射能量的向外扩散：来自pg内部农作物或土壤反射的太阳光能量或地物发射的长波辐射能量会部分被反射回pg底部，经多次散射或反射而最终被吸收，从而减少pg内能量向大气中的扩散。3)隔断蒸散发向大气中的扩散：塑料膜具有很好的气密性，能够完全阻止植被蒸腾或地表蒸发作用产生的水分向大气中的扩散，从而达到保湿的效果。基于以上几方面的作用，pg农业覆盖区上空极易形成干燥的上升气流，从而影响当地的气候(降水、温度等)，参见campra,p.,garcia,m.,canton,y.,palacios-orueta,a.,2008.surfacetemperaturesoutheasterncoolingtrendsandnegativeradiativeforcingduetolandusechangetowardgreenhousefarminginspain.journalofgeophysicalresearch:atmospheres.113(d18)(《西班牙大棚农业发展引起地表温度东南冷却及负面辐射强迫》，campra,p等，大气圈杂志，2008年)。因此，准确的掌握pg农业的空间分布以及时间变化对于最大化其经济效益同时最小化其环境影响具有深远意义。

近年来，利用高分辨率影像进行pg农业的监测与制图受到了一定程度的重视。譬如，levin等(2007)利用实测光谱数据对pg农业材料的光谱特征进行了分析；并前瞻性的研究了利用遥感技术对pg农业区进行检测与制图的可行性。参见levin,n.,lugassi,r.,ramon,u.,braun,o.,ben-dor,e.,2007.remotesensingasatoolformonitoringplasticultureinagriculturallandscapes.internationaljournalofremotesensing.28(1),183-202(《以遥感为手段对塑料农业进行检测》，levin等，国际遥感杂志，2007年)。koc-san等(2013)基于高分辨率影像，评价了不同分类方法对于pg农业提取的效果。参见koc-san,d.,2013.evaluationofdifferentclassificationtechniquesforthedetectionofglassandplasticgreenhousesfromworldview-2satelliteimagery.journalofappliedremotesensing.7(1),073553-073553(《评价不同分类技术对于从worldview-2遥感影像中检测玻璃与塑料大棚的效果》，koc-san，应用遥感杂志，2013年)。aguilar等(2014)提出了一种面向对象的方法用于高分辨率影像的pg检测与制图。参见aguilar,m.a.,bianconi,f.,aguilar,f.j.,fernández,i.,2014.object-basedgreenhouseclassificationfromgeoeye-1andworldview-2stereoimagery.remotesensing.6(5),3554-3582(《一种从geoeye-1和worldview-2立体影像对中提取大棚的面向对象分类方法》，aguilar等，遥感杂志，2014年)。但是目前高分辨率遥感影像的获取较为困难，数据公开程度低，空间覆盖范围有限，导致难以对pg区域进行实时的大范围检测。因此，具有更大空间覆盖范围与更高时间观测频率的中分辨影像被认为更有利于及时准确的大范围pg监测。例如，zhang等(2004)基于landsat影像发展了一个pg农业提取指数(vi)，并利用该指数对山东青州的塑料农业区进行了提取。参见zhao,g.,li,j.,li,t.,warner,t.,2004utilizinglandsattmimagerytomapgreenhousesinqingzhou,shandongprovince,china.pedosphere,14(3),363-369。(《在中国山东青州区域，利用landsattm影像提取大棚》，zhao等，土壤圈杂志，2004年)。此外aguilar等(2015)通过结合面向对象方法与决策树影像分类方法，实现了对pg农业的大区域检测与制图。参见aguilar,m.a.,vallario,a.,aguilar,f.j.,lorca,a.g.,parente,c.,2015.object-basedgreenhousehorticulturalcropidentificationfrommulti-temporalsatelliteimagery:acasestudyinalmeria,spain.remotesensing.7(6),7378-7401。(《利用多时期遥感影像进行面向对象的大棚农作物提取》，aguilar等，遥感杂志，2015年)。然而，现存的方法缺乏对pg光谱变异性的充分考虑；且都只能得到二值化的pg提取结果(即影像像元内是否存在pg)。考虑到中分辨率遥感影像空间分辨率的不足，估测影像像元内pg盖度显得异常重要。

技术实现要素：

为此本发明提供了一种中分辨率遥感影像中塑料大棚提取及盖度反演的方法，该方法可以有效减少或避免前面提到的问题。

为解决上述问题，本发明提供的一种中分辨率遥感影像中塑料大棚提取方法，所述方法包括如下步骤：

步骤a，对遥感影像进行辐射定标与大气纠正，获取地表的真实反射率。

步骤b，利用ndvi(归一化植被指数，normalizeddifferencevegetationindex)与ndbi(归一化城建指数，normalizeddifferencebuildingindex)阈值进行影像的初步处理，去除影像中的浓密植被(例如：森林)及大部分人造地表区域(城市、道路等)。其中所使用的阈值为ndvi>0.73，ndbi>0.005，两者其一符合上述阈值的像元即被认为是非pg覆盖区域；ndvi与ndbi的定义如下：

其中，rred、rnir、rsw1分别为红光(red)、近红外(nearinfrared,nir)与短波红外(short-waveinfrared-1,sw1)波段的真实地表反射率。

步骤c，对每一个像元定义一个pgi指数，利用所述像元的pgi指数对步骤b处理后的遥感影像进行塑料大棚区域的提取，所述pgi指数定义如下：

其中，rblue、rgreen、rred、rnir分别为蓝光(blue)、绿光(green)、红光、近红外段的真实地表反射率。当像元的pgi∈[1.3,6.7]时，认为该像元内存在pg；否则认为该像元内不存在pg。

优选地，所述像元的上述阈值所对应的可检测塑料大棚的盖度范围为12％～100％。

另外，本发明还提供了一种中分辨率遥感影像中塑料大棚盖度反演的方法，所述方法包括如下步骤：

步骤a，对遥感影像进行辐射定标与大气纠正，获取地表的真实反射率；

步骤b，利用ndvi与ndbi阈值进行影像的初步处理，去除影像中的浓密植被及大部分人造地表区域；其中所使用的阈值为ndvi>0.73，ndbi>0.005，两者其一符合上述阈值的像元即被认为是非塑料大棚覆盖区域；ndvi与ndbi的定义如下：

其中，rred、rnir、rsw1分别为红光、近红外与短波红外波段的真实地表反射率；

步骤c，对每一个像元定义一个pgi指数，利用所述像元的pgi指数对步骤b处理后的遥感影像进行塑料大棚区域的提取，所述pgi指数定义如下：

其中，rblue、rgreen、rred、rnir分别为蓝光、绿光、红光、近红外段的真实地表反射率；

当像元的pgi∈[1.3,6.7]时，认为该像元内存在塑料大棚；否则认为该像元内不存在塑料大棚；

步骤d，在步骤c判断像元内存在塑料大棚的图像区域，利用下述盖度反演公式计算像元内塑料大棚的盖度fpg，所述盖度反演公式定义如下：

fpg＝0.2405×ln(rpgi)+0.7687

本发明具有的有益效果是:

基于本发明的pgi指数与pg盖度反演模型，本发明可以有效利用中分辨率遥感影像进行pg覆盖区的提取以及像元内pg盖度的反演，实现更准确的pg农业监测与制图。

本发明的其他特征和优点将在随后的说明书中阐述，并且部分的从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。

附图说明

以下附图仅旨在与对本发明做示意性说明和解释，并不限定本发明的范围。其中，

图1为根据本发明的一种具体实施的利用中分辨率遥感影像进行pg提取及盖度反演方法的流程示意图；

图2a～2e为本发明用于指数设计及验证的样本研究区域数据；

图3a～3b为确定ndvi与ndbi阈值的模拟实验结果；

图4为不同地物的光谱特征曲线；

图5a～5h为pg覆盖区像元各波段的反射率对像元内pg盖度的敏感性；

图6a～6i为不同地物类型间的反射率可分离性；

图7a～7b为利用混合像元模拟实验得到的pgi有效阈值；

图8a～8d为本发明的盖度反演模型，以及利用统计回归对模型表现的评价结果；

图9a～9c为本发明的pgi指数与pg盖度反演模型在样本区的表现；

图10a～10e为本发明的pgi指数与pg盖度反演模型在大试验区域的表现。

具体实施方式:

为了对本发明的技术特征、目的和效果有更加清楚的理解，现说明本发明的具体实施方式。但本领域的技术人员应该知道，以下实施例并不是对本发明技术方案的唯一限定，凡是在本发明技术方案精神实质下所做的任何同等变换或改动，均应视为属于本发明的保护范围。

图1为根据本发明的一个具体实施例的利用中分辨率遥感影像进行pg提取及盖度反演方法的流程示意图；参照图1所示，下面详细说明根据本发明的中分辨率遥感影像中塑料大棚提取及盖度反演的方法，即，所述方法涉及两个方面的内容，即方法一：塑料大棚的提取方法，以及方法二：塑料大棚的盖度反演的方法。其中，方法二的塑料大棚的盖度反演的方法是在方法一的塑料大棚的提取方法实施之后，在提取出来塑料大棚的像素内进行的盖度反演，图1显示的是方法一和方法二整体流程示意图。

即方法一的塑料大棚的提取方法包括三个步骤，分别为：

步骤a，影像辐射定标与大气纠正，获取地表真实反射率；

步骤b，ndvi与ndbi阈值初步处理，去除浓密植被及大部分人造地表；

步骤c，pgi指数处理，提取pg区域。

相对于方法一，方法二的塑料大棚的盖度反演的方法在方法一的上述三个步骤之后，增加了盖度反演的步骤，即，方法二的塑料大棚的盖度反演的方法包括四个步骤，分别为：

步骤a，影像辐射定标与大气纠正，获取地表真实反射率；

步骤b，ndvi与ndbi阈值初步处理，去除浓密植被及大部分人造地表；

步骤c，pgi指数处理，提取pg区域。

步骤d，pg盖度反演。

本发明旨在利用pg内部农作物生长茂盛时期的遥感影像进行pg提取及盖度反演。因为此时期的pg具有与其它地表覆盖类型更可分的光谱特征。在本研究案例中，我们选取2012年4月17日的一幅美国陆地卫星(landsat，参见图10a)影像进行指数发展的阐述及其效果的检验。

下面分别对上述步骤a，b，c，d分别予以详细介绍。

步骤a，影像辐射定标与大气纠正，获取地表真实反射率。

对于遥感影像而言，在进行pg提取之前，首先要对遥感影像进行辐射定标与大气纠正。本发明采用envi软件的辐射定标与大气纠正模块(flaash)对遥感影像进行辐射定标及大气纠正处理，去除大气影响，并获得地表真实反射率。

步骤b，ndvi与ndbi阈值初步处理，去除浓密植被及大部分人造地表。

本发明利用预设的ndvi与ndbi阈值(ndvi>0.73，ndbi>0.005)分别去除浓密植被及大部分的人造地物。当像元内的ndvi或ndbi符合上述阈值范围时，即认为像元内不存在pg，不再参与后续步骤的处理(说明:只要像元内的ndvi与ndbi二者其一符合阈值范围，即认为像元内不存在pg)。对于ndvi与ndbi阈值的确定，本发明采用三端元的混合像元模拟计算多种地物类型组合下的像元ndvi与ndbi，然后根据pgi指数的最小可检测盖度(12％)确定ndvi与ndbi的阈值。详细步骤阐述如下：

1.选取典型pg覆盖区域(图2a)，获取该区域的中分辨率遥感影像(landsat，30米空间分辨率)与高分辨率遥感影像(quickbird-2全色影像，0.6米空间分辨率，图2b)。基于以上两幅影像，提取pg、植被、裸土、人造地表四种地物类型的纯像元(端元)光谱。

2.利用线性光谱混合模型，分别根据公式1、2模拟以下两种端元组合类型下的混合光谱:

端元组合(1)pg+植被+裸土：

r1(λ)＝fv×rv(λ)+fs×rs(λ)fv+fs＝1,fv＝0,0.1,......,0.9,1.0

rmix(λ)＝f1×r1(λ)+f2×rpg(λ)f1+f2＝1,f1＝0,0.1,......,0.9,1.0.(1)

其中r1为植被(rv)与裸土(rs)组合的光谱；fv与fs为植被与裸土各自在r1中所占的比例；rmix为设计场景下的混合光谱；f1与f2为混合场景下r1与pg各自所占的比例；λ为波长。

端元组合(2)pg+植被+人造地表：

r1(λ)＝fv×rv(λ)+fm×rm(λ)fv+fm＝1,fv＝0,0.1,......,0.9,1.0

rmix(λ)＝f1×r1(λ)+f2×rpg(λ)f1+f2＝1,f1＝0,0.1,......,0.9,1.0(2)

其中，rm为人造地表光谱；fm为人造地表在r1中所占的比例。

3.计算端元组合(1)所有模拟场景下的ndvi以及端元组合(2)所有模拟场景下的ndbi，结果如图3所示。

4.根据图3，当像元内pg盖度大于pgi指数最小可检测盖度(12％)时，ndvi的最大值为0.73；ndbi的最大值为0.005。基于此可知，当像元pg盖度大于12％时，其ndvi与ndbi不可能大于以上两个阈值。因此可以利用以上两个阈值分别去除浓密植被以及人造地表覆盖。

步骤c，pgi指数处理，提取pg区域

本发明根据pg独特的光谱特征发展了pgi指数，用于从步骤b处理后的遥感影像中识别出pg区域；并使得具有不同pg盖度的像元显示不同的pgi值(像元内pg盖度越大，pgi值越大)。在本发明中，当像元pgi∈[1.3,6.7]时，认为该像元内存在pg；否则认为该像元内不存在pg。该pgi阈值对应的可检测pg盖度范围为12％—100％。

对于pgi指数的设计以及pgi阈值的确定过程，本发明阐述如下：

1.选取典型pg覆盖区域(图2a)，获取该区域的中分辨率遥感影像(美国陆地卫星，landsat，30米空间分辨率)与高分辨率遥感影像(快鸟卫星，quickbird-2全色影像，0.6米空间分辨率，图3b)。在quickbird-2影像上人工提取pg区域(图2c)，并制作30米空间分辨率的pg盖度参考图(图2d)。

具体来说，对于参考影像中盖度大于0的像元，认为像元中存在pg；否则认为像元中不存在pg。对于影像中的非pg像元，利用监督分类器(譬如支持向量机)进行分类，分类类别主要包含：露天农田、裸体、人造地表及水体。(图2e)。该分类结果为后边的可分离性分析提供基础。可分离性分析中统计的各类别即为该处得到的分类结果。监督分类器是利用人工选取的样本对影像进行分类的一类方法，对应于非监督分类器。支持向量机是监督分类器的一种方法，是公知的分类方法。

2.基于以上数据，提取不同地物类型的光谱曲线(图4，说明：对于pg区域分，别提取像元内pg盖度为30％，50％，70％，90％的光谱曲线)，并进行光谱特征分析。得出pg的光谱特征：①塑料地膜的覆盖可以整体提高地表的反射率；②pg像元的反射率随pg盖度的增大而增大。

3.基于典型pg区域提取的pg盖度图(图2d)与landsat影像(图2a)，建立像元内pg盖度与光谱反射率之间的统计关系，分析光谱反射率对pg盖度的敏感性。结果如图5所示，由此得出pg的光谱特征：③相较于其它波段，蓝光、绿光、近红外对pg盖度具有较高的敏感性(三者的线性回归斜率分别为0.1309，0.1302，0.1214)。

4.基于典型区域的分类图(图2e)与landsat影像，统计不同地物类型的波段反射率分布，进行不同地物类型间的光谱可分离性分析(图6)。得出pg的光谱特征：④近红外与红光波段的差值(nir-red)可以有效的将pg与植被从裸土中分离出来(图6g)；⑤蓝光波段对区分pg与植被具有最好的效果(图6a和b)。

说明：鉴于蓝光与绿光对于区分pg与植被效果相似，我们采用m指数准确判别两者之中表现更好的一个。m指数越大，可分离性越高。其定义如下：

其中μc1与μc2分别为c1与c2两个类别统计分布的均值，σc1与σc2分别为两个类别统计分布的方差。参见kaufman,y,j.,remerl,a.,1994.detectionofforestsusingmid-irreflectance:anapplicationforaerosolstudies.ieeetransactionsongeoscienceandremotesensing.32(3),672-683。(《利用中红外波段检测森林：一种气溶胶研究应用》，kaufman和remer,地球科学与遥感学报，1994年)。

5.基于以上分析得到的光谱特征，定义pgi指数如下：

其定义的原理如下：

在上述公式的分子中，pg像元的rblue要大于露天农田、裸土以及人造地物(如图4所示)；对于rnir-rred，pg像元与露天农田像元具有相似的值，且要远大于裸土与人造地物。因此pg像元的rblue与rnir-rred乘积要大于其它地物类型。对于上述公式中的分母，由于pg覆盖可以整体提高地表反射率，pg像元的rblue+rgreen+rnir要大于其它地物类型。与1差值之后，pg像元即在分母部分表现出小值。综合分子、分母可知相对于其它地物类型，pg像元具有较高的pgi值；且像元内pg盖度越大，像元的pgi值越大。

6.利用公式1所述的方法模拟的不同混合场景下的光谱，计算所有场景下的pgi，并与pg盖度建立相关关系。根据图7可知，像元pgi随pg盖度的增大而增大。当像元pgi出现在0.6-1.3范围内时，该值的出现可能由于像元内存在pg，也可能由于露天农田与裸土的混合导致。因此，当像元pgi处在该范围内时，认为不确定像元内是否存在pg。但当像元pgi处于1.3-6.7范围内时，可认为像元内一定存在pg。因为露天农田与裸土的混合不可能导致pgi大于1.3。由此可以得出pgi的阈值范围为pgi∈[1.3,6.7]，其对应的可检测像元内pg的最小盖度为12％。

经过c步骤的处理即可获知哪些像元内一定存在pg；并去除不存在pg的像元以及不确定是否存在pg的像元之后，可以进行下一步的盖度反演。应当指出的是，通过步骤a～步骤c，可以从中分辨率遥感影像中，去除不存在pg的像元以及不确定是否存在pg的像元之后，剩下的像元即为确定存在pg的区域，从而可以完成方法一定塑料大棚的全部提取过程。

步骤d，pg盖度反演

受中分辨率遥感影像分辨率不足的影响，像元内经常会出现多种地物的混合。因此，为了实现更准确的pg覆盖检测，本发明发展了一个回归模型用以预测像元内的pg盖度。该模型的建立过程阐述如下：

1.将典型pg覆盖区域均分为两部分(如图2d所示)。将上半部的区域ⅰ用于建立盖度反演模型；下半部的区域ⅱ用于模型效果检验。计算典型区域的pgi值，并进行阈值处理，提取pg区域(图9b)。

2.pg内部植被的生长状况在时间与空间上都会存在一定的差异。为了减弱pg内部植被生长状况差异对pg盖度反演结果的影响，我们移除了pgi指数中的nir-red(近红外与红光波段的差值，rnir-rred)及常数项部分(pgi公式中的常数100)，建立了一个回退的pgi指数rpgi(retrogressivepgi)。此后，在区域ⅰ上分别建立pgi指数、rpgi指数与参考pg盖度的统计关系；并分别利用线性模型与对数模型进行拟合(图8a与8b)。以期在两个模型中选择表现较好的模型用于pg盖度的预测。得到的定量关系模型如下：

模型一：

fpg＝0.1546×pgi-0.0759

模型二：

fpg＝0.2405×ln(rpgi)+0.7687

3.利用区域ⅱ对以上两个模型的pg盖度预测效果进行检验。首先利用上述两个模型分别预测区域ⅱ的pg盖度(图9c)，其后建立预测pg盖度与参考pg盖度的回归关系，并与1：1线进行比较(图8c和d)。结果显示相对于基于pgi的模型一(r²＝0.7624，rmse＝0.088)，基于prgi模型二的pg盖度预测结果与真实的pg盖度更接近(r²＝0.7874，rmse＝0.083)。由此可知，rpgi模型对于pg盖度的预测更为准确。因此，本发明选择基于rpgi的模型二用于pg盖度的反演。

根据本发明提供的以上步骤，即可对中分辨率遥感影像的pg区域进行提取，并对像元内的pg盖度进行估算。

为了更好的说明本发明的技术效果，针对一幅landsatenhancedthematicmapperplus(etm+)影像，分别利用本发明提出的一种利用中分辨率遥感影像进行塑料大棚提取及盖度反演的方法以及zhao等(参见zhao,g.,li,j.,li,t.,warner,t.,2004utilizinglandsattmimagerytomapgreenhousesinqingzhou,shandongprovince,china.pedosphere,14(3),363-369。(《在中国山东青州区域，利用landsattm影像提取大棚》，zhao等，土壤圈杂志，2004年)提出的vi指数对该区域的pg进行了提取。其中vi指数定义如下：

结果如图10所示，可以看出，pgi指数在pg区域表现出较高的值；而露天农田、裸土、人造地表以及水体都表现出相对较低的值。图10c为利用本发明提出的盖度反演模型估算得到的像元内pg盖度，与真实的塑料农田覆盖区域具有较好的一致性。相反，vi指数(图10d)未能有效的将pg从其它地物类型当中凸显出来，利用其方法提出的pg区域出现了较大的误差(图10e)。这充分说明本发明提供的方法可以有效的利用中分辨率影像对pg覆盖区域进行提取及盖度估算，实现更准确的pg覆盖检测。

本领域技术人员应当理解，虽然本发明是按照多个实施例的方式进行描述的，但是并非每个实施例仅包含一个独立的技术方案。说明书中如此叙述仅仅是为了清楚起见，本领域技术人员应当将说明书作为一个整体加以理解，并将各实施例中所涉及的技术方案看作是可以相互组合成不同实施例的方式来理解本发明的保护范围。

以上所述仅为本发明示意性的具体实施方式，并非用以限定本发明的范围。任何本领域的技术人员，在不脱离本发明的构思和原则的前提下所作的等同变化、修改与结合，均应属于本发明保护的范围。