基于Logistic回歸分析的塑料大棚遙感指數構建

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

(1.南京信息工程大學地理科學學院，南京 210044； 2.中國科學院南京土壤研究所土壤與農業可持續發展國家重點實驗室，南京 210008)

0 引言

根據第三次全國農業普查結果，2016年末全國塑料大棚占地面積為98.1×104hm2，比2006年增長了111.0%[1]。塑料大棚大幅提高了蔬菜作物產量，但也加劇了“白色污染”。因此，準確監測塑料大棚的分布范圍，對區域農業的可持續發展和對環境影響的分析至關重要。

利用不同空間分辨率的遙感影像檢測設施菜地、地膜覆蓋等土地利用，已經引起越來越多的關注[2]?；诟呖臻g分辨率遙感影像(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影像作為提取樣本區域塑料大棚的參考影像?；贓NVI5.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個像元?；诒?所示的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個像元，組建訓練樣本?；赟PSS19.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指數更適用于復雜地形條件下的塑料大棚提取。