面向土地覆蓋變化信息局域精度預(yù)測的邏輯回歸協(xié)變量選用的直接策略與組合策略

張景雄，萬 月，梅瑩瑩

(1. 武漢大學(xué) 測繪學(xué)院，湖北 武漢 430079；2. 武漢大學(xué) 遙感信息工程學(xué)院，湖北 武漢 430079；3.地球空間信息技術(shù)協(xié)同創(chuàng)新中心，湖北 武漢 430079)

0 引 言

反映地表覆蓋類型的空間分布及其動(dòng)態(tài)變化特征的土地覆蓋（或稱地表覆蓋）信息產(chǎn)品日益積累，其在資源環(huán)境規(guī)劃和管理中具有越來越重要的作用。得益于時(shí)間維度信息的富集，地表覆蓋變化信息產(chǎn)品對于地理國情監(jiān)測、生態(tài)多樣性等具有重要意義[1]。精確的地表覆蓋變化信息有助于人們更好地理解地表過程，開展地球系統(tǒng)的科學(xué)研究，并為地理國情監(jiān)測和諸多地學(xué)應(yīng)用提供信息服務(wù)。地理國情監(jiān)測通過對地理環(huán)境進(jìn)行感知、統(tǒng)計(jì)和分析，服務(wù)于國計(jì)民生[2-3]。

有關(guān)土地覆蓋或類似的專題分類信息精度的研究也稱為真實(shí)性驗(yàn)證，在這一領(lǐng)域國內(nèi)外學(xué)者做出了許多重要工作[4]。1991年，Congalton[5]提出了分類精度的度量指標(biāo)PCC(percentage correctly classified)和誤差矩陣。為了將精度指標(biāo)應(yīng)用到具體區(qū)域，即從全部問題域到局部，McGwire[6]提出報(bào)告用戶精度時(shí)應(yīng)考慮具體小區(qū)域的精度。為計(jì)算用戶精度，De Wit[7]使用設(shè)置參考點(diǎn)的方式，以荷蘭為研究區(qū)域進(jìn)行真實(shí)性驗(yàn)證。為了避免設(shè)置參考點(diǎn)方式對精度計(jì)算的限制，一些學(xué)者探討了基于模型的精度空間分布估計(jì)方法。Steele等[8]使用克里金方法從樣本點(diǎn)空間內(nèi)插得到了誤分類概率分布圖。Smith等[9]使用邏輯回歸方法評價(jià)了地塊大小和鄰域異質(zhì)性對單一像素分類精度的影響。2004年，Van oort等[10]考慮了地表類別、鄰域異質(zhì)性、鄰域同質(zhì)性、地塊大小、熵和蔓延度等多種因素，分別計(jì)算了對應(yīng)的邏輯回歸模型參數(shù)，并評價(jià)了各因素對地表精度估計(jì)的影響程度。

土地覆蓋變化信息的精度研究也越來越受到重視[11]。Burnicki[12]于2011年提出了一種方法，通過綜合分析時(shí)間序列圖和景觀結(jié)構(gòu)來估計(jì)像素級的分類錯(cuò)誤概率空間分布。2016年，陳軍等[13]基于景觀指數(shù)，對精度驗(yàn)證使用的抽樣方法加以改進(jìn)，分別基于驗(yàn)證區(qū)域、區(qū)內(nèi)地類(土地景觀的空間異質(zhì)性)和抽樣格網(wǎng)三個(gè)層級，為大范圍精度驗(yàn)證提供了定量化的自適應(yīng)抽樣方法。Zhang and Mei (2016)[14]應(yīng)用邏輯回歸方法估計(jì)了湖北省洪湖市某子區(qū)2012～2013地表覆蓋變化的像素級的分類精度。

本文對地表覆蓋變化信息的局域精度預(yù)測方法作了進(jìn)一步的探究。前已述及，變化類別的像素級的正確分類概率可以通過邏輯回歸來計(jì)算,選用基于像素點(diǎn)所在3×3鄰域的類別的空間發(fā)生模式所計(jì)量的景觀指數(shù)和基于區(qū)內(nèi)地類鄰域的景觀指數(shù)（如面積）作為回歸建模的協(xié)變量。這里，景觀指數(shù)可以有兩種方式作為協(xié)變量，一是直接基于待驗(yàn)證的地表覆蓋變化地圖，二是組合利用變化前后兩時(shí)相對應(yīng)的地表覆蓋地圖。本研究通過實(shí)驗(yàn)，比較協(xié)變量選用的直接和組合策略的精度預(yù)測效果，評判兩種策略的預(yù)測效果和適用性。

論文第一部分說明研究區(qū)域和實(shí)驗(yàn)數(shù)據(jù)。第二部分描述邏輯回歸估計(jì)局域土地覆蓋精度的基本原理。第三部分分析實(shí)驗(yàn)結(jié)果。最后是結(jié)論和展望。

1 研究區(qū)域與研究數(shù)據(jù)

1.1 研究區(qū)域

本文選擇武漢市中心城區(qū)約10 km×10 km的區(qū)域作為研究區(qū)域，經(jīng)度范圍為114°1′40″E～114°5′40″E，緯度范圍為30°3′35″N～30°8′39″。選取的遙感影像為ETM+數(shù)據(jù)，像素尺寸為30 m×30 m，研究區(qū)域包括313×312個(gè)像素，時(shí)間為2011年12月和2013年8月兩期影像。

1.2 類別體系

考慮到研究區(qū)域的景觀特點(diǎn)，覆蓋類型大致包括植被（農(nóng)田、草地、樹木）、裸土（未開發(fā)土壤、河岸/湖岸灘地）、不透水表面（道路、橋梁、公共設(shè)施等建設(shè)用地、建筑物、居民地）和水體（水稻田、湖泊、河流、魚塘等）。因此，為選擇代表性和分離性較強(qiáng)的地表覆蓋類別，選用分類體系（不透水表面、水體、裸土、植被）在后文中分別以英文首字母I、W、S、V代替。

1.3 土地覆蓋變化信息提取

利用分類后處理的方法（post-classification）獲得土地覆蓋變化地圖，如圖1所示?？紤]到2011年11月～2013年8月研究區(qū)域可能發(fā)生的地理景觀變化，去除了不可能發(fā)生的(變化)類別和所占像素點(diǎn)過少的類別，比如水體->不透水表面，保留共9種類別（包括4種未變化類別和5種變化類別）。4類未變化類別為5種變化類別包括：IS：不透水表面->裸地；IV：不透水表面->植被；SI：裸地->不透水表面；SV：裸地->植被；VI：植被->不透水表面。

圖1 土地覆蓋與土地覆蓋變化圖：(a)-(b)時(shí)相1與時(shí)相2的土地覆蓋圖，(c)土地覆蓋變化Fig.1 Maps showing: (a)-(b) land cover at Time 1 and Time 2,respectively, and (c) land cover change

1.4 采樣方案與樣本數(shù)據(jù)采集

根據(jù)空間采樣原理[15]，在地表覆蓋變化圖中抽取樣本數(shù)據(jù)，包括訓(xùn)練樣本和測試樣本。前者用于建立邏輯回歸模型，后者用于比較兩種策略預(yù)測結(jié)果好壞。為了能準(zhǔn)確表達(dá)用戶精度，訓(xùn)練樣本和測試樣本的選取均須反映土地覆蓋變化類別的分布。本文采用基于分層隨機(jī)抽樣的原則，抽取符合要求的訓(xùn)練樣本：

式中，p表示期望達(dá)到的精度，Z為正態(tài)分布95%置信區(qū)間的標(biāo)準(zhǔn)差，Z=1.96?，F(xiàn)取p=77.12% p=0.7712，設(shè)定允許誤差E=0.05，由（1）式計(jì)算得到N為304。但在實(shí)際應(yīng)用中，每類別的樣本數(shù)不能過少，根據(jù)經(jīng)驗(yàn)設(shè)置最小類別數(shù)為20，調(diào)整后的訓(xùn)練樣本總數(shù)為322，根據(jù)同樣的抽樣原理，抽取同樣數(shù)量的測試樣本。表1給出訓(xùn)練樣本和測試樣本中各個(gè)變化類別的個(gè)數(shù)，圖2是訓(xùn)練樣本和測試樣本的分布。

表1 樣本點(diǎn)數(shù)量Tab.1 The number of sample points

圖2 訓(xùn)練和測試樣本分布圖Fig.2 Distribution map of training and testing sample points

2 方法

2.1 景觀指數(shù)計(jì)算

采用下述景觀指數(shù)：

CLASS: 像素所屬類別，使用二值變量來表示。變化類別共有9種類別，用8個(gè)二值變量表示，因?yàn)橛幸蝗哂?；如類別1可表示為[10 000 000]；單時(shí)相影像上分別有4種類別，各自用3個(gè)二值變量表示，則組合變量由6個(gè)二值變量表示。如若time1的類別為1，time2的類別為4，它們的組合可表示為[100 000]。

HOM(homogeneity): 局域同質(zhì)性，像素所在3×3鄰域內(nèi)與中心像素類別相同的像素個(gè)數(shù)。

HET(heterogeneity): 局域異質(zhì)性，即像素所在3×3鄰域所包含的不同類別數(shù)量。

LPS(log10(Patch size)): patch size為像素所在相連區(qū)域包含的元素個(gè)數(shù)（此處相連為4鄰域相連），再對其取對數(shù)。此處的計(jì)算方法利用了計(jì)算機(jī)形態(tài)學(xué)相關(guān)知識(shí)，首先從分類圖中分離出各類別的二值圖像，再使用函數(shù)bwconncomp分別計(jì)算各二值圖像相連區(qū)域的像素個(gè)數(shù)和各斑塊包含的像素地址列表，最后確定各元素所屬斑塊面積。本部分操作使用matlab完成。

鑒于本文研究的方法是面向局域的精度預(yù)測，上述景觀指數(shù)的量化的空間尺度包括單像素、像素所在3×3鄰域以及像素所在相連通的區(qū)域，類似于Zhang and Mei (2016)[14]。

2.2 邏輯回歸

建立分類精度與各景觀指數(shù)之間的邏輯回歸模型，從而得到逐像素的分類正確概率，可顯示為分類精度表面。回歸關(guān)系式可由式（2）表示，式中xi為各協(xié)變量，βi為對應(yīng)于各協(xié)變量的系數(shù)，β0為截距。本文選擇CLASS、HOM、HET、LPS等景觀指數(shù)作為協(xié)變量，P為回歸精度。

選取統(tǒng)計(jì)量Deviance來進(jìn)行模型選擇。它通過比較包含或不包含待檢驗(yàn)協(xié)變量的模型的對數(shù)似然函數(shù)，計(jì)算待評估協(xié)變量對模型的有用性。當(dāng)樣本量較大時(shí)，Deviance近似滿足自由度為上述待檢驗(yàn)變量個(gè)數(shù)（相對于不包含這些變量的基準(zhǔn)模型，這里的待檢驗(yàn)變量個(gè)數(shù)是待檢驗(yàn)?zāi)Ｐ退黾拥淖兞總€(gè)數(shù)）的χ2分布。

2.3 直接策略與組合策略模型介紹

為了探究土地覆蓋變化信息局部精度與景觀指數(shù)之間的定量關(guān)系，可考慮前述的兩種協(xié)變量選用策略。直接策略是直接利用所測試的地表覆蓋變化地圖的CLASS、HOM、HET、LPS作為協(xié)變量，見表2。組合策略利用兩個(gè)單時(shí)相地表覆蓋地圖的協(xié)變量的組合，即CLASS1、CLASS2、HOM1、HOM2、HET1、HET2、LPS1、LPS2，見表3。

表2 直接策略候選模型描述Tab.2 Descriptions of candidate models in the direct strategy

表3 組合策略候選模型描述Tab.3 Descriptions of candidate models in the composite strategy

在表2中模型0不包含任何協(xié)變量，只有截距。模型1a、1b、1c、1d中分別包含CLASS、HOM、HET、LPS。模型2a、2b、2c在模型1a的基礎(chǔ)上，分別增加協(xié)變量HOM、HET、LPS。模型3a，3b則相比于模型2a，分別增加HET、LPS。表3與表2類似。

2.4 模型驗(yàn)證

本文選擇以下評價(jià)指標(biāo)來比較模型的預(yù)測精度：均值ME(mean error), 平均絕對誤差MAE(mean absolute error)均方根誤差RMSE(root mea〈n square error)。i(xj)代表測試樣本點(diǎn)xj實(shí)際精度，(xj)為點(diǎn)xj估計(jì)值。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 邏輯回歸建模

根據(jù)表2、表3分別建立直接策略和組合策略的邏輯回歸模型；利用1.4節(jié)抽取的訓(xùn)練樣本估計(jì)回歸模型的系數(shù)，直接策略和組合策略邏輯回歸模型系數(shù)分別見表4、表5。表4、表5中的每一行分別對應(yīng)表2、表3中對應(yīng)項(xiàng)建立模型。以表4為例，第一行對應(yīng)模型0，β0為截距。第二行對應(yīng)模型1a, β1～β8為協(xié)變量CLASS中每一位二值變量的系數(shù)。

表4 直接策略候選模型系數(shù)Tab.4 Estimated regression coefficients for candidate models, the direct strategy

表5 組合策略候選模型系數(shù)Tab.5 Estimated regression coefficients for candidate models,the composite strategy

3.2 模型選擇

為了比較不同模型擬合的優(yōu)劣，根據(jù)（3）式選取統(tǒng)計(jì)量Deviance，采用似然比檢驗(yàn)法對模型進(jìn)行顯著性檢驗(yàn)，用于選擇含有最多有效協(xié)變量的模型。直接策略和組合策略邏輯回歸模型選擇的分析過程分別見表6、表7。

表6 直接策略模型選擇的x2檢驗(yàn)Tab.6 Chi-square tests for model selection,the direct strategy

表7 組合策略模型選擇的x2檢驗(yàn)Tab.7 Chi-square tests for model selection, the composite strategy

由表6可知，模型2a、2b、2c、2d中分別相對于模型0增加了CLASS、HOM、HET、LPS4個(gè)變量，這4個(gè)模型相對于基準(zhǔn)模型的改善作用可以分別用Deviance差值來衡量。由表6可以看出，CLASS的增加對改善模型的作用最大，且滿足自由度為8，α=0.01的χ2檢驗(yàn)。因此，增加CLASS作為優(yōu)選模型的第一個(gè)協(xié)變量。同理，相對于模型1a的比較，增加HOM作為第二個(gè)協(xié)變量。在模型3a、3b相對于模型2a的比較中，可看出變量HET、LPS的增加均對改善模型不具顯著作用。因此可得直接策略中協(xié)變量CLASS、HOM對改善模型有顯著作用，故直接策略中優(yōu)選模型為2a。

同樣，由表7可知，模型2a，2b，2c,2d中分別相對于模型0增加了CLASS1,2、HOM1,2、HET1,2、LPS1,2這幾個(gè)變量，這4個(gè)模型相對于基準(zhǔn)模型的改善作用分別可以用Deviance差值來衡量。可看出，CLASS1,2的增加對改善模型的作用最大，且CLASS1,2的增加滿足自由度為6，α=0.01的χ2檢驗(yàn)。因此增加CLASS1,2作為優(yōu)選模型的第一組協(xié)變量。同理，相對于模型1a的比較過程，增加HOM1,2作為第二組協(xié)變量。在模型3a、3b相對于模型2a的比較中，可看出變量HET1,2、LPS1,2的增加均對改善模型不具顯著作用。因此，組合策略協(xié)變量CLASS1,2、HOM1,2對改善模型有顯著作用，故得到組合策略優(yōu)選模型為2a。

3.3 策略比較

本部分根據(jù)2.4節(jié)方法對兩策略精度估計(jì)結(jié)果進(jìn)行定量評價(jià)。首先，分別根據(jù)直接策略2a模型和組合策略2a模型計(jì)算表示變化類型的分類精度的概率值。然后，基于322個(gè)測試樣本點(diǎn)，由式（4），（5），（6）計(jì)算精度預(yù)測值與樣本真值之間的平均誤差ME(mean error)、平均絕對誤差MAE(mean absolute error)、均方根誤差RMSE(root mean square error)，見表8。

表8 直接策略與組合策略模型2a預(yù)測精度比較Tab.8 Comparison accuracy estimation by model 2a's of the direct and composite strategies

RMSE、MAE、ME3個(gè)指標(biāo)值愈小，模型預(yù)測精度愈高。由表8可看出，兩策略下土地覆蓋變化信息的預(yù)測精度差別不大。由于直接策略相對于組合策略計(jì)算量較小，故推薦選擇直接策略。

4 結(jié)束語

本文比較了預(yù)測地表覆蓋變化精度的邏輯回歸方法中使用協(xié)變量的兩種策略，探究兩者在預(yù)測精度上是否有顯著差異。直接策略直接使用待驗(yàn)證的地表覆蓋變化地圖中提取的CLASS、HOM、HET、LPS等景觀指數(shù)作為協(xié)變量，組合策略以所涉及的兩個(gè)時(shí)相的土地覆蓋分類圖的景觀指數(shù)的組合（即CLASS1、CLASS2、HOM1、HOM2、HET1、HET2、LPS1、LPS2）作為協(xié)變量。以似然比檢驗(yàn)法對兩種策略可選模型進(jìn)行了顯著性檢驗(yàn)，分別獲得兩種策略下含有最多有效協(xié)變量的模型：前者包含CLASS和HOM，后者包括CLASS1、CLASS2、HOM1和HOM2。通過獨(dú)立的測試樣本，計(jì)算比較了上述最優(yōu)模型的ME、MAE和RMSE等評價(jià)指標(biāo)，發(fā)現(xiàn)兩策略下的模型預(yù)測能力差異較小。但是，直接基于土地覆蓋變化地圖進(jìn)行局域精度預(yù)測建模和計(jì)算的效率更高。因此，基于運(yùn)算效率的考慮，推薦在實(shí)際精度預(yù)測工作中使用直接策略。

本文的研究成果可以應(yīng)用于地表覆蓋變化圖的局域精度估計(jì)。后續(xù)研究可探討邏輯回歸與地質(zhì)統(tǒng)計(jì)方法的聯(lián)合使用，用以深入探討地理研究領(lǐng)域廣泛存在的空間-時(shí)間依賴特性，精化土地覆蓋變化信息精度預(yù)測方法。其次，應(yīng)研究發(fā)展面向土地覆蓋變化的精度驗(yàn)證的空間抽樣方法[13]，結(jié)合地圖精度-土地景觀格局的耦合關(guān)系，優(yōu)化抽樣方案，滿足精度驗(yàn)證的可靠性要求和特定的信息服務(wù)和應(yīng)用目的。另外，由于不同傳感器具有不同的空間分辨率，當(dāng)像素尺寸增大或減小時(shí)，景觀指數(shù)將隨之變化，最優(yōu)邏輯回歸模型及其預(yù)測值也將不同；探究邏輯回歸預(yù)測精度與所研究變化信息的尺度效應(yīng)將具有重要意義。