基于Logistic回歸分析的塑料大棚遙感指數(shù)構(gòu)建

陳俊， 沈潤平， 李博倫， 遆超普， 顏曉元， 周旻悅， 王紹武

(1.南京信息工程大學(xué)地理科學(xué)學(xué)院，南京 210044； 2.中國科學(xué)院南京土壤研究所土壤與農(nóng)業(yè)可持續(xù)發(fā)展國家重點實驗室，南京 210008)

0 引言

根據(jù)第三次全國農(nóng)業(yè)普查結(jié)果，2016年末全國塑料大棚占地面積為98.1×104hm2，比2006年增長了111.0%[1]。塑料大棚大幅提高了蔬菜作物產(chǎn)量，但也加劇了“白色污染”。因此，準(zhǔn)確監(jiān)測塑料大棚的分布范圍，對區(qū)域農(nóng)業(yè)的可持續(xù)發(fā)展和對環(huán)境影響的分析至關(guān)重要。

利用不同空間分辨率的遙感影像檢測設(shè)施菜地、地膜覆蓋等土地利用，已經(jīng)引起越來越多的關(guān)注[2]?；诟呖臻g分辨率遙感影像(0.5～2 m)進(jìn)行塑料大棚提取是目前常用的檢測方法，例如，Agüera等[3]利用QuickBird影像，確定塑料大棚提取的最佳波段組合為綠光、藍(lán)光和近紅外波段； Aguilar等[4]根據(jù)0.5 m空間分辨率GeoEye-1和WorldView-2影像，提出一種面向?qū)ο蟮乃芰洗笈锓诸惙椒āａ槍Υ竺娣e的塑料大棚提取，國內(nèi)外學(xué)者探討了中等空間分辨率影像(2～30 m)用于塑料大棚提取的可行性[5]。例如，國內(nèi)最早進(jìn)行嘗試的Zhao等[6]基于Landsat TM影像，利用指數(shù)方法在山東省進(jìn)行塑料大棚的制圖； Novelli等[7]比較了Sentinel-2 MSI和Landsat8 OLI影像的溫室檢測性能； Yang等[8]基于Landsat ETM+影像，利用指數(shù)法提取了山東省濰坊市塑料大棚面積。

目前，我國塑料大棚遙感提取的研究主要位于塑料大棚集中分布的北方地區(qū)。相比之下，太湖流域河網(wǎng)密布，景觀破碎化程度高，土地覆蓋類型復(fù)雜，塑料大棚分布相對離散。而現(xiàn)有塑料大棚指數(shù)計算時所需的波段數(shù)量較少，難以滿足復(fù)雜地表覆蓋下準(zhǔn)確識別塑料大棚的要求。Logistic回歸分析是一種根據(jù)單個或多個自變量，分析和預(yù)測因變量的多元分析方法，是目前常用的處理分類因變量的統(tǒng)計分類模型[9]。因此，本文基于Landsat8影像，結(jié)合4種遙感指數(shù)，即歸一化植被指數(shù)(normalized difference vegetation index，NDVI)、歸一化建筑指數(shù)(normalized difference building index，NDBI)、歸一化裸土指數(shù)(normalised difference bareness index，NDBaI)以及改進(jìn)的歸一化水體指數(shù)(modified normalized difference water index，MNDWI)，經(jīng)可分離性篩選后，運用Logistic回歸模型，設(shè)計新塑料大棚遙感指數(shù)(new plastic greenhouse index, NewPGI)，以期更好地提取塑料大棚的分布信息。

1 研究區(qū)概況和數(shù)據(jù)源

1.1 研究區(qū)概況

本文選取的研究區(qū)為位于太湖流域西北部的常州市，如圖1所示，地理坐標(biāo)為N31°49′，E119°50′，行政區(qū)劃包括新北區(qū)、鐘樓區(qū)、天寧區(qū)、武進(jìn)區(qū)、金壇區(qū)和溧陽市。研究區(qū)屬于亞熱帶季風(fēng)氣候區(qū)，盛夏高溫多雨，嚴(yán)冬寒冷干燥，以塑料大棚為主導(dǎo)的設(shè)施菜地種植模式，成為該地區(qū)農(nóng)業(yè)生產(chǎn)最主要的發(fā)展方向。經(jīng)過實地調(diào)研確定該地區(qū)蔬菜大棚的建材以透明塑料薄膜為主，根據(jù)Yang等[8]研究，選在作物生長季節(jié)提取塑料大棚信息。

(a) 研究區(qū)相對位置示意圖 (b) Landsat8 B5(R),B4(G),B3(B)假彩色合成影像

圖1 研究區(qū)位置及遙感影像

Fig.1Locationandimageofstudyarea

1.2 數(shù)據(jù)源及其預(yù)處理

采用覆蓋整個研究區(qū)的2014年3月16日Landsat8影像(空間分辨率為30 m)提取常州市塑料大棚分布信息。利用ENVI 軟件中FLAASH模塊對OLI影像進(jìn)行輻射校準(zhǔn)和大氣校正預(yù)處理，以消除大氣影響； 對TIRS影像的熱紅外波段值進(jìn)行歸一化處理，使其與OLI影像反射率范圍一致，以便后續(xù)處理。借助資源三號(ZY-3)影像，選取常州市塑料大棚相對集中區(qū)域作為樣本區(qū)域，用于構(gòu)建和驗證塑料大棚指數(shù)(圖2(a)，大小為9 km × 9 km)。采用2014年4月6日ZY-3影像作為提取樣本區(qū)域塑料大棚的參考影像?；贓NVI5.3軟件，利用ZY-3全色影像(圖2(b))及30 m空間分辨率數(shù)字高程模型(digital elevation model,DEM)數(shù)據(jù)，對ZY-3多光譜影像(圖2(c))進(jìn)行正射校正，以糾正因系統(tǒng)因素或地形因子引起的幾何畸變。對ZY-3多光譜及全色影像采用Gram-Schmidt光譜銳化法進(jìn)行影像融合，融合后影像空間分辨率為2 m。

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

圖2 樣本區(qū)域遙感影像數(shù)據(jù)

Fig.2Satellitedatainthesamplearea

1.3 驗證數(shù)據(jù)處理

1.3.1 樣本區(qū)域分類參考圖制作

利用eCognition影像分析軟件，基于地理空間信息和專家知識，采用面向?qū)ο蠓诸惙椒▽θ诤虾蟮腪Y-3影像分類，并通過目視解譯修正分類結(jié)果作為樣本區(qū)域分類參考(圖3(a))。將樣本區(qū)域典型土地覆蓋類型分為5類，包括塑料大棚、人造地表、裸地和休耕地、水體以及植被，得到樣本區(qū)域土地覆蓋分類結(jié)果(圖3(b))。

(a) 塑料大棚分類參考(b) 土地覆蓋分類

圖3 樣本區(qū)域分類參考

Fig.3Referenceofclassificationinthesamplearea

將塑料大棚類型所占像元比例超過50%的像元作為塑料大棚像元，通過IDL語言按照30 m×30 m像元尺寸,對圖3(b)進(jìn)行重采樣，生成塑料大棚分類參考圖，空間分辨率為30 m，用于驗證基于Landsat8影像的塑料大棚提取結(jié)果。

1.3.2 研究區(qū)驗證樣本

基于Google Earth影像進(jìn)行整個研究區(qū)塑料大棚分類的精度驗證。隨機(jī)選取2014年3月16日常州市塑料大棚、人造地表、裸地和休耕地、水體以及植被等土地覆蓋類型(圖4)，共抽取2 466個參考像元并歸類為“塑料大棚像元”和“非塑料大棚像元”，生成感興趣區(qū)域(region of interest,ROI)，構(gòu)建塑料大棚分類的驗證樣本，用于驗證研究區(qū)塑料大棚的分類精度。

(a) 研究區(qū)驗證樣本 (b) 驗證樣本示例1(c) 驗證樣本示例2

圖4 研究區(qū)驗證樣本

Fig.4Validationdataofstudyarea

2 研究方法

2.1 塑料大棚光譜特征分析

識別塑料大棚的光譜特征對于設(shè)計合理的塑料大棚指數(shù)至關(guān)重要。本文基于Landsat8影像和樣本區(qū)域土地覆蓋分類，使用樣本量計算工具(http: //fluidsurveys.com/university/survey-sample-size-calculator)，按照95%置信水平，隨機(jī)抽取5種地物類型共385個像元作為研究樣本，每種地物類型各選取77個像元?；诒?所示的Landsat8影像7個OLI多光譜波段(B1—B7)、1個TIR熱紅外波段(B10)以及常用于反映植被、人造地表、裸土和水體信息的4個遙感指數(shù)，即NDVI，NDBI，NDBaI和MNDWI，得出不同土地覆蓋類型光譜的平均值(圖5)。

表1 基于Landsat8影像的光譜數(shù)據(jù)及多種遙感指數(shù)Tab.1 Spectral data based on Landsat8 image and various remote sensing indexes

圖5 不同土地覆蓋類型光譜曲線(平均值)Fig.5 Spectral curves of differentland cover types (mean values)

塑料大棚作為一種人造設(shè)施，因在農(nóng)作物之上覆蓋一層白色塑料薄膜，削弱了作物的植被信息，使其同時兼具覆蓋少量植被的土壤和人造地表雙重的光譜特征。根據(jù)圖5不難發(fā)現(xiàn)，塑料大棚與裸地和休耕地、人造地表的光譜特征十分相似，區(qū)分難度比較大，佐證了復(fù)雜土地覆蓋類型下，塑料大棚遙感提取的影響因子眾多，需要考慮加入更多的波段信息的設(shè)想。

2.2 塑料大棚光譜可分離性分析

塑料大棚指數(shù)的設(shè)計應(yīng)該考慮最大限度地區(qū)分塑料大棚與其他土地覆蓋類別。為實現(xiàn)這一目標(biāo)，本文基于Kaufman等[13]研究，采用可分離性指標(biāo)M，逐一比較塑料大棚與其他土地覆蓋類型(人造地表、裸地和休耕地、植被以及水體)的分離度。該指標(biāo)被定義為2種地物類型光譜曲線平均值μ之間的差異，按照標(biāo)準(zhǔn)差δ之和進(jìn)行歸一化處理，即

(1)

式中：μPG為塑料大棚反射率均值；μ1～μ4分別為人造地表、裸地和休耕地、植被以及水體的反射率均值；δPG為塑料大棚光譜標(biāo)準(zhǔn)差；δ1～δ4分別為人造地表、裸地和休耕地、植被以及水體的光譜標(biāo)準(zhǔn)差。M≥1，表示分離性好；M<1，表示分離性較差。

塑料大棚與其他4種土地覆蓋類型之間的可分離性指標(biāo)M如表2所示。根據(jù)分離標(biāo)準(zhǔn)(M≥1)，確定區(qū)分塑料大棚與人造地表的最佳波段是B3和B5； 區(qū)分塑料大棚與裸地和休耕地的最佳波段是B1和B2； 區(qū)分塑料大棚與植被的最佳波段和指數(shù)是B1，B2，B3，B4，B7，NDVI和MNDWI； 塑料大棚與水體可通過除NDBI外的其他波段和指數(shù)進(jìn)行區(qū)分。綜合考慮塑料大棚與其他地物類型區(qū)分度，加入更多波段及指數(shù)以涵蓋不同土地覆蓋類型，最終選定B1—B7，B10，NDVI，NDBaI和MNDWI共計11個參數(shù)，用于構(gòu)建NewPGI。

表2 塑料大棚與典型土地覆蓋類型可分離性指標(biāo)MTab.2 Separation degree of plastic greenhousesfrom typical land cover types

2.3 Logistic回歸分析

本文引入的Logistic回歸模型屬于一種機(jī)器學(xué)習(xí)分類算法。進(jìn)行統(tǒng)計分析時，自變量可以是連續(xù)變量，也可以是離散變量，且不需要呈正態(tài)分布，增強(qiáng)了模型的應(yīng)用范圍和靈活性[14]，其最大優(yōu)勢是能夠?qū)Ρ姸嘤绊懸蛩剡M(jìn)行擬合分析，通過機(jī)器學(xué)習(xí)確定眾多自變量最佳回歸系數(shù)。參考王濟(jì)川等[15]和Cox[16]的研究，采用二分類因變量Logistic回歸模型。

作為一種概率型的非線性回歸模型，Logistic回歸是一種研究二分類觀察結(jié)果y與諸多影響因素Xk(k=1,2,…，n)之間關(guān)系的多變量分析方法。當(dāng)輸入測試樣本數(shù)據(jù)時，Logistic回歸分類器就是一組權(quán)值(a0,a1,…,an)，這組權(quán)值與測試數(shù)據(jù)按照線性加和得到預(yù)測值z，即

(2)

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

Logistic回歸模型的關(guān)鍵在于通過預(yù)測值z定義判定邊界，以此確定樣本的類型。通過單位階躍函數(shù)，可表示預(yù)測值z與二分類觀察結(jié)果y之間的關(guān)系，即

(3)

z為模型的判定邊界，用于確定哪些樣本是正樣本，哪些為負(fù)樣本。當(dāng)z> 0時，函數(shù)值為1，判定為正樣本； 當(dāng)z< 0時，函數(shù)值為0，判定為負(fù)樣本； 當(dāng)z= 0時，函數(shù)值為0.5，表示樣本為正樣本或負(fù)樣本概率相同，可任意判斷。但是單位階躍函數(shù)是分段式非連續(xù)性函數(shù)，無法應(yīng)用于實際問題。因此，需要引入一個連續(xù)性函數(shù)——Sigmoid函數(shù)。Sigmoid函數(shù)在一定程度上近似于單位階躍函數(shù)，同時單調(diào)可微，函數(shù)公式為

(4)

函數(shù)圖像如圖6所示。

圖6 Sigmoid函數(shù)圖像Fig.6 Image of Sigmoid function

因此，在線性回歸模型基礎(chǔ)上耦合Sigmoid函數(shù)，得到Logistic回歸模型，可應(yīng)用于二分類問題。

2.4 NewPGI模型構(gòu)建

通過塑料大棚光譜特征分析和可分離性分析，最終選取Landsat8 影像7個OLI多光譜波段、1個TIR熱紅外波段以及NDVI，NDBaI和MNDWI這3個遙感指數(shù)，基于Logistic回歸模型，構(gòu)建NewPGI，計算公式為

(5)

(6)

式中：Xk(k=1,2,…,11)為上文選擇的11個參數(shù)變量；ak(k=1,2,…,11)為11個參數(shù)相對應(yīng)的回歸系數(shù)，具體見表3，其中Sig.為顯著性。如表3所示,Xk(k=1,2,…,11)(a1,…,a11)表示11個參數(shù)的Sig.值都小于0.05，通過置信度為95%的顯著性檢驗。根據(jù)Logistic回歸模型可知，z> 0，判定為正樣本，即NewPGI> 0.5時，判定為塑料大棚像元；z< 0，判定為負(fù)樣本，即NewPGI< 0.5時，判定為非塑料大棚像元。

表3 NewPGI中的參數(shù)Tab.3 Parameter variables in NewPGI

使用Logistic回歸模型分析需滿足以下條件： ①因變量為二分類變量； ②樣本不能完全線性可分； ③樣本數(shù)量不能太少(一般不少于200)。 本文設(shè)計的NewPGI基于Landsat8影像判定塑料大棚區(qū)域和非塑料大棚區(qū)域，所選樣本區(qū)域的土地覆蓋類型復(fù)雜，無法完全線性可分。在樣本區(qū)域抽取385個像元，組建訓(xùn)練樣本。基于SPSS19.0分析軟件，采用Backward法篩選分類變量。模型系數(shù)及擬合度檢驗如表4所示，模型的X2=247.03，Sig.=0.000, Logistic回歸模型具有顯著性； 偽決定系數(shù)Cox&Snell和Nagelkerke值分別為0.414和0.697，模型的擬合度高，表明選取的11個參數(shù)(解釋變量)對于塑料大棚(因變量)提取效果顯著。

表4 Logistic回歸模型系數(shù)的綜合檢驗及擬合度檢驗Tab.4 Omnibus test of Logistic regression modelcoefficients and test of R-squared

3 結(jié)果與分析

3.1 樣本區(qū)域塑料大棚信息提取與檢驗

基于NewPGI提取樣本區(qū)域塑料大棚的分布信息(圖7)。

(a) 樣本區(qū)域NewPGI分布 (b) 設(shè)定閾值后NewPGI分布(c) 樣本區(qū)域塑料大棚提取結(jié)果

圖7 樣本區(qū)域塑料大棚信息提取

Fig.7Plasticgreenhouseinformationextractioninthesamplearea

圖7(a)為樣本區(qū)域NewPGI分布圖，相比其他土地覆蓋類型，塑料大棚區(qū)域擁有更高的NewPGI值，與植被、人造地表等土地覆蓋類型區(qū)別明顯。通過設(shè)定最低閾值(0.5)和最高閾值(1)(圖7(b))，去除其他土地覆蓋類型，掩模得到樣本區(qū)域二值化的塑料大棚分類結(jié)果(圖7(c))。對比塑料大棚分類參考圖(圖3(a))，除去城市工業(yè)園區(qū)部分金屬廠房的干擾，整體上具有高度的一致性。

為了評估NewPGI的分類精度，采用混淆矩陣，通過對比塑料大棚分類結(jié)果(圖7(c))與塑料大棚分類參考圖(圖3(a))，引入Kappa系數(shù)，驗證分類精度，結(jié)果如表5所示，NewPGI提取結(jié)果的總體精度為94.9%、塑料大棚用戶精度為80.50%，Kappa系數(shù)為0.74，表明擬建的塑料大棚指數(shù)具有較高的分類精度。

表5 樣本區(qū)域塑料大棚/非塑料大棚混淆矩陣Tab.5 PGs/No PGs confusion matrix of the sample area

3.2 研究區(qū)塑料大棚信息提取與檢驗

基于新構(gòu)建的NewPGI提取整個研究區(qū)塑料大棚的分布信息。作為對比，同時運用Yang等[8]構(gòu)建的適用于中國北方地區(qū)的遙感指數(shù)(plastic greenhouse index,PGI)提取常州市塑料大棚的分布信息，以驗證2種遙感指數(shù)在太湖流域的適用性。PGI的公式為

(7)

該方法先利用NDVI和NDBI這2種遙感指數(shù)進(jìn)行掩模處理，去除常綠植被和人造地表的影響； 然后利用B2，B3，B4，B5以及“B5-B2”波段組合，構(gòu)建PGI，并確定PGI的下限和上限閾值分別為1.3和6.7。

圖8為整個研究區(qū)塑料大棚的提取結(jié)果。

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

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

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

Fig.8PlasticgreenhouseinformationextractioninChangzhouCity

如圖8(a)所示，新構(gòu)建的NewPGI擴(kuò)大了塑料大棚與人造地表、植被、水體等地表覆蓋類型的光譜差異，相比于其他土地覆蓋類型(綠色區(qū)域)，塑料大棚區(qū)域擁有更高的NewPGI值(紅色區(qū)域)，與其他土地覆蓋類型區(qū)分明顯。相比之下， PGI對于區(qū)分常州市塑料大棚與其他土地覆蓋的能力相對不足(圖8(c))。

為了進(jìn)一步量化NewPGI與現(xiàn)有遙感指數(shù)PGI監(jiān)測塑料大棚的性能，本文基于Google Earth影像，隨機(jī)抽取2 466個像元并歸類為“塑料大棚像元”和“非塑料大棚像元”作為分類驗證的樣本點，用于驗證2種指數(shù)的分類精度,分類精度評價如表6所示。

表6 研究區(qū)域塑料大棚/非塑料大棚混淆矩陣Tab.6 PGs/No PGs confusion matrix of the study area

由表6可知，與現(xiàn)有遙感指數(shù)PGI相比，NewPGI擁有更好的分類精度，總體分類精度為91.28%，Kappa系數(shù)為0.78，說明NewPGI指數(shù)在常州市擁有更好的適用性。

4 討論

現(xiàn)有的塑料大棚遙感指數(shù)構(gòu)建主要基于像元的光譜特征，通過選取針對塑料大棚較為敏感的波段，利用數(shù)學(xué)方法擴(kuò)大塑料大棚與其他土地覆蓋類型之間的光譜差異。但是，擁有復(fù)雜土地覆蓋類型的太湖流域，存在混合像元且遙感“同物異譜，異物同譜”現(xiàn)象嚴(yán)重，單純依靠少量的多光譜波段，無法有效擴(kuò)大塑料大棚與其他地物類型的光譜差異。如若考慮加入更多遙感數(shù)據(jù)波段，簡單的數(shù)學(xué)方程又難以確定遙感指數(shù)的數(shù)學(xué)形式，只能局限于利用少量的多光譜波段。本文基于塑料大棚可分離性分析，綜合塑料大棚與各地物類型之間的最佳區(qū)分波段及遙感指數(shù)，通過Logistic回歸模型確定了11個自變量(多光譜數(shù)據(jù)及遙感指數(shù))及其最佳回歸系數(shù)，并將塑料大棚比例大于0.5的像元歸類為塑料大棚像元。因此，相較于已有的塑料大棚遙感指數(shù)，基于Logistic回歸分析構(gòu)建的NewPGI遙感指數(shù)能夠在較復(fù)雜條件下有效提取塑料大棚的分布信息。

5 結(jié)論

運用遙感技術(shù)大范圍監(jiān)測塑料大棚的空間范圍對于估算農(nóng)作物產(chǎn)量以及預(yù)測塑料大棚對環(huán)境的影響至關(guān)重要。本文以地處太湖流域的常州市為例,得到如下結(jié)論：

1)在南方地區(qū)(如太湖流域)利用中等空間分辨率Landsat8影像進(jìn)行塑料大棚提取，遙感“同物異譜，異物同譜”的現(xiàn)象嚴(yán)重，通過塑料大棚光譜特征分析，發(fā)現(xiàn)塑料大棚的光譜特征與裸地和休耕地、人造地表十分相似，區(qū)分難度比較大。

2)通過塑料大棚光譜的可分離性分析，選取Landsat8 影像7個OLI多光譜波段(B1—B7)、1個TIR1熱紅外波段(B10)以及3個遙感指數(shù)(NDVI，NDBaI和MNDWI)，共計11個參數(shù)，運用Logistic回歸分析法，構(gòu)建新的塑料大棚指數(shù)NewPGI，用于擴(kuò)大塑料大棚與其他土地覆蓋類型之間的光譜差異，使其與植被、裸地和休耕地、水體以及復(fù)雜的人造地表等土地覆蓋類型分離。

3)精度驗證結(jié)果表明，基于塑料大棚分類參考圖，按照“逐個像元比對”原則，NewPGI在樣本區(qū)域的Kappa系數(shù)為0.74，塑料大棚用戶精度為80.50%，總體精度為94.9%； 基于Google Earth影像構(gòu)建的驗證樣本，NewPGI在整個常州市的Kappa系數(shù)為0.78，塑料大棚用戶精度為86.59%，總體精度為91.28%。與現(xiàn)有塑料大棚指數(shù)相比，本文構(gòu)建的NewPGI指數(shù)更適用于復(fù)雜地形條件下的塑料大棚提取。