王書賢,張加龍,鮑 瑞,韓東陽,廖 易

(1.西南林業(yè)大學(xué) 林學(xué)院,云南 昆明 650224;2.云南省遙感中心,云南 昆明 650034;3.國家林業(yè)和草原局 西南調(diào)查規(guī)劃院,云南 昆明 650031)

森林地上生物量(abovegroundbiomass,AGB)是估算森林生態(tài)系統(tǒng)碳儲量、碳通量的重要參數(shù)[1-2],森林AGB空間分布變化能夠更好地量化碳儲量的空間分異特征,對于研究森林生態(tài)系統(tǒng)碳循環(huán)具有重要意義?，F(xiàn)階段主要運(yùn)用ArcGIS軟件中的空間插值法來對森林AGB的空間異質(zhì)性進(jìn)行分析研究,而用遙感影像數(shù)據(jù)反演森林AGB,然后對其進(jìn)行空間異質(zhì)性分析的研究較少[3-6]?；谶b感與地統(tǒng)計(jì)方法,能夠有效揭示屬性變量在空間上的分布特征及其變異性[7-13],當(dāng)前運(yùn)用空間地統(tǒng)計(jì)學(xué)方法的林業(yè)研究主要集中于森林土壤空間變異、物種空間分布、森林干擾、林分因子和種子資源這五個(gè)方面[7,14],而森林AGB空間的統(tǒng)計(jì)分析研究還較少。在森林AGB的估算中,AGB的空間關(guān)系未引起研究者的重視。當(dāng)前估算森林AGB的方法還是以統(tǒng)計(jì)分析為主,這種方法不但不能準(zhǔn)確、及時(shí)地得到森林AGB的空間分布信息,而且也不能有效反映森林AGB的空間自相關(guān)性和異質(zhì)性[13,15]。

香格里拉地處青藏高原地區(qū)東南部,其氣候變化是全球氣候變化的敏感響應(yīng)區(qū)[16]。高山松(Pinusdensata)是該地區(qū)的優(yōu)勢樹種之一,探究其AGB的時(shí)空變化特征對研究高原地區(qū)生態(tài)系統(tǒng)的碳循環(huán),提高對該類型地區(qū)森林生態(tài)價(jià)值的認(rèn)識具有重要意義[15]。本研究以香格里拉市高山松區(qū)域7個(gè)時(shí)期(1987-2017年)的Landsat遙感影像為基礎(chǔ),對不同階段高山松AGB采用地理信息系統(tǒng)(GIS)和地統(tǒng)計(jì)學(xué)2種方法進(jìn)行時(shí)空變異分析,探究該地區(qū)高山松AGB的空間變化特征,為日后監(jiān)測AGB空間變化和模擬預(yù)測提供理論依據(jù)。

1 研究區(qū)概況

研究區(qū)位于云南省迪慶藏族自治州香格里拉市(99°20′-100°19′E,26°52′-28°52′N),地處云南省西北部,迪慶州東北部。境內(nèi)海拔最高5 545 m,最低1 503 m,相對高差4 042 m。香格里拉市位于高海拔低緯度地帶,屬山地寒溫帶季風(fēng)氣候。境內(nèi)地貌復(fù)雜多樣,地勢北高南低,立體氣候明顯,植被南北分布差異明顯。研究區(qū)植被類型主要為寒溫性針葉林,其樹種有高山松、麗江云杉(Picealikiangensis)等,其中高山松森林覆蓋面積占優(yōu)勢樹種面積的22.7%[15,17]。

2 材料與方法

2.1 數(shù)據(jù)來源

2.1.1 樣地?cái)?shù)據(jù) 樣地?cái)?shù)據(jù)來自國家森林資源連續(xù)清查項(xiàng)目,時(shí)間為1987-2017年,每隔5 a進(jìn)行一次樣地實(shí)測,共7期。每塊調(diào)查樣地面積為0.08 hm2(28.28 m×28.28 m),共有136塊高山松樣地,各年份樣地?cái)?shù)量各異,部分樣地為復(fù)測樣地。

2.1.2 遙感數(shù)據(jù) 采用美國地質(zhì)勘探局(United states geological survey,USGS)和美國國家航空航天局(national aeronautics and space administration,NASA)已經(jīng)過校準(zhǔn)的Landsat地表反射率數(shù)據(jù)(SR)作為時(shí)間序列遙感數(shù)據(jù),運(yùn)用可提高數(shù)據(jù)精度的最優(yōu)斜率算法(the best slope temporal segmentation filtering,BSTS)對數(shù)據(jù)進(jìn)行了時(shí)間序列濾波重構(gòu)[17]。

2.1.3 其他數(shù)據(jù) DEM(數(shù)字高程模型)數(shù)據(jù)來自地理空間數(shù)據(jù)云(http://www.gscloud.cn/),分辨率為30 m×30 m。

2.2 基于隨機(jī)森林回歸的生物量反演結(jié)果

2.2.1 遙感特征變量的選擇 遙感特征變量的選取參考了香格里拉高山松AGB估測研究中常用的遙感因子,紋理因子是基于ENVI 5.3軟件對Landsat的6個(gè)波段分別計(jì)算了7種紋理因子,總共計(jì)算了3×3、5×5、7×7、9×9、11×11、13×13、15×15、17×17、19×19共 9 種不同的窗口大小。由于紋理特征種類過多,單獨(dú)對 378種紋理特征進(jìn)行了篩選。將全部的 378個(gè)特征因子數(shù)據(jù)輸入隨機(jī)森林模型中,通過迭代篩選,選出15個(gè)遙感與紋理特征(表1)。

表1 15個(gè)遙感因子特征變量Table 1 Characteristic variables of 15 remote sensing factors

2.2.2 多期高山松AGB估測數(shù)據(jù) 使用篩選出的15個(gè)遙感與紋理特征變量,基于隨機(jī)森林估測模型對1987、1992、1997、2002、2007、2012、2017年香格里拉高山松地上生物量進(jìn)行估測反演,得到各年份的高山松地上生物量反演數(shù)據(jù)結(jié)果。

2.3 地統(tǒng)計(jì)分析方法

在ArcGIS 10.7中對7期的高山松AGB隨機(jī)取點(diǎn),樣本量取值為200個(gè)隨機(jī)點(diǎn),進(jìn)行標(biāo)準(zhǔn)化異常值處理后,各年份樣本量見表2。應(yīng)用GS+9.0地理空間分析軟件對高山松AGB的隨機(jī)點(diǎn)計(jì)算Moran'sI系數(shù)和半變異函數(shù),然后對其進(jìn)行空間自相關(guān)性和空間異質(zhì)性分析。

表2 隨機(jī)點(diǎn)在各年份中的分配Table 2 Distribution of randomplots in each year

2.3.1 空間自相關(guān)性分析方法 通過Moran'sI系數(shù)計(jì)算高山松AGB和遙感因子的空間自相關(guān)度,并用Moran'sI系數(shù)圖進(jìn)行分析,計(jì)算公式[5,9,11,18-20]如下。

(1)

Moran'sI系數(shù)的取值范圍在( 0,1],說明高山松AGB和遙感因子的空間分布呈正相關(guān)是聚集化狀態(tài),有整體性分布; Moran'sI系數(shù)的取值范圍在[-1,0),說明高山松AGB和遙感因子的空間分布呈負(fù)相關(guān),是破碎化狀態(tài),沒有整體性分布。

2.3.2 空間異質(zhì)性分析方法 運(yùn)用半變異方差函數(shù)對高山松AGB和遙感因子進(jìn)行空間異質(zhì)性分析,計(jì)算公式[5,9,11,18-20]如下。

(2)

式中:N(h)表示某一空間方位上相距距離為h的固定樣地對數(shù);Z(xi)和Z(xi+h)分別表示固定樣地xi和與xi相隔距離h處的生物量值。

采用4個(gè)參數(shù)進(jìn)行空間異質(zhì)性分析,分別為塊金值C0、基臺值(C0+C)、變程A0、結(jié)構(gòu)比C0/(C0+C)。塊金值C0為區(qū)域變量的方差;基臺值(C0+C)反映區(qū)域化變量在研究范圍內(nèi)變異的強(qiáng)度[21];變程A0表示變量空間自相關(guān)的有效范圍,距離越近,空間特征的相似程度越高,超出該值范圍則變量沒有自相關(guān)關(guān)系;結(jié)構(gòu)比C0/(C0+C)判斷隨機(jī)因子對空間總變異的影響程度[9],也可表示變量的空間相關(guān)性程度,若比例<25%,說明變量具有強(qiáng)烈的空間自相關(guān);若比例在25%～75%,變量具有中等空間自相關(guān);比例>75%時(shí),變量空間自相關(guān)很弱[22]。

3 結(jié)果與分析

3.1 高山松AGB的時(shí)間變化

根據(jù)AGB反演結(jié)果(表3),1987-1997年的高山松分布面積減小,由17.156萬hm2減少到17.059萬hm2,期間的AGB也在快速減少,由951.911萬t減少到810.696萬t;1997年后香格里拉市高山松的分布面積開始增加,截至2017年,高山松分布面積達(dá)18.484萬hm2,但該區(qū)域的高山松AGB并未隨之相應(yīng)增加,而是呈減少趨勢,2017年的高山松AGB總量為741.084萬 t。

表3 1987-2017年高山松AGB反演結(jié)果Table 3 Inversion results of P.densata AGB during 1987-2017

3.2 高山松AGB的空間自相關(guān)性

從圖1可以看出,1987年高山松AGB在55 000 m范圍內(nèi)為聚集狀態(tài),呈連片化、整體性空間分布,在南北方向的聚集狀態(tài)為最好,空間格局尺度達(dá)85 000 m左右;1992年高山松AGB在60 000 m內(nèi)為聚集狀態(tài),在南北方向的聚集狀態(tài)為最好,空間格局尺度為75 000 m,但總體的AGB空間分布較為破碎;1997年高山松AGB空間自相關(guān)性較差,整體呈破碎化空間分布,不具有整體性;2002年和2007年的高山松AGB在20 000～50 000 m都呈現(xiàn)破碎化空間分布,在南北方向和東西方向也在該空間距離內(nèi)呈現(xiàn)破碎化的空間分布;2012年在20 000 m內(nèi)和50 000～70 000 m具有聚集狀態(tài)的高山松AGB,其南北方向、東北-西南方向和西北-東南方向也在該空間距離內(nèi)具有同樣空間分布狀態(tài)的高山松AGB;2017年的高山松AGB在不同空間距離上具有不同的空間分布,全局和各方向的空間分布趨勢大體一致。

圖1 1987-2017年高山松AGB的Moran's I系數(shù)Fig.1 Moran's I coefficient of P.densata AGB from 1987 to 2017

3.3 高山松AGB的空間異質(zhì)性

從表4可以看出:經(jīng)過GS+9.0的篩選,1987-2012年的AGB的最優(yōu)半變異函數(shù)模型均為指數(shù)模型,其結(jié)構(gòu)比C0/(C0+C)都接近0。除1992年和1997年的變程A0>10 km外,其余年份的變程A0都在10 000 m內(nèi)。經(jīng)過變程A0篩選出與各年份AGB具有空間分異特征相似的遙感因子,分別是NBR、B3/Albedo、FR11B2VA和FR17B7DS,其中,NBR和B3/Albedo為出現(xiàn)頻率最高的遙感因子;由于2007年AGB與遙感因子的變程A0相差較大,因此該年沒有與其具有相似空間分異特征的遙感因子;2017年高山松AGB的半變異函數(shù)模型的最優(yōu)模型為球形模型,其結(jié)構(gòu)比C0/(C0+C)是近30 a來最接近0的,說明2017年香格里拉高山松AGB的空間自相關(guān)程度非常高。根據(jù)變程A0的相近程度,具有與2017年AGB相似空間分異特征的遙感因子為B3/Albedo。

表4 1987-2017年半變異參數(shù)Table 4 Semi-variation parameter table in 1987-2017

4 結(jié)論與討論

1987-1992年,南北方向的高山松AGB在60 000 m內(nèi)的空間自相關(guān)性最好,呈聚集狀態(tài),高山松為連片化分布;1997-2007年的高山松AGB全局空間自相關(guān)性差,呈破碎化空間狀態(tài);而2017年的高山松AGB在不同空間距離上呈現(xiàn)出不同的空間分布狀態(tài)。

森林AGB的空間異質(zhì)性產(chǎn)生的主要原因可能是立地條件的空間異質(zhì)性[9-10]。1987年和2002-2017年的結(jié)構(gòu)比C0/(C0+C)都接近0,變程A0都在10 000 m內(nèi),說明該時(shí)期的高山松AGB和遙感因子變量的異質(zhì)性變化是由空間自相關(guān)引起的,即其異質(zhì)性變化是由結(jié)構(gòu)性因素(立地條件,海拔、坡度等)及水熱條件等自然因素引起的。該時(shí)期的AGB空間變異受人為干擾因素(砍伐、污染等)的影響較小,也說明1987年和2002-2017年香格里拉高山松地區(qū)的森林保護(hù)良好。

1992年和1997年的C0/(C0+C)雖然接近于0,但變程A0>10 000 m,說明這2個(gè)時(shí)期的高山松AGB和遙感因子變量的異質(zhì)性變化可能是由人為干擾等外界因素引起的。結(jié)合空間自相關(guān)可判斷該時(shí)期的AGB不具有空間整體性。

追溯香格里拉市各年高山松變化情況可知:20世紀(jì)80-90年代,當(dāng)?shù)貫榱税l(fā)展經(jīng)濟(jì)盲目新建了許多森林木材工廠企業(yè),從而加劇了森林資源的消耗,因此可解釋引起該變化的因素主要是人為干擾[16]。1997年后國家出臺全面禁止天然林采伐的政策,同時(shí)實(shí)施了退耕還林等一系列森林生態(tài)保護(hù)工程和森林保護(hù)措施,使得香格里拉高山松開始恢復(fù),從2002-2017年的研究結(jié)果也反映了該恢復(fù)情況[15,22,23]。研究結(jié)果表明,基于遙感與地統(tǒng)計(jì)學(xué)相結(jié)合的方法能夠高效地獲得AGB及其相關(guān)的遙感因子的時(shí)空變異特征,即可在一定變程A0范圍內(nèi)判斷高山松AGB的分布,以及可大體判斷受何種類型因素的影響,該方法可為日后高山松AGB的時(shí)空模擬預(yù)測和反演制圖奠定一定的理論依據(jù)。

相較于空間插值法受樣本點(diǎn)的影響[22],本研究采用經(jīng)過BSTS濾波算法重構(gòu)時(shí)間序列數(shù)據(jù)來估測反演高山松AGB提供了更可靠的數(shù)據(jù)源[3,4,11,25-26];受地統(tǒng)計(jì)方法樣本量的取值大小,以及獲取樣地點(diǎn)實(shí)測數(shù)據(jù)局限性的影響,樣地點(diǎn)的樣本量采用隨機(jī)取點(diǎn)的方式滿足研究條件。樣本量的取值多少在此類研究中暫無具體的科學(xué)依據(jù),后續(xù)對該部分內(nèi)容值得探究;本研究運(yùn)用了30 m分辨率的遙感影像數(shù)據(jù),若用分辨率更高的遙感影像數(shù)據(jù),可能會獲得更加精確的結(jié)果。