郭含茹，張茂震，徐麗華，袁振花，陳田閣

（1.浙江農(nóng)林大學(xué) 浙江省森林生態(tài)系統(tǒng)碳循環(huán)與固碳減排重點(diǎn)實(shí)驗(yàn)室，浙江 臨安 311300；2.浙江農(nóng)林大學(xué) 環(huán)境與資源學(xué)院，浙江 臨安311300；3.河南省濮陽市華龍區(qū)科學(xué)技術(shù)局，河南濮陽457001）

基于地理加權(quán)回歸的區(qū)域森林碳儲(chǔ)量估計(jì)

郭含茹1,2，張茂震1,2，徐麗華1,2，袁振花1,2，陳田閣3

（1.浙江農(nóng)林大學(xué) 浙江省森林生態(tài)系統(tǒng)碳循環(huán)與固碳減排重點(diǎn)實(shí)驗(yàn)室，浙江 臨安 311300；2.浙江農(nóng)林大學(xué) 環(huán)境與資源學(xué)院，浙江 臨安311300；3.河南省濮陽市華龍區(qū)科學(xué)技術(shù)局，河南濮陽457001）

森林碳儲(chǔ)量與其調(diào)查因子之間的關(guān)系及其空間相關(guān)性特征是區(qū)域森林碳儲(chǔ)量及其分布估計(jì)模型建立的基礎(chǔ)，然而某一調(diào)查因子在不同空間位置對森林碳儲(chǔ)量估計(jì)的影響程度存在差異。地理加權(quán)回歸方法考慮了調(diào)查因子作用的空間異質(zhì)性，進(jìn)行局部回歸估計(jì)。以浙江省臺(tái)州市仙居縣森林資源二類調(diào)查樣地實(shí)測數(shù)據(jù)為數(shù)據(jù)源，利用地理加權(quán)回歸方法結(jié)合陸地探測衛(wèi)星系統(tǒng)Landsat TM影像數(shù)據(jù)進(jìn)行森林碳儲(chǔ)量及其密度的分布估計(jì)，并檢驗(yàn)地理和海拔加權(quán)回歸在地勢平緩的區(qū)域是否可用。結(jié)果表明：①仙居縣森林地上部分總碳儲(chǔ)量3.132×106Mg，與樣地實(shí)測統(tǒng)計(jì)得到的全縣碳儲(chǔ)量3.192×106Mg相差1.880%；地理加權(quán)回歸模型估計(jì)結(jié)果與實(shí)際碳密度分布情況相一致，研究區(qū)內(nèi)碳密度的取值范圍為0～89.964 Mg·hm-2，保留了70%以上的空間異質(zhì)性特征；基于地理加權(quán)回歸的森林地上部分碳儲(chǔ)量估計(jì)方法是有效的，地理加權(quán)回歸在區(qū)域碳儲(chǔ)量方面估計(jì)結(jié)果合理且精度較高。②在地勢較為平緩的地區(qū)，海拔對植被的影響不顯著，地理和海拔加權(quán)回歸并不適用；若將海拔作為解釋變量加入建模，能夠提高估計(jì)精度，但存在多重共線性問題。圖6表9參26

森林生態(tài)學(xué)；森林碳儲(chǔ)量;空間異質(zhì)性;地理加權(quán)回歸;地理和海拔加權(quán)回歸

森林是陸地生態(tài)系統(tǒng)中最大的碳庫，在全球碳循環(huán)過程中占據(jù)重要地位。Phillips等［1］對熱帶森林生物量的研究表明：熱帶森林生物量的較小波動(dòng)即可引起全球碳循環(huán)的顯著變化?？梢?，森林生態(tài)系統(tǒng)碳儲(chǔ)量和動(dòng)態(tài)的研究直接關(guān)系到降低大氣中二氧化碳濃度、減緩全球氣候變暖的重大議題。中國作為一個(gè)發(fā)展中國家，能源的消耗和二氧化碳的排放量相對較大，且中國在全球固碳減排行動(dòng)中充當(dāng)著重要的角色。精準(zhǔn)地估算中國森林植被碳儲(chǔ)量，是進(jìn)行碳收支和碳交易評估的前提［2］。區(qū)域森林碳儲(chǔ)量的估測普遍采用的方法是由直接或間接測定森林植被的生產(chǎn)量與生物現(xiàn)存量再乘以生物量中碳元素的含量，推算得到［3］。森林植被的含碳率一般為0.45～0.50［4］。本研究采用政府間氣候變化專門委員會(huì)（IPCC）提供的默認(rèn)參數(shù)0.50計(jì)［5］。目前，國內(nèi)外應(yīng)用最為廣泛的生物量估測方法有3種：樣地清查法、渦度相關(guān)法、應(yīng)用遙感等新技術(shù)的模型模擬法［6］。近年來，隨著遙感、地理信息系統(tǒng)（GIS）等新技術(shù)在森林生物量估測方面的應(yīng)用，第3種方法在森林生物量估算研究中應(yīng)用最為廣泛?；贕IS的地統(tǒng)計(jì)學(xué)分析為研究森林調(diào)查因子的空間相關(guān)性和依賴性現(xiàn)象提供了科學(xué)方法［7］。李明詩等［8］的研究表明：利用光譜、地形及紋理特征對森林生物量進(jìn)行建?？梢蕴岣呱稚锪康墓罍y精度，且不同樹種與不同的特征具有強(qiáng)相關(guān)性，特征貢獻(xiàn)為0.044～0.205。常規(guī)基于普通最小二乘法（OLS）的生物量模型估測，忽略了調(diào)查因子作用的空間異質(zhì)性。而實(shí)際上，某一調(diào)查因子在不同空間位置對生物量估測的影響程度是有差異的。在空間建模領(lǐng)域，一些統(tǒng)計(jì)技術(shù)已經(jīng)發(fā)展到模擬空間變量之間的空間關(guān)系。在探索空間關(guān)系的異質(zhì)性方面，最有力的工具之一就是地理加權(quán)回歸（GWR）［9］。地理加權(quán)回歸方法是傳統(tǒng)回歸分析方法的擴(kuò)展，進(jìn)行局部而非全局的參數(shù)估計(jì)，GWR方法在數(shù)據(jù)集的所有位置點(diǎn)的參數(shù)進(jìn)行估計(jì)時(shí)考慮了非平穩(wěn)性的關(guān)系，針對每一個(gè)坐標(biāo)位置點(diǎn)都有相對應(yīng)的參數(shù)。GWR方法被廣泛地運(yùn)用到各領(lǐng)域研究中，且均得到較好的研究結(jié)果。邵一希等［10］將地理加權(quán)回歸方法應(yīng)用在區(qū)域土地利用格局模擬建模中，與OLS的Logistic回歸模型進(jìn)行比較，得出GWR模型可以得到更加精確的擬合優(yōu)度和擬合準(zhǔn)確率。郭龍等［11］的研究表明：在土壤屬性空間預(yù)測領(lǐng)域，GWR方法比協(xié)同克里格方法具有更高的預(yù)測精度。然而，林業(yè)遙感應(yīng)用領(lǐng)域，GWR建模是一個(gè)相對未開發(fā)的區(qū)域［9］。Propastin［12］針對多山的熱帶雨林地區(qū)地上生物量的估測進(jìn)行了地理加權(quán)模型的改進(jìn)研究，考慮到海拔對森林生物量的分布影響，將高度影響作為第三維度加入到GWR加權(quán)矩陣來擴(kuò)展GWR模型，得到地理和海拔加權(quán)回歸（GAWR）模型。結(jié)果表明：針對海拔高度影響較大的數(shù)據(jù)集，GAWR模型比GWR模型效果更好［9］。但對于不同于多山的熱帶雨林的其他地形地貌植被類型等來說，在地上生物量的估測研究中，尚未有明確的結(jié)論證明地理加權(quán)回歸方法是否也可以獲得較高的估測精度；如果應(yīng)用數(shù)據(jù)來自于不太明顯的海拔梯度的區(qū)域，第三維度的加入是否會(huì)有更好的效果？或者使用不同模型與GWR模型的比較研究，也是值得深入研究的。本研究以浙江省臺(tái)州市仙居縣為研究區(qū)，嘗試地理加權(quán)回歸方法估計(jì)區(qū)域森林碳儲(chǔ)量，利用實(shí)測地上部分碳儲(chǔ)量數(shù)據(jù)和美國陸地探測衛(wèi)星系統(tǒng)專題繪圖儀（Landsat TM）數(shù)據(jù)建立局部回歸模型，并檢驗(yàn)在海拔較為平緩的地區(qū)，地理和海拔加權(quán)回歸模型是否可用，以期提高森林碳儲(chǔ)量的估計(jì)精度，對進(jìn)一步量化區(qū)域森林碳匯能力具有重要現(xiàn)實(shí)意義。

1 材料和數(shù)據(jù)

1.1 研究區(qū)概況

仙居縣位于臺(tái)州市西部（28°28′～28°59′N，120°17′～120°55′E）（圖1），總面積為20.131 8×104hm2。全縣地形從外向內(nèi)傾斜，略向東傾，其間有大小不等、錯(cuò)落相間的谷地和盆地。屬于亞熱帶季風(fēng)氣候，雨水充沛。仙居縣是浙江省重點(diǎn)林區(qū)縣之一，2011年全縣森林覆蓋率為77.9%，最高海拔1 384.4 m。研究區(qū)中部狹長地帶為平緩地帶和城鎮(zhèn)聚居區(qū)，西南部是國家級森林公園。研究區(qū)內(nèi)杉木Cunninghamia lanceolata，馬尾松Pinus massonania和闊葉樹占主要地位。

圖1 研究區(qū)位置和樣地分布示意圖Figure 1 Location of the study area and the plots distribution

1.2 數(shù)據(jù)準(zhǔn)備

1.2.1 地面調(diào)查數(shù)據(jù) 研究所使用的森林資源數(shù)據(jù)來源于仙居縣2008年森林資源二類清查中的抽樣調(diào)查數(shù)據(jù)。按等距抽樣方法進(jìn)行抽樣調(diào)查，抽樣總體為仙居縣行政區(qū)范圍。樣地間距2.0 km×3.0 km，共有樣地307個(gè)，均為28.28 m×28.28 m的正方形，面積0.08 hm2（圖1）。野外樣地測設(shè)由全球定位系統(tǒng)（GPS）定位，水平位置誤差±（5～8）m。樣地調(diào)查內(nèi)容包括樹種、平均年齡、平均胸徑、平均樹高、地類、地貌、海拔、坡向、坡位、坡度等。樣木起測胸徑為5 cm。調(diào)查樣木分為杉木、馬尾松、硬闊、軟闊等4個(gè)樹種（組）。根據(jù)利用浙江省2009年森林資源連續(xù)清查體系（簡稱CFI體系）樣地?cái)?shù)據(jù)建立的單株立木生物量模型計(jì)算地上部分生物量（表1）［13-15］；求和得到各樣地總的地上部分生物量，并根據(jù)生物量/碳儲(chǔ)量轉(zhuǎn)換系數(shù)和樣地面積將地上部分生物量換算成碳密度（Mg·hm-2）；其中生物量/碳儲(chǔ)量轉(zhuǎn)換系數(shù)選用采用政府間氣候變化專門委員會(huì)（IPCC）提供的默認(rèn)參數(shù)0.50計(jì)［5］。通過每木檢尺數(shù)據(jù)計(jì)算出每個(gè)樣地的地上部分生物量，其碳密度值為0～87.005 Mg·hm-2，平均值為15.854 Mg·hm-2，標(biāo)準(zhǔn)差為16.747 Mg·hm-2，變異系數(shù)達(dá)到1.056，屬于強(qiáng)變異類型。中位數(shù)為10.208 Mg·hm-2，呈現(xiàn)出明顯的負(fù)偏（表2）。對碳密度及相關(guān)地形因子進(jìn)行描述性統(tǒng)計(jì)的結(jié)果如表2，其相關(guān)性分析如表3。這幾種統(tǒng)計(jì)量均存在一定的空間變異性，其中碳密度、坡度屬于強(qiáng)變異類型，其他均屬于中等變異類型。在0.01顯著水平下，森林碳密度與坡位、海拔等的相關(guān)系數(shù)為-0.337～0.294。

表1 研究區(qū)碳密度計(jì)算Table 1 Carbon density calculated in the study area

表2 樣本碳密度統(tǒng)計(jì)特征Table 2 Statistical description of the plot carbon density

表3 地形因子與碳密度的相關(guān)性分析Table 3 Correlation between carbon density and the inventory factors

1.2.2 遙感數(shù)據(jù) 研究使用2007年10月仙居縣全境Landsat TM影像數(shù)據(jù)，空間分辨率為30 m×30 m，與樣地大小基本吻合。圖像的預(yù)處理包括常見的衛(wèi)星數(shù)據(jù)處理步驟，如幾何校正和輻射校正等處理，總誤差小于1個(gè)像元。利用軟件Erdas 9.2提取307個(gè)樣地所對應(yīng)的遙感圖像6個(gè)波段灰度值，即band1，band2，band3，band4，band5和band7。通過波段/波段組合與碳密度相關(guān)性分析選擇建模變量。如表4，其中Band3與碳密度呈負(fù)相關(guān)，在0.01顯著水平下達(dá)到-0.436，為最高；植被指數(shù)NDVI，SAVI與碳密度的相關(guān)性在0.01顯著性水平下達(dá)到0.35以上，呈正相關(guān)。

表4 遙感波段及其組合與碳密度相關(guān)性分析Table 4 Correlation between carbon density and the band combination of the remote sensed data

1.2.3 變量篩選 通過相關(guān)分析發(fā)現(xiàn)，坡位，Landsat TM band1，band2，band3，band7，NDVI，SAVI與森林碳密度均在0.01顯著水平下顯著相關(guān)。考慮到其中band3是葉綠素的主要吸收波段，NDVI在森林碳估測中被普遍使用，SAVI為土壤調(diào)整植被指數(shù)，所以將坡位，band3，NDVI，SAVI作為建模的備選因子。由于參與模擬建模的解釋變量不僅要與因變量顯著相關(guān)，而且解釋變量間要求相互獨(dú)立，所以對以上4個(gè)備選因子和碳密度進(jìn)行相關(guān)分析（表5）。結(jié)果表明：解釋變量之間同樣具有顯著相關(guān)性，甚至更高，不適合聯(lián)合建模，若聯(lián)合建模解釋變量之間存在共線性，最終確定選用band3與碳密度進(jìn)行獨(dú)立建模。為了與地理和海拔加權(quán)回歸模型進(jìn)行對比，分別嘗試band3獨(dú)立建模、band3和海拔聯(lián)合建模、地理和海拔加權(quán)回歸建模。

表5 解釋變量與因變量的相關(guān)性分析Table 5 Correlation between the explanatory variables and the dependent variable

2 研究方法

2.1 地理加權(quán)回歸

地理加權(quán)回歸（GWR）方法是傳統(tǒng)回歸方法的擴(kuò)展，可以進(jìn)行局部參數(shù)估計(jì)。其本質(zhì)思想是為數(shù)據(jù)集中各要素建立獨(dú)立的方程［18］。研究區(qū)內(nèi)任何位置的參數(shù)估計(jì)，由1個(gè)因變量和1個(gè)或多個(gè)自變量結(jié)合已知樣本信息的空間關(guān)系構(gòu)成的空間權(quán)重實(shí)現(xiàn)。根據(jù)地理學(xué)第一定律，相互接近的觀測點(diǎn)數(shù)據(jù)貢獻(xiàn)更多的影響，GWR模型以一個(gè)關(guān)于距離的函數(shù)進(jìn)行加權(quán)，函數(shù)的回歸系數(shù)是觀測點(diǎn)地理位置的位置函數(shù)。全局空間回歸假定回歸參數(shù)在整個(gè)研究區(qū)域內(nèi)是一致的，GWR模型的參數(shù)不再保持不變，而是隨著地理位置的變化而變化。研究區(qū)內(nèi)抽樣得到若干觀測點(diǎn)，可視為將整個(gè)研究區(qū)分為這若干區(qū)域，針對各觀測點(diǎn)分別建立模型，作為該觀測點(diǎn)所在區(qū)域的局部回歸方程來進(jìn)行區(qū)域化變量的估計(jì)。GWR的模型結(jié)構(gòu)可表達(dá)為：

式（1）中：（ui,vi）表示點(diǎn)處的坐標(biāo)，表示i點(diǎn)處的因變量，本研究中為i點(diǎn)處碳密度值；n表示變量的數(shù)目，則x1i～xni表示第n個(gè)變量在i點(diǎn)的值；β0表示截距，β0～βn表示第n個(gè)變量的估計(jì)回歸系數(shù)。ε是誤差項(xiàng)。對于傳統(tǒng)全局回歸模型，通過矩陣方程對參數(shù)進(jìn)行估計(jì)：

式（4）中：wij表示附近點(diǎn)j相對于位置i的空間權(quán)重。在空間權(quán)重函數(shù)的選擇上，最常被作為距離的連續(xù)函數(shù)來計(jì)算各點(diǎn)權(quán)重的是Gaussian函數(shù)：

然而由LeSage編寫的GWR程序中的Gaussian函數(shù)表達(dá)式為：

另外，LeSage編寫的GWR程序中還提供了一種Tricube權(quán)重函數(shù)：

本研究將分別嘗試這3種權(quán)重函數(shù)。3個(gè)函數(shù)中，dij表示附近點(diǎn)j相對于位置i的距離，距離通常被定義為歐氏距離：

式（5～7）3個(gè)函數(shù)式中b是一個(gè)非負(fù)參數(shù)，表示帶寬，是GWR模型中最重要的參數(shù)，模型參數(shù)的估計(jì)和模型預(yù)測的準(zhǔn)確度很大程度上取決于帶寬的選擇，因?yàn)閰?shù)的估計(jì)只與帶寬范圍內(nèi)的樣本數(shù)據(jù)有關(guān)，帶寬范圍決定。帶寬通常有2種表現(xiàn)方式：一種是最佳固定距離，另一種是最佳相鄰點(diǎn)數(shù)。就Gaussian函數(shù)而言，不同的權(quán)重函數(shù)表達(dá)式又呈現(xiàn)出不同的帶寬表示形式，式（5）中帶寬b為某一固定距離，式（6）中帶寬b與距離標(biāo)準(zhǔn)差的乘積代表某固定距離，所以，式（6）中的帶寬b應(yīng)是一個(gè)數(shù)值。下文將詳細(xì)對帶寬選擇的結(jié)果進(jìn)行展示；Tricube函數(shù)帶寬則以鄰近點(diǎn)個(gè)數(shù)表示，式（7）中d（b）表示帶寬為b時(shí)，即包含b個(gè)鄰近點(diǎn)，沿橫坐標(biāo)方向到點(diǎn)i最遠(yuǎn)的距離。帶寬的選擇通過Cleveland［19］和Bowman［20］提出的交叉驗(yàn)證（CV）方法［21-22］實(shí)現(xiàn)：

式（9）中：y＾≠i（b）表示帶寬為b時(shí)yi的模型預(yù)測值，其中≠i表示在對yi進(jìn)行模擬的過程中，把yi視為未知，用帶寬b范圍內(nèi)不包括樣本i的其他已知樣本信息對yi模擬建模。最終找到一個(gè)最優(yōu)的帶寬b值，使得函數(shù)CV（b）降到最低。根據(jù)Fotheringham等［20］提出的方法作為準(zhǔn)則：Akaike信息準(zhǔn)則（AAICc）值最小時(shí)，帶寬b為最佳帶寬。

2.2 地理和海拔加權(quán)回歸

地理和海拔加權(quán)回歸（GAWR）與地理加權(quán)回歸（GWR）僅存在一點(diǎn)區(qū)別：如何定義空間關(guān)系。GWR模型通過二維地理空間坐標(biāo)u，v來定義，而GAWR則通過三維地理空間坐標(biāo)來定義［9］：

式（10）中：z代表點(diǎn)i處的海拔，作為第三維度參與計(jì)算空間距離。但多數(shù)情況下，海拔的差異性遠(yuǎn)遠(yuǎn)沒有橫縱坐標(biāo)的差異性明顯。就本研究來說，2個(gè)采樣點(diǎn)間橫坐標(biāo)最小間距為3.0 km，縱坐標(biāo)最小間距為2.0 km，而整個(gè)研究區(qū)的海拔差異最大只有1.0 km左右，此時(shí)式（10）中海拔只提供極小的貢獻(xiàn)。為了改善這種情況，將海拔乘以一個(gè)擴(kuò)大系數(shù)α以強(qiáng)調(diào)垂直方向上的影響因素：

最優(yōu)擴(kuò)大系數(shù)的選擇與帶寬的選擇類似，同樣采用交叉驗(yàn)證的方法，即變異系數(shù)（CV）值最小時(shí)，α的取值為最優(yōu)。在本研究中擴(kuò)大系數(shù)選用α=5。

2.3 回歸參數(shù)的非平穩(wěn)性測試

地理加權(quán)回歸模型是局部回歸模型，不同地理位置上模型的參數(shù)是變化的，于是需要判斷參數(shù)是否在整個(gè)研究區(qū)內(nèi)有顯著變化，即參數(shù)的非平穩(wěn)性測試。本研究采用Fotheringham等［23］提出的方法，將GWR模型的局部參數(shù)與普通最小二乘法回歸的全局固定參數(shù)作對比，若局部估計(jì)的參數(shù)4分位值波動(dòng)大于全局估計(jì)參數(shù)的一倍標(biāo)準(zhǔn)差范圍，則認(rèn)為模型的參數(shù)是顯著非平穩(wěn)的。

2.4 模型性能指標(biāo)

為了檢驗(yàn)GWR模型估計(jì)地上部分生物量的精度，本研究通過決定系數(shù)（R2），殘差平方和（RRSS），均方根誤差（RRMSE）等標(biāo)準(zhǔn)對不同試驗(yàn)結(jié)果進(jìn)行統(tǒng)計(jì)描述。每一預(yù)測點(diǎn)的估計(jì)值與實(shí)際值的差的平方之和，即殘差平方和，用來度量回歸方程和數(shù)據(jù)組之間的總誤差，此測量值越小，模型擬合度越高。均方根誤差為估計(jì)值與實(shí)測值之差的平方的算術(shù)平均值再開方，也稱為標(biāo)準(zhǔn)偏差。它反映了估計(jì)值偏離實(shí)測值的程度，在0.01或0.05顯著性水平下，標(biāo)準(zhǔn)偏差越小，表示估計(jì)精度越高。式（13）和式（14）中，y＾i，yi分別是因變量的估計(jì)值和實(shí)測值；總體平方和（RTSS）由回歸平方和、殘差平方和構(gòu)成。

3 結(jié)果與分析

3.1 模型擬合

本研究借助GWR 4.0和Matlab軟件平臺(tái)，采用Gaussian和Tricube權(quán)重函數(shù)，對307個(gè)樣地碳密度進(jìn)行評估，以估計(jì)值相對于實(shí)測值的RRSS和RRMAE，來檢驗(yàn)GWR模型的估計(jì)精度。軟件GWR 4.0的估計(jì)結(jié)果為GWR（G）′和GWRa（G）′；由Matlab實(shí)現(xiàn)的GWR程序估計(jì)結(jié)果為GWR（G），GWR（T），GWRa（G），GWRa（T），GAWR（G）和GAWR（T），其中下標(biāo)（G）表示Gaussian權(quán)函數(shù)，下標(biāo)（T）表示Tricube權(quán)函數(shù)，下標(biāo)a表示海拔作為解釋變量之一的雙變量GWR模型。

GWR模型參數(shù)的估計(jì)和模型預(yù)測的準(zhǔn)確度很大程度上取決于帶寬的選擇，因?yàn)閰?shù)的估計(jì)只與帶寬范圍內(nèi)的樣本數(shù)據(jù)有關(guān)，所以需對帶寬進(jìn)行優(yōu)化選擇（圖2）。GWR（G）的帶寬默認(rèn)為最小0.1，最大20.0，通過非線性優(yōu)化進(jìn)行交叉驗(yàn)證計(jì)算得分確定最優(yōu)帶寬（圖2A）。GWR（G）′的帶寬為某一固定距離，根據(jù)Akaike信息準(zhǔn)則（AAICc），交叉驗(yàn)證（CV）準(zhǔn)則確定最佳固定距離（圖2B），隨著帶寬的增大，AAICc值和CV值有大幅度的減小，在帶寬為2 000 m左右時(shí)，AAICc值和CV值的變化趨勢趨于平緩，變化幅度減小。對于模型GWR（T），帶寬為待測位置鄰近點(diǎn)的個(gè)數(shù)q，程序默認(rèn)鄰近點(diǎn)個(gè)數(shù)最少應(yīng)為變量數(shù)加2，即CVmin= nvar+2；qmax缺省設(shè)置為4倍的變量數(shù)，可根據(jù)實(shí)際情況改變。

圖2 GWR模型的最優(yōu)帶寬確定Figure 2 Selection of the optimal bandwidth for the GWR model

各模型的結(jié)果均在表6中有所呈現(xiàn)。結(jié)果表明：研究區(qū)內(nèi)碳密度可以通過地理加權(quán)回歸建模得出；單一波段能夠反映森林地上部分碳儲(chǔ)量，適當(dāng)增加解釋變量的個(gè)數(shù)可以得到更好的估計(jì)結(jié)果，GWRa（G）′和GWRa（T）模型給出最小的殘差平方和、均方根誤差以及最高的R2，但由于存在多重共線性問題，海拔系數(shù)的t檢驗(yàn)結(jié)果不顯著；其次，GWR（G）′模型、GWR（T）模型給出較好的估計(jì)結(jié)果，決定系數(shù)可達(dá)到0.65以上，且估計(jì)結(jié)果在0.01顯著水平下顯著。

對于所有GWR模型和GAWR模型而言，選用的權(quán)重函數(shù)不同，帶寬的表述形式也不同，建模的結(jié)果也不盡相同。GWR（G）′模型的估計(jì)結(jié)果與GWR（T）的估計(jì)結(jié)果較為相近。基于Tricube權(quán)重函數(shù)的模型總能模擬出最大的碳密度范圍，GWR（T）估計(jì)結(jié)果的變異系數(shù)達(dá)到0.779，相較實(shí)測碳密度保留了原數(shù)據(jù)70%以上的變異特征。值得注意的是：GWR（G）模型估計(jì)結(jié)果與另外2個(gè)估計(jì)結(jié)果相差甚遠(yuǎn)。GWRa（G）模型比相對應(yīng)的單變量模型估計(jì)精度不增反減，加入海拔作為解釋變量之后，模型的估計(jì)精度降低了2.2%。另外2種模型精度略有提升，GWRa（T）模型比相對應(yīng)的單變量模型提升了14.0%，GWRa（G）′模型提升9.0%。

將海拔信息加入到待測點(diǎn)的空間權(quán)重計(jì)算中，可以避免多元回歸建模中的多重共線性；并且針對海拔因素貢獻(xiàn)較小的問題，可通過加入一個(gè)擴(kuò)大系數(shù)進(jìn)行改善，即GWRS（G)和GAWRS（T）模型。然而針對本研究區(qū)，擴(kuò)大系數(shù)的加入并沒有起到明顯的作用，對于GAWR（T）模型來說也只有微弱的改善，擴(kuò)大系數(shù)的貢獻(xiàn)僅有0.046。再對比GAWR與GWR的建模，結(jié)果表明：空間權(quán)重計(jì)算中海拔信息的加入，一定程度上降低了GWR回歸模型的預(yù)測精度。GAWR5（T）模型比GWR（G）模型降低了0.6%；GAWR5（T）模型比GWR（T）模型降低10.4%，比GWRa（T）模型降低24.4%。對于本研究區(qū)來說海拔影響不明顯，GAWR模型不適用。

表6 GWR模型、GAWR模型森林碳密度估計(jì)統(tǒng)計(jì)Table 6 Descriptive statistics of the predictive forest carbon density models by GWR and GAWR

3.2 回歸系數(shù)的非平穩(wěn)性

GWR模型參數(shù)波動(dòng)范圍的劇烈程度是評判模型參數(shù)在整個(gè)研究區(qū)內(nèi)是否有顯著變化的標(biāo)準(zhǔn)。若模型參數(shù)非平穩(wěn)性顯著，則表明該模型能反映更多的空間變異信息。如圖3，以 GWR（G）和GWR（T）模型的參數(shù)范圍為例，GWR（G）模型截距和自變量系數(shù)曲線是平穩(wěn)的，2個(gè)參數(shù)的變化范圍分別為16.430～35.320和-0.195～0.026；說明此GWR模型的系數(shù)局部差異不明顯，幾乎與全局回歸模型無異。這也就解釋了GWR（G）模型在估計(jì)研究區(qū)森林碳密度時(shí)精度較低的原因。

圖3 GWR（G）和GWR（T）回歸參數(shù)曲線Figure 3 Local regression parameter curve for GWR（G）and GER（T）

相對來說GWR（T）模型參數(shù)變化幅度明顯優(yōu)于GWR（G）模型（圖3）。由表7可知：該模型截距的取值范圍為-50.595～147.780，而不是全局回歸模型中的恒定值26.718；band3的系數(shù)取值范圍是-0.912～1.224，全局回歸模型中為恒定值-0.124。對比表7和表8得出結(jié)論：GWR（T）截距的第1四分位與第3四分位的間距27.177小于全局回歸截距的1倍標(biāo)準(zhǔn)差間距54.169；band3系數(shù)的第1四分位與第3四分位的間距0.225小于全局回歸Band3系數(shù)的1倍標(biāo)準(zhǔn)差間距0.514。所以針對本研究區(qū)，GWR（T）模型參數(shù)存在一定的空間非平穩(wěn)特征，但不顯著，相對于全局回歸能反映更多的空間變異信息。

表7 GWR（T）回歸參數(shù)描述Table 7 Descriptive statistics of the local regression parameters from GWR（T）

表8 全局回歸參數(shù)描述Table 8 Descriptive statistics of the global regression parameters

3.3 仙居縣全縣總碳儲(chǔ)量及其分布估計(jì)

按照樣地?cái)?shù)據(jù)統(tǒng)計(jì)，仙居全縣森林地上部分碳儲(chǔ)量為3.192×106Mg，平均碳密度為15.854 Mg·hm-2。以GWR（T）模型為例估計(jì)仙居縣全縣森林地上部分碳儲(chǔ)量及其分布（圖4），由于地上部分森林碳密度值≥0，而由GWR（T）模型估計(jì)結(jié)果中出現(xiàn)極少數(shù)負(fù)值現(xiàn)象，估計(jì)結(jié)果為負(fù)的樣地主要集中于仙居縣城鎮(zhèn)聚居區(qū)（圖1），實(shí)測碳密度接近0，于是將估計(jì)結(jié)果中負(fù)值替換為0。還原0值后研究區(qū)碳密度為0～89.964 Mg·hm-2，分布趨勢與實(shí)際情況相符：全縣碳密度整體偏低，低值區(qū)域以城鎮(zhèn)聚居區(qū)為中心向周邊擴(kuò)散，高值區(qū)域出現(xiàn)在研究區(qū)東北、西北及西南部；碳密度平均值為 15.555 Mg·hm-2，標(biāo)準(zhǔn)差為11.332 Mg·hm-2。研究區(qū)內(nèi)森林地上部分總碳儲(chǔ)量估計(jì)值為3.132×106Mg。總碳儲(chǔ)量估計(jì)結(jié)果低于樣地?cái)?shù)據(jù)統(tǒng)計(jì)結(jié)果1.880%，吻合度較高。

圖4 碳密度分布估計(jì)圖Figure 4 Carbon density estimation

以307個(gè)樣地碳密度估計(jì)值與實(shí)測值的殘差平方和（RRSS）和均方根誤差（RRMSE），來檢驗(yàn)地理加權(quán)回歸模型的估計(jì)精度（圖5，表9）。圖5為樣地碳儲(chǔ)量實(shí)測值與相對應(yīng)估計(jì)結(jié)果的對比，可以看出擬合程度好，估計(jì)結(jié)果與實(shí)際情況相一致，碳密度分布差異大，且大部分地區(qū)碳密度較低；碳密度估計(jì)值相對于實(shí)測值的RRSS為29 497.265 4 Mg·hm-2，RRMSE為9.802 Mg· hm-2。由圖6估計(jì)偏差曲線可以看出：GWR模型單個(gè)點(diǎn)的殘差波動(dòng)范圍較大，為-32.812～32.280 Mg· hm-2。對比圖5和圖6可以發(fā)現(xiàn)：當(dāng)實(shí)測值較大時(shí)，更容易得到較大的估計(jì)誤差，是由于在局部建模回歸過程中，帶寬范圍內(nèi)的樣本數(shù)據(jù)對預(yù)測點(diǎn)的局部平滑作用，研究區(qū)內(nèi)碳儲(chǔ)量普遍偏低，個(gè)別具有較大碳密度的樣地被低估。同理，具有較小碳密度的樣地也在一定程度上被高估，但程度較低。

3.4 與其他方法的對比

為了驗(yàn)證地理加權(quán)回歸方法在估計(jì)區(qū)域森林碳儲(chǔ)量時(shí)的優(yōu)勢，我們用相同樣本數(shù)據(jù)分別進(jìn)行傳統(tǒng)全局回歸估計(jì)、協(xié)同克里格插值，將樣地的預(yù)測結(jié)果對比分析見表9。傳統(tǒng)的全局回歸在反映遙感數(shù)據(jù)與生物量的關(guān)系中會(huì)丟失部分局部信息，通過全局統(tǒng)計(jì)導(dǎo)致較差的生物量估計(jì)。表9中全局回歸、地理加權(quán)回歸估計(jì)結(jié)果均值與實(shí)測碳密度值相一致，但變異系數(shù)差異明顯，表明全局回歸反映了研究區(qū)內(nèi)區(qū)域化變量的平均狀態(tài)；相對而言地理加權(quán)回歸在保證平均水平的同時(shí)，更能反映空間數(shù)據(jù)變化規(guī)律的真實(shí)特征，保留了70%以上的空間異質(zhì)性。在空間異質(zhì)性分析方面，空間插值方法應(yīng)用較廣泛，主要依賴量化區(qū)域化變量的自相關(guān)性，通過相似的變異函數(shù)特征利用變程范圍內(nèi)的數(shù)據(jù)對未知位置進(jìn)行預(yù)測。協(xié)同克里格插值結(jié)果雖保留了近70%的空間異質(zhì)性特征，但從碳密度取值范圍、標(biāo)準(zhǔn)差、誤差分析等來看，預(yù)測精度較低。地理加權(quán)回歸的預(yù)測結(jié)果誤差小，取值范圍跨度大，更適合與對森林碳儲(chǔ)量空間變異性的分析和預(yù)測。地理加權(quán)回歸模型具有較高的空間模擬預(yù)測精度。

圖5 樣地碳密度估計(jì)值與實(shí)測值比較Figure 5 Carbon densities from GWR estimation and forest inventory

圖6 樣地碳密度估計(jì)偏差曲線Figure 6 Deviation curves of the forest carbon density estimate

表9 不同方法估計(jì)樣地碳密度的對比Table 9 Comparison of carbon density estimation in the study area by different methods

4 結(jié)論與討論

基于GWR（T）的2008年仙居縣森林碳儲(chǔ)量估計(jì)值為3.132×106Mg，平均碳密度為15.555 Mg·hm-2，比實(shí)測碳密度低1.880%，總量及其分布與實(shí)際一致；研究區(qū)內(nèi)碳密度的取值范圍為0～89.964 Mg·hm-2，保留了實(shí)測碳密度70%以上的空間異質(zhì)性特征；決定系數(shù)0.65以上，且估計(jì)結(jié)果在0.01顯著水平下顯著?；诘乩砑訖?quán)回歸的森林地上部分碳儲(chǔ)量估計(jì)方法是有效的，其估計(jì)結(jié)果精度高于傳統(tǒng)回歸和協(xié)同克里格插值方法。

本研究對地理與海拔加權(quán)回歸模型（GAWR）在地勢平緩地區(qū)的應(yīng)用進(jìn)行檢測；作為與GAWR模型的對比，嘗試了波段、海拔雙變量的GWR模型建模。結(jié)果表明：GAWR模型的應(yīng)用具有地域的特殊性，并不是針對所有地形地貌的區(qū)域GAWR模型都存在優(yōu)勢。這與海拔對植被的生長、分布的影響密切相關(guān)。在海拔對植被的影響不顯著的平緩地區(qū)，GAWR并不適用；將海拔作為解釋變量加入建模，能夠提高碳儲(chǔ)量估計(jì)精度，但存在多重共線性問題，使得雙變量模型t檢驗(yàn)不顯著。

地理加權(quán)回歸模型在估計(jì)區(qū)域森林碳儲(chǔ)量時(shí)仍有一些值得注意和思考的問題：①在GWR估計(jì)中，帶寬是影響估計(jì)結(jié)果精度的主要因素之一。盡管有Akaike信息準(zhǔn)則和交叉驗(yàn)證準(zhǔn)則在理論上指導(dǎo)尋找最優(yōu)帶寬，但實(shí)際操作中仍需結(jié)合具體情況綜合考慮?？臻g權(quán)重計(jì)算方法的差異也會(huì)導(dǎo)致帶寬的差異，GWR模型對帶寬的變化非常敏感，所以針對不同研究區(qū)、研究數(shù)據(jù)選擇合適的權(quán)重函數(shù)以確定最佳帶寬是保證GWR模型精度的根本。②森林碳儲(chǔ)量的大小與空間分布特征直接影響森林的碳匯潛力。目前，面臨的一個(gè)挑戰(zhàn)是，怎樣從點(diǎn)測量尺度上推到更大的國家、區(qū)域和全球尺度上［24］。已有研究表明，地理加權(quán)回歸分析對粒度變化呈現(xiàn)出一定的穩(wěn)健性［25］，因此，將其用于有明顯地理差異的大尺度區(qū)域森林碳儲(chǔ)量及其分布估計(jì)中，能取得更顯著的效果。③目前。基于GWR模型的擴(kuò)展模型還有：半?yún)?shù)地理加權(quán)回歸（semiparametric GWR）模型，一部分自變量的參數(shù)為局部回歸，另一部分自變量的參數(shù)為全局回歸；地理加權(quán)泊松回歸（GWPR）模型、地理加權(quán)Logistic回歸（GWLR）模型等［26］。這些模型是否也適用于林業(yè)遙感領(lǐng)域，用以估計(jì)森林碳儲(chǔ)量，是否能得到更好的估計(jì)精度，都值得進(jìn)一步探索。

［1］ PHILIPS O L,MALHI Y,HIGUCHI N,et al.Changes in the carbon balance of tropical forests：evidence from longterm plots［J］.Science,1998,282（5388）：439-442.

［2］ 李?？?，雷淵才，曾偉生.基于森林清查資料的中國森林植被碳儲(chǔ)量［J］.林業(yè)科學(xué)，2011，47（7）：7-12.

LI Haikui,LEI Yuancai,ZENG Weisheng.Forest carbon storage in china estimated using forestry inventory data［J］. Sci Silv Sin,2011,47（7）：7-12.

［3］ 王宏鈺，沈歡.國內(nèi)外森林碳儲(chǔ)量估測研究現(xiàn)狀［J］.科技信息，2012（18）：77.

WANG Hongyu,SHEN Huan.Estimation of forest carbon stocks both at home and abroad［J］.Sci Technol Inf,2012（18）：77.

［4］ BROOK E J,SOWERS T,ORCHARDO J.Rapid variations in atmospheric methane concentration during the past 110 000 years［J］.Science,1996,273（5278）：1087-1091

［5］ WATSON R T,NOBLE I R,BOLIN B,et al.Land Use,Land-Use Change and Forestry:Special Report of the IPCC［R］.Cambridge：Cambridge University Press,2000：30-31.

［6］ 曹吉鑫，田赟，王小平，等.森林碳匯的估算方法及其發(fā)展趨勢［J］.生態(tài)環(huán)境學(xué)報(bào)，2009，18（5）：2001-2005.

CAO Jixin,TIAN Yun,WANG Xiaoping,et al.Estimation methods of forest sequestration and their prospects［J］. Ecol Environ Sci,2009,18（5）：2001-2005.

［7］ 李明陽，劉敏，劉米蘭.基于GIS的森林調(diào)查因子地統(tǒng)計(jì)學(xué)分析 ［J］.南京林業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2010，34（6）：66-70.

LI Mingyang,LIU Min,LIU Milan.GIS-based geostatistical analysis of forest inventory factors［J］.J Nanjing For Univ Nat Sci Ed,2010,34（6）：66-70.

［8］ 李明詩，譚瑩，潘潔，等.結(jié)合光譜、紋理及地形特征的森林生物量建模研究［J］.遙感信息，2006（6）：6-9.

LI Mingshi,TAN Ying,PAN Jie,et al.Modeling forest aboveground biomass by combining the spectrum,textures with topographic features［J］.Rem Sens Inf,2006（6）：6-9.

［9］ PROPASTIN P.Modifying geographically weighted regression for estimating aboveground biomass in tropical rainforests by multispectral remote sensing data［J］.Int J Appl Earth Observ Geoinform,2012,18：82-90.

［10］ 邵一希，李滿春，陳振杰，等.地理加權(quán)回歸在區(qū)域土地利用格局模擬中的應(yīng)用：以常州市孟河鎮(zhèn)為例［J］.地理科學(xué)，2010，30（1）：92-97.

SHAO Yixi,LI Manchun,CHEN Zhenjie,et al.Simulation on regional spatial land use patterns using geographically weighted regression：a case study of Menghe Town,Changzhou［J］.Sci Geogr Sin,2010,30（1）：92-97.

［11］ 郭龍，張海濤，陳家贏，等.基于協(xié)同克里格插值和地理加權(quán)回歸模型的土壤屬性空間預(yù)測比較［J］.土壤學(xué)報(bào)，2012，49（5）：1037-1042.

GUO Long,ZHANG Haitao,CHEN Jiaying,et al.Comparison between co-kriging model and geographocally weighted regression model in spatial arediction of soil attributes［J］.Acta Pedol Sin,2012,49（5）：1037-1042.

［12］ PROPASTIN P.Modifying geographically weighted regression for estimating aboveground biomass in tropical rainforests by multispectral remote sensing data［J］.Int J Appl Earth Observ Geoinf,2012,18（1）：82-90.

［13］ 朱湯軍，沈楚楚，季碧勇，等.基于LULUCF溫室氣體清單編制的浙江省杉木林生物量換算因子［J］.生態(tài)學(xué)報(bào)，2013，33（13）：3925-3932.

ZHU Tangjun,SHEN Chuchu,JI Biyong,et al.Research on biomass expansion factor of Chinese fir forest in Zhejiang Province based on LULUCF greenhouse gas inventory［J］.Acta Ecol Sin,2013,33（13）：3925-3932.

［14］ 沈楚楚，朱湯軍，季碧勇，等.浙江省馬尾松林生物量換算因子研究［J］.浙江林業(yè)科技，2013,33（3）：39-42.

SHEN Chuchu,ZHU Tangjun,JI Biyong,et al.Research on biomass expansion factor of Pinus massoniana forest in Zhejiang［J］.J Zhejiang For Sci Technol,2013,33（3）：39-42.

［15］ 沈楚楚.浙江省主要樹種（組）生物量轉(zhuǎn)換因子研究［D］.臨安：浙江農(nóng)林大學(xué)，2013.

SHEN Chuchu.The Research on Biomass Expansion Factors of the Dominant Tree Species in Zhejiang Province［D］.Lin’an：Zhejiang A&F University,2013.

［16］ WEISTER R L,ASRAR G,MILLER G P,et al.Assessing grassland biophysical characteristics from spectral measurements［J］.Rem Sens Environ,1986,20（2）：141-152.

［17］ JOHN G L,DING Y,ROSS S L,et al.A change detection experiment using vegetation indices［J］.Photogr Eng Rem Sens,1998,64（2）：143-150.

［18］ CHARLTON M,FOTHERINGHAM S.Geographically Weighted Regression:White Paper［M］.Maynooth：National Centre for Geocomputation,National University of Ireland Maynooth,2009：5-11.

［19］ CLEVELAND W S.Robust locally weighted regression and smoothing scatterplots［J］.J Amer Stat Assoc,1979,74（368）：829-836.

［20］ BOWMAN A W.An alternative method of cross-validation for the smoothing of density estimates［J］.Biometrika, 1984,71（2）：353-360.

［21］ SALVATI N,TZAVIDIS N,PRATESI M,et al.Small area estimation via M-quantile geographically weighted regression［J/OL］.TEST,2012,21（1）：1-28.

［22］ 艾福利，龐西磊，湯慶園，等.基于地理加權(quán)回歸模型的浙江省臺(tái)風(fēng)經(jīng)濟(jì)損失影響因子分布規(guī)律研究［J］.北京師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2013，49（1）：61-67.

AI Fuli,PANG Xilei,TANG Qingyuan,et al.Impact factor distribution for economic losses due to Typhoon in Zhejiang province by geographical weighted regression［J］.J Beijing Norm Univ Nat Sci,2013,49（1）：61-67.

［23］ FOTHERINGHAM A S,BRUNSDON C,CHARLTON M.Geographically Weighted Regression：The Analysis of Spatially Varying Relationship［M］.West Sussex：Wiley,2002.

［24］ WANG Guangxing,OYANA T,ZHANG Maozhen,et al.Mapping and spatial uncertainty analysis of forest vegetation carbon by combining national forest inventory data and satellite images［J］.For Ecol Manage,2009,258（7）：1275-1283.

［25］ 覃文忠.地理加權(quán)回歸基本理論與應(yīng)用研究［D］.上海：同濟(jì)大學(xué)，2007.

TAN Wenzhong.The Basic Theoretics and Application Research on Geographically Weighted Regression［D］. Shanghai：Tongji University,2007.

［26］ NAKAYA T.GWR4：Windows Application for Geographically Weighted Regression Modeling［R］.Liverpool：GWR 4 Development Team,2012.

Geographically weighted regression based on estimation of regional forest carbon storage

GUO Hanru1,2,ZHANG Maozhen1,2,XU Lihua1,2,YUAN Zhenhua1,2,CHEN Tiange3
（1.Zhejiang Provincial Key Laboratory of Carbon Cycling in Forest Ecosystems and Carbon Sequestration, Zhejiang A&F University,Lin’an 311300,Zhejiang,China;2.School of Environmental&Resource Sciences, Zhejiang A&F University,Lin’an 311300,Zhejiang,China;3.Hualong District Science&Technology Bureau of Puyang City,Puyang 457001,Henan,China）

Global climate issues have confirmed the irreplaceable role forest carbon stocks play in the global carbon cycle.To research whether the geographically weighted regression（GWR）model method which considers the role of survey factors’spatial heterogeneity and establish the local regression model,can improve the estimation accuracy of forest carbon stocks,instead of the more commonly used methods of global regression model such as ordinary least squares analysis（OLS）,we used forest management inventory data in Xianju County,Zhejiang Province,combined with Landsat TM image data developing local models using GWR to estimate forest carbon stock and its density.Available of geographically and altitudinal weighted regression（GAWR）model was then tested in smooth terrain.Analysis is included comparison to traditional regression and co-kriging interpolation.Results showed that the total forest aboveground carbon stocks estimated by the GWR（T）model for Xianju County were 3.132×106Mg,and carbon density ranged from 0 to 89.964 Mg·hm-2with a mean value of 15.555 Mg·hm-2.Meanwhile,the total forest aboveground carbon stocks calculated from diameter measurements were 3.192×106Mg with a mean value of 15.854 Mg·hm-2.The overall result from GWR（T）model was lower than diameter measured by 1.880%,R2=0.654（P＜0.01）,and carbon density distribution was consistent with the actual situation.The estimated results also had a higher accuracy with the RRMSE=9.802（P＜0.01）than traditional regression method with the RRMSE=15.033（P＜0.01）and co-kriging interpolation method with the RRMSE=16.427（P＜0.01）.GWR method can effectively estimate the regional forest aboveground carbon stocks reasonably and accurately,however,the GAWR model is not applicable for the areas with smooth terrain.Adding altitude as an explanatory variable in the modeling could improve estimation accuracy but would in turn create a multi-collinearity problem.［Ch,6 fig.9 tab.26 ref.］

forest ecology;forest carbon storage;spatial heterogeneity;geographically weighted regression（GWR）;geographically and altitudinal weighted regression（GAWR）

S718.55,S757.2

A

2095-0756（2015）04-0497-12

10.11833/j.issn.2095-0756.2015.04.002

2014-09-19；

2014-11-01

國家自然科學(xué)基金資助項(xiàng)目（30972360，41201563）；浙江省林業(yè)碳匯與計(jì)量創(chuàng)新團(tuán)隊(duì)資助項(xiàng)目（2012R10030-01）

郭含茹，從事森林碳匯計(jì)量和監(jiān)測技術(shù)研究。E-mail：15869132363@163.com。通信作者：張茂震，教授，博士，從事森林碳匯計(jì)量和監(jiān)測技術(shù)研究。E-mail：zhangmaozhen@163.com