基于Logistic回归分析的塑料大棚遥感指数构建

陈俊， 沈润平， 李博伦， 遆超普， 颜晓元， 周旻悦， 王绍武

(1.南京信息工程大学地理科学学院，南京 210044； 2.中国科学院南京土壤研究所土壤与农业可持续发展国家重点实验室，南京 210008)

0 引言

根据第三次全国农业普查结果，2016年末全国塑料大棚占地面积为98.1×104hm2，比2006年增长了111.0%[1]。塑料大棚大幅提高了蔬菜作物产量，但也加剧了“白色污染”。因此，准确监测塑料大棚的分布范围，对区域农业的可持续发展和对环境影响的分析至关重要。

利用不同空间分辨率的遥感影像检测设施菜地、地膜覆盖等土地利用，已经引起越来越多的关注[2]。基于高空间分辨率遥感影像(0.5～2 m)进行塑料大棚提取是目前常用的检测方法，例如，Agüera等[3]利用QuickBird影像，确定塑料大棚提取的最佳波段组合为绿光、蓝光和近红外波段； Aguilar等[4]根据0.5 m空间分辨率GeoEye-1和WorldView-2影像，提出一种面向对象的塑料大棚分类方法。针对大面积的塑料大棚提取，国内外学者探讨了中等空间分辨率影像(2～30 m)用于塑料大棚提取的可行性[5]。例如，国内最早进行尝试的Zhao等[6]基于Landsat TM影像，利用指数方法在山东省进行塑料大棚的制图； Novelli等[7]比较了Sentinel-2 MSI和Landsat8 OLI影像的温室检测性能； Yang等[8]基于Landsat ETM+影像，利用指数法提取了山东省潍坊市塑料大棚面积。

目前，我国塑料大棚遥感提取的研究主要位于塑料大棚集中分布的北方地区。相比之下，太湖流域河网密布，景观破碎化程度高，土地覆盖类型复杂，塑料大棚分布相对离散。而现有塑料大棚指数计算时所需的波段数量较少，难以满足复杂地表覆盖下准确识别塑料大棚的要求。Logistic回归分析是一种根据单个或多个自变量，分析和预测因变量的多元分析方法，是目前常用的处理分类因变量的统计分类模型[9]。因此，本文基于Landsat8影像，结合4种遥感指数，即归一化植被指数(normalized difference vegetation index，NDVI)、归一化建筑指数(normalized difference building index，NDBI)、归一化裸土指数(normalised difference bareness index，NDBaI)以及改进的归一化水体指数(modified normalized difference water index，MNDWI)，经可分离性筛选后，运用Logistic回归模型，设计新塑料大棚遥感指数(new plastic greenhouse index, NewPGI)，以期更好地提取塑料大棚的分布信息。

1 研究区概况和数据源

1.1 研究区概况

本文选取的研究区为位于太湖流域西北部的常州市，如图1所示，地理坐标为N31°49′，E119°50′，行政区划包括新北区、钟楼区、天宁区、武进区、金坛区和溧阳市。研究区属于亚热带季风气候区，盛夏高温多雨，严冬寒冷干燥，以塑料大棚为主导的设施菜地种植模式，成为该地区农业生产最主要的发展方向。经过实地调研确定该地区蔬菜大棚的建材以透明塑料薄膜为主，根据Yang等[8]研究，选在作物生长季节提取塑料大棚信息。

(a) 研究区相对位置示意图 (b) Landsat8 B5(R),B4(G),B3(B)假彩色合成影像

图1 研究区位置及遥感影像

Fig.1Locationandimageofstudyarea

1.2 数据源及其预处理

采用覆盖整个研究区的2014年3月16日Landsat8影像(空间分辨率为30 m)提取常州市塑料大棚分布信息。利用ENVI 软件中FLAASH模块对OLI影像进行辐射校准和大气校正预处理，以消除大气影响； 对TIRS影像的热红外波段值进行归一化处理，使其与OLI影像反射率范围一致，以便后续处理。借助资源三号(ZY-3)影像，选取常州市塑料大棚相对集中区域作为样本区域，用于构建和验证塑料大棚指数(图2(a)，大小为9 km × 9 km)。采用2014年4月6日ZY-3影像作为提取样本区域塑料大棚的参考影像。基于ENVI5.3软件，利用ZY-3全色影像(图2(b))及30 m空间分辨率数字高程模型(digital elevation model,DEM)数据，对ZY-3多光谱影像(图2(c))进行正射校正，以纠正因系统因素或地形因子引起的几何畸变。对ZY-3多光谱及全色影像采用Gram-Schmidt光谱锐化法进行影像融合，融合后影像空间分辨率为2 m。

(a) Landsat8 B5(R),B4(G),(b) ZY-3全色影像(c) ZY-3 B3(R),B2(G),B1(B) B3(B)假彩色合成影像彩色合成影像

图2 样本区域遥感影像数据

Fig.2Satellitedatainthesamplearea

1.3 验证数据处理

1.3.1 样本区域分类参考图制作

利用eCognition影像分析软件，基于地理空间信息和专家知识，采用面向对象分类方法对融合后的ZY-3影像分类，并通过目视解译修正分类结果作为样本区域分类参考(图3(a))。将样本区域典型土地覆盖类型分为5类，包括塑料大棚、人造地表、裸地和休耕地、水体以及植被，得到样本区域土地覆盖分类结果(图3(b))。

(a) 塑料大棚分类参考(b) 土地覆盖分类

图3 样本区域分类参考

Fig.3Referenceofclassificationinthesamplearea

将塑料大棚类型所占像元比例超过50%的像元作为塑料大棚像元，通过IDL语言按照30 m×30 m像元尺寸,对图3(b)进行重采样，生成塑料大棚分类参考图，空间分辨率为30 m，用于验证基于Landsat8影像的塑料大棚提取结果。

1.3.2 研究区验证样本

基于Google Earth影像进行整个研究区塑料大棚分类的精度验证。随机选取2014年3月16日常州市塑料大棚、人造地表、裸地和休耕地、水体以及植被等土地覆盖类型(图4)，共抽取2 466个参考像元并归类为“塑料大棚像元”和“非塑料大棚像元”，生成感兴趣区域(region of interest,ROI)，构建塑料大棚分类的验证样本，用于验证研究区塑料大棚的分类精度。

(a) 研究区验证样本 (b) 验证样本示例1(c) 验证样本示例2

图4 研究区验证样本

Fig.4Validationdataofstudyarea

2 研究方法

2.1 塑料大棚光谱特征分析

识别塑料大棚的光谱特征对于设计合理的塑料大棚指数至关重要。本文基于Landsat8影像和样本区域土地覆盖分类，使用样本量计算工具(http: //fluidsurveys.com/university/survey-sample-size-calculator)，按照95%置信水平，随机抽取5种地物类型共385个像元作为研究样本，每种地物类型各选取77个像元。基于表1所示的Landsat8影像7个OLI多光谱波段(B1—B7)、1个TIR热红外波段(B10)以及常用于反映植被、人造地表、裸土和水体信息的4个遥感指数，即NDVI，NDBI，NDBaI和MNDWI，得出不同土地覆盖类型光谱的平均值(图5)。

表1 基于Landsat8影像的光谱数据及多种遥感指数Tab.1 Spectral data based on Landsat8 image and various remote sensing indexes

图5 不同土地覆盖类型光谱曲线(平均值)Fig.5 Spectral curves of differentland cover types (mean values)

塑料大棚作为一种人造设施，因在农作物之上覆盖一层白色塑料薄膜，削弱了作物的植被信息，使其同时兼具覆盖少量植被的土壤和人造地表双重的光谱特征。根据图5不难发现，塑料大棚与裸地和休耕地、人造地表的光谱特征十分相似，区分难度比较大，佐证了复杂土地覆盖类型下，塑料大棚遥感提取的影响因子众多，需要考虑加入更多的波段信息的设想。

2.2 塑料大棚光谱可分离性分析

塑料大棚指数的设计应该考虑最大限度地区分塑料大棚与其他土地覆盖类别。为实现这一目标，本文基于Kaufman等[13]研究，采用可分离性指标M，逐一比较塑料大棚与其他土地覆盖类型(人造地表、裸地和休耕地、植被以及水体)的分离度。该指标被定义为2种地物类型光谱曲线平均值μ之间的差异，按照标准差δ之和进行归一化处理，即

(1)

式中：μPG为塑料大棚反射率均值；μ1～μ4分别为人造地表、裸地和休耕地、植被以及水体的反射率均值；δPG为塑料大棚光谱标准差；δ1～δ4分别为人造地表、裸地和休耕地、植被以及水体的光谱标准差。M≥1，表示分离性好；M<1，表示分离性较差。

塑料大棚与其他4种土地覆盖类型之间的可分离性指标M如表2所示。根据分离标准(M≥1)，确定区分塑料大棚与人造地表的最佳波段是B3和B5； 区分塑料大棚与裸地和休耕地的最佳波段是B1和B2； 区分塑料大棚与植被的最佳波段和指数是B1，B2，B3，B4，B7，NDVI和MNDWI； 塑料大棚与水体可通过除NDBI外的其他波段和指数进行区分。综合考虑塑料大棚与其他地物类型区分度，加入更多波段及指数以涵盖不同土地覆盖类型，最终选定B1—B7，B10，NDVI，NDBaI和MNDWI共计11个参数，用于构建NewPGI。

表2 塑料大棚与典型土地覆盖类型可分离性指标MTab.2 Separation degree of plastic greenhousesfrom typical land cover types

2.3 Logistic回归分析

本文引入的Logistic回归模型属于一种机器学习分类算法。进行统计分析时，自变量可以是连续变量，也可以是离散变量，且不需要呈正态分布，增强了模型的应用范围和灵活性[14]，其最大优势是能够对众多影响因素进行拟合分析，通过机器学习确定众多自变量最佳回归系数。参考王济川等[15]和Cox[16]的研究，采用二分类因变量Logistic回归模型。

作为一种概率型的非线性回归模型，Logistic回归是一种研究二分类观察结果y与诸多影响因素Xk(k=1,2,…，n)之间关系的多变量分析方法。当输入测试样本数据时，Logistic回归分类器就是一组权值(a0,a1,…,an)，这组权值与测试数据按照线性加和得到预测值z，即

(2)

式中：Xk(k=1,2,…，n)为每个样本的n个特征； a0为常量。

Logistic回归模型的关键在于通过预测值z定义判定边界，以此确定样本的类型。通过单位阶跃函数，可表示预测值z与二分类观察结果y之间的关系，即

(3)

z为模型的判定边界，用于确定哪些样本是正样本，哪些为负样本。当z> 0时，函数值为1，判定为正样本； 当z< 0时，函数值为0，判定为负样本； 当z= 0时，函数值为0.5，表示样本为正样本或负样本概率相同，可任意判断。但是单位阶跃函数是分段式非连续性函数，无法应用于实际问题。因此，需要引入一个连续性函数——Sigmoid函数。Sigmoid函数在一定程度上近似于单位阶跃函数，同时单调可微，函数公式为

(4)

函数图像如图6所示。

图6 Sigmoid函数图像Fig.6 Image of Sigmoid function

因此，在线性回归模型基础上耦合Sigmoid函数，得到Logistic回归模型，可应用于二分类问题。

2.4 NewPGI模型构建

通过塑料大棚光谱特征分析和可分离性分析，最终选取Landsat8 影像7个OLI多光谱波段、1个TIR热红外波段以及NDVI，NDBaI和MNDWI这3个遥感指数，基于Logistic回归模型，构建NewPGI，计算公式为

(5)

(6)

式中：Xk(k=1,2,…,11)为上文选择的11个参数变量；ak(k=1,2,…,11)为11个参数相对应的回归系数，具体见表3，其中Sig.为显著性。如表3所示,Xk(k=1,2,…,11)(a1,…,a11)表示11个参数的Sig.值都小于0.05，通过置信度为95%的显著性检验。根据Logistic回归模型可知，z> 0，判定为正样本，即NewPGI> 0.5时，判定为塑料大棚像元；z< 0，判定为负样本，即NewPGI< 0.5时，判定为非塑料大棚像元。

表3 NewPGI中的参数Tab.3 Parameter variables in NewPGI

使用Logistic回归模型分析需满足以下条件： ①因变量为二分类变量； ②样本不能完全线性可分； ③样本数量不能太少(一般不少于200)。 本文设计的NewPGI基于Landsat8影像判定塑料大棚区域和非塑料大棚区域，所选样本区域的土地覆盖类型复杂，无法完全线性可分。在样本区域抽取385个像元，组建训练样本。基于SPSS19.0分析软件，采用Backward法筛选分类变量。模型系数及拟合度检验如表4所示，模型的X2=247.03，Sig.=0.000, Logistic回归模型具有显著性； 伪决定系数Cox&Snell和Nagelkerke值分别为0.414和0.697，模型的拟合度高，表明选取的11个参数(解释变量)对于塑料大棚(因变量)提取效果显著。

表4 Logistic回归模型系数的综合检验及拟合度检验Tab.4 Omnibus test of Logistic regression modelcoefficients and test of R-squared

3 结果与分析

3.1 样本区域塑料大棚信息提取与检验

基于NewPGI提取样本区域塑料大棚的分布信息(图7)。

(a) 样本区域NewPGI分布 (b) 设定阈值后NewPGI分布(c) 样本区域塑料大棚提取结果

图7 样本区域塑料大棚信息提取

Fig.7Plasticgreenhouseinformationextractioninthesamplearea

图7(a)为样本区域NewPGI分布图，相比其他土地覆盖类型，塑料大棚区域拥有更高的NewPGI值，与植被、人造地表等土地覆盖类型区别明显。通过设定最低阈值(0.5)和最高阈值(1)(图7(b))，去除其他土地覆盖类型，掩模得到样本区域二值化的塑料大棚分类结果(图7(c))。对比塑料大棚分类参考图(图3(a))，除去城市工业园区部分金属厂房的干扰，整体上具有高度的一致性。

为了评估NewPGI的分类精度，采用混淆矩阵，通过对比塑料大棚分类结果(图7(c))与塑料大棚分类参考图(图3(a))，引入Kappa系数，验证分类精度，结果如表5所示，NewPGI提取结果的总体精度为94.9%、塑料大棚用户精度为80.50%，Kappa系数为0.74，表明拟建的塑料大棚指数具有较高的分类精度。

表5 样本区域塑料大棚/非塑料大棚混淆矩阵Tab.5 PGs/No PGs confusion matrix of the sample area

3.2 研究区塑料大棚信息提取与检验

基于新构建的NewPGI提取整个研究区塑料大棚的分布信息。作为对比，同时运用Yang等[8]构建的适用于中国北方地区的遥感指数(plastic greenhouse index,PGI)提取常州市塑料大棚的分布信息，以验证2种遥感指数在太湖流域的适用性。PGI的公式为

(7)

该方法先利用NDVI和NDBI这2种遥感指数进行掩模处理，去除常绿植被和人造地表的影响； 然后利用B2，B3，B4，B5以及“B5-B2”波段组合，构建PGI，并确定PGI的下限和上限阈值分别为1.3和6.7。

图8为整个研究区塑料大棚的提取结果。

(a) NewPGI分布(b) 基于NewPGI塑料大棚分类

(c) PGI分布 (d) 基于PGI塑料大棚分类

图8 常州市塑料大棚信息提取

Fig.8PlasticgreenhouseinformationextractioninChangzhouCity

如图8(a)所示，新构建的NewPGI扩大了塑料大棚与人造地表、植被、水体等地表覆盖类型的光谱差异，相比于其他土地覆盖类型(绿色区域)，塑料大棚区域拥有更高的NewPGI值(红色区域)，与其他土地覆盖类型区分明显。相比之下， PGI对于区分常州市塑料大棚与其他土地覆盖的能力相对不足(图8(c))。

为了进一步量化NewPGI与现有遥感指数PGI监测塑料大棚的性能，本文基于Google Earth影像，随机抽取2 466个像元并归类为“塑料大棚像元”和“非塑料大棚像元”作为分类验证的样本点，用于验证2种指数的分类精度,分类精度评价如表6所示。

表6 研究区域塑料大棚/非塑料大棚混淆矩阵Tab.6 PGs/No PGs confusion matrix of the study area

由表6可知，与现有遥感指数PGI相比，NewPGI拥有更好的分类精度，总体分类精度为91.28%，Kappa系数为0.78，说明NewPGI指数在常州市拥有更好的适用性。

4 讨论

现有的塑料大棚遥感指数构建主要基于像元的光谱特征，通过选取针对塑料大棚较为敏感的波段，利用数学方法扩大塑料大棚与其他土地覆盖类型之间的光谱差异。但是，拥有复杂土地覆盖类型的太湖流域，存在混合像元且遥感“同物异谱，异物同谱”现象严重，单纯依靠少量的多光谱波段，无法有效扩大塑料大棚与其他地物类型的光谱差异。如若考虑加入更多遥感数据波段，简单的数学方程又难以确定遥感指数的数学形式，只能局限于利用少量的多光谱波段。本文基于塑料大棚可分离性分析，综合塑料大棚与各地物类型之间的最佳区分波段及遥感指数，通过Logistic回归模型确定了11个自变量(多光谱数据及遥感指数)及其最佳回归系数，并将塑料大棚比例大于0.5的像元归类为塑料大棚像元。因此，相较于已有的塑料大棚遥感指数，基于Logistic回归分析构建的NewPGI遥感指数能够在较复杂条件下有效提取塑料大棚的分布信息。

5 结论

运用遥感技术大范围监测塑料大棚的空间范围对于估算农作物产量以及预测塑料大棚对环境的影响至关重要。本文以地处太湖流域的常州市为例,得到如下结论：

1)在南方地区(如太湖流域)利用中等空间分辨率Landsat8影像进行塑料大棚提取，遥感“同物异谱，异物同谱”的现象严重，通过塑料大棚光谱特征分析，发现塑料大棚的光谱特征与裸地和休耕地、人造地表十分相似，区分难度比较大。

2)通过塑料大棚光谱的可分离性分析，选取Landsat8 影像7个OLI多光谱波段(B1—B7)、1个TIR1热红外波段(B10)以及3个遥感指数(NDVI，NDBaI和MNDWI)，共计11个参数，运用Logistic回归分析法，构建新的塑料大棚指数NewPGI，用于扩大塑料大棚与其他土地覆盖类型之间的光谱差异，使其与植被、裸地和休耕地、水体以及复杂的人造地表等土地覆盖类型分离。

3)精度验证结果表明，基于塑料大棚分类参考图，按照“逐个像元比对”原则，NewPGI在样本区域的Kappa系数为0.74，塑料大棚用户精度为80.50%，总体精度为94.9%； 基于Google Earth影像构建的验证样本，NewPGI在整个常州市的Kappa系数为0.78，塑料大棚用户精度为86.59%，总体精度为91.28%。与现有塑料大棚指数相比，本文构建的NewPGI指数更适用于复杂地形条件下的塑料大棚提取。