王 鑫 戴倩倩 徐夢潔 莊舜堯
(1.中國科學院南京土壤研究所,江蘇 南京210008;2.南京農(nóng)業(yè)大學 公共管理學院,江蘇 南京210095)
我國是竹林資源最豐富的國家,竹林不僅具有巨大的經(jīng)濟效益和社會效益,而且其碳匯能力在應對氣候變化中也扮演著舉足輕重的作用[1-2]。土壤養(yǎng)分是竹林植被賴以發(fā)育的基礎(chǔ),土壤提供了植物生長必要的營養(yǎng)元素,是植被生長所需營養(yǎng)的重要來源[3-4]。對竹林土壤變異性的研究有助于改善竹林的綜合效益,從而推動竹產(chǎn)業(yè)及上下游相關(guān)產(chǎn)業(yè)的發(fā)展[5]。目前許多學者對竹林土壤屬性的空間分布特征以及竹林土壤養(yǎng)分的影響因素進行了相關(guān)研究,如高強偉等[6]對蜀南竹海景區(qū)毛竹Phyllostachys edulis 林土壤物理性質(zhì)和空間異質(zhì)性進行了分析,劉軍等[7]在杭州市臨安區(qū)和余杭區(qū)、蘇木榮等[8]在佛山市以及沈德才等[9]在東莞市也開展了類似的研究。與此同時,隨著遙感技術(shù)的快速發(fā)展,利用遙感技術(shù)提取森林信息成為獲取林業(yè)基礎(chǔ)數(shù)據(jù)的常用途徑之一[10],如許章華等[11]、官鳳英等[12]利用遙感技術(shù)提取竹林分布信息,對福建省順昌縣竹林的空間分布特征和動態(tài)變化開展了研究,李陽光等[13]基于Landsat數(shù)據(jù)提取了浙江省的竹林信息并分析了其時空演變規(guī)律。因此為了準確、高效的獲取竹林資源的變化信息,通過結(jié)合遙感技術(shù)和地理信息技術(shù),豐富數(shù)據(jù)來源,及時更新林業(yè)信息,為林業(yè)資源動態(tài)監(jiān)測和管理提供更加高效和多元化的管理措施,契合林業(yè)資源研究的前景[14]。
福建省建甌市是全國重點林業(yè)縣,2018 年森林覆蓋率達79.6%,是竹林富集區(qū)之一,不少研究人員以此為研究區(qū)域開展了竹林林地土壤的研究。黃啟堂等[15]分析了當?shù)夭煌窳至值赝寥赖睦砘再|(zhì),董蔚[16]總結(jié)了毛竹林土壤的相關(guān)屬性的特點及其空間分布特征,林振清[17]著重研究了毛竹林土壤速效氮磷鉀含量及其空間分布,Zhang等[18]關(guān)注土壤容重、礫石含量、陽離子交換量、土壤有機碳等屬性的空間變異性,鮑思屹等[19]納入時間維度,探討了近40 a 來竹林土壤有機質(zhì)的時空變異特征。這些研究通常在GIS 軟件輔助下展開,使我們深入了解當?shù)孛窳滞寥缹傩缘目臻g分布及其伴隨時間的動態(tài)變化,但是基于遙感影像的竹林土壤屬性研究卻很少,所使用的底圖中竹林的分布信息往往相對滯后,因而影響了數(shù)據(jù)的即時性。與此同時,土壤有機質(zhì)是土壤的重要組成部分,既能為植物生長提供各種營養(yǎng)成分,又在提升土壤肥力和維持土壤生態(tài)系統(tǒng)穩(wěn)定方面發(fā)揮著重要作用[20-21]。因此本文以建甌市為研究區(qū)域,基于遙感影像提取建甌市竹林分布信息并形成研究底圖,以鄉(xiāng)(鎮(zhèn))為單位進行竹林土壤采樣,對測定后的土壤數(shù)據(jù)進行分析,通過半方差分析和空間插值來揭示研究區(qū)竹林土壤有機質(zhì)空間變異特征。研究能提供更新后的竹林分布專題圖,研究結(jié)果也能為當?shù)刂窳之a(chǎn)業(yè)的管理與規(guī)劃提供參考。
建甌市位于福建省北部,是福建省面積最大、閩北人口最多的縣級市,素有“中國竹子之鄉(xiāng)”的美譽。地處東經(jīng)117°58′45″~118°57′11″,北緯26°38′54″~27°20′26″的東南沿海低山丘陵區(qū),屬亞熱帶海洋性季風氣候,年平均氣溫19 ℃、降雨量1 600~1 800 mm,日照1 612 h,無霜期286 d,氣候溫和,雨量充沛,境內(nèi)毛竹林面積、立竹量、竹材和鮮筍產(chǎn)量均居全國縣(市)之首,竹產(chǎn)業(yè)已經(jīng)成為當?shù)氐闹еa(chǎn)業(yè)[22]。
2.1.1 數(shù)據(jù)來源 本文選用的遙感影像來自于2017 年的Landsat-8 OLI 影像和同期30 m 分辨率的DEM 數(shù)據(jù),下載自中國科學院計算機網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺(http://www.gscloud.cn/)。
2.1.2 數(shù)據(jù)處理 由于遙感影像的精度會受到如衛(wèi)星姿態(tài)、衛(wèi)星運行軌道和大氣反射等因素的影響,因此在進行遙感解譯前,需要對原始影像進行大氣校正、幾何精校正、圖像裁剪和波段融合等預處理[23],其中要用到的方法有臨近點插值法[24]、直方圖匹配法[25]、卷積算法[26]。
對于預處理后的遙感影像,需要篩選信息提取的最佳波段組合,本文使用最佳指數(shù)法[27](OIF)進行篩選,通過計算出影像數(shù)據(jù)各個波段間的標準差,構(gòu)建相關(guān)關(guān)系矩陣。在相關(guān)關(guān)系矩陣的基礎(chǔ)上,計算出三波段組合的OIF 值。OIF 值越大,說明相應組合獲得的圖像信息量就越大,組合效果趨于優(yōu)化。
篩選出竹林提取的最佳波段組合后,再采用最大似然分類法(Maximum Likelihood Classification)[28]進行監(jiān)督分類。最大似然分類法基于貝葉斯準則,假設(shè)遙感影像的每個波段數(shù)據(jù)都呈正態(tài)分布,通過訓練樣本計算出其平均值及方差、協(xié)方差等特征參數(shù),從而得出總體的先驗概率密度函數(shù),然后將待分類像元的均值、協(xié)方差等特征參數(shù)帶入判別函數(shù)計算其屬于各類別的概率,最后將像元分到歸屬概率最大的類別中[29]。為了驗證該方法的有效性,我們還采用了支持向量機法進行分類[30]。遙感數(shù)據(jù)預處理等程序利用ENVI 5.3 和ERADS 9.2 完成。
2.2.1 采樣與測定方法 2018 年4 月在研究區(qū)內(nèi)以鄉(xiāng)鎮(zhèn)為單位,根據(jù)竹林面積大小共布設(shè)了209個竹林土壤采樣點,研究人員在采樣時利用手持GPS 確定采樣預設(shè)點的位置,并導出其經(jīng)緯度坐標,然后采用五點采樣法在以預設(shè)點為中心的100 m×100 m 矩形4 個頂點和中心點處各采集一份0~20 cm 竹鞭分布層土樣,混合均勻后帶回實驗室風干處理。土壤有機質(zhì)采用低溫外熱重鉻酸鉀氧化-比色法測定[31]。
2.2.2 數(shù)據(jù)處理 (1)描述性統(tǒng)計分析 對研究區(qū)209 個樣點的土壤屬性數(shù)據(jù)進行異常值剔除,對剩余的205 個樣點統(tǒng)計各屬性的最大值、最小值、平均值、中位數(shù)和變異系數(shù)等指標,統(tǒng)計分析利用SPSS 22.0 完成。
(2)地統(tǒng)計學方法 對于服從近似正態(tài)分布的土壤屬性數(shù)據(jù)進行半方差函數(shù)擬合,半方差函數(shù)包括球狀 (Spherical) 、高斯 (Gaussian) 、指數(shù) (Exponential) 和線性 (Linear) 等模型[32],通常根據(jù)決定系數(shù)(R2)最大,殘差(RSS)最小原則確定最佳擬合模型[33],并獲取塊金值、基臺值、變程等參數(shù),以進一步描述區(qū)域化變量空間變異性的強弱、結(jié)構(gòu)等;以此為基礎(chǔ),采用克里格方法對土壤屬性進行空間插值,生成土壤屬性的空間分布圖[34]。我們還采用了一種確定性插值方法——反距離權(quán)重法(IDW),以插值點與樣本點的距離為權(quán)重進行插值[35]。半方差函數(shù)及理論模型擬合利用地統(tǒng)計學軟件GS+9.0,克里格插值在ArcGIS 10.2 下完成。
(3)精度檢驗 利用ArcGIS 10.2 軟件從原始樣點中隨機抽取20%作為獨立樣本驗證點,剩下的80%作為實驗數(shù)據(jù)參與空間插值過程。通過計算平均誤差(ME)、平均絕對誤差(MAE)和均方根誤差(RMSE) 來評價模型的預測精度,其中ME、MAE 和RMSE 越小,預測精度越高[36]。
表1 各波段組合OIF 指數(shù)Table 1 The OIF index table of band combination
在ERADS 9.2 軟件支持下,選取同期Google Earth 影像(2 m)分辨率作為標準地圖,結(jié)合實地野外調(diào)查,共選擇了15 個地面點為控制點。采用臨近點插值法對圖像進行重采樣操作,對重采樣后影像進行幾何精糾正操作,再對糾正后的單波段原始影像進行精度檢驗[37]。通過精度檢驗后,使用建甌市行政邊界依次裁剪三幅影像各波段文件,用直方圖匹配方法對圖像進行勻色處理。圖像增強則采用卷積算法,以3×3 分析窗口進行,得到預處理后的影像以提取竹林分布信息。
B1-B7 波段中,將7 個波段進行三三組合,計算三波段組合的OIF 指數(shù),結(jié)果如表1 所示,根據(jù)各組合的OIF 指數(shù)以及典型地物在各個波段的亮度均值,篩選出最佳波段組合方案為356、236 和367,結(jié)果如圖1 所示。
建甌市竹林分布情況較為復雜,竹林總面積近8.93×105hm2,遍布10 個鎮(zhèn)、4 個鄉(xiāng)和4 個街道。為了更高精度進行竹林的判讀和提取,采用最大似然分類法。按照監(jiān)督分類的一般步驟[38],首先將研究區(qū)內(nèi)地物劃分為竹林、水體、非竹林和其他四類。通過直觀觀察和輔助驗證點判定發(fā)現(xiàn),在356波段組合處理下竹林為綠色明顯區(qū)別于236 和367波段組合,因此以356 波段組合為竹林分類提取的波段組合,以356 波段組合圖像作為分類基礎(chǔ)圖像,每個大類選擇至少100 個以上訓練區(qū),合成分類地物色彩特征,構(gòu)建地物分類模板。利用混淆矩陣對分類模板的精度進行檢驗,如表2 所示竹林解譯的Kappa 系數(shù)為89.40%[39-40],滿足竹林信息提取要求。其Kappa 系數(shù)為88.78%,因此仍以最大似然分類法分類結(jié)果為基準,提取出竹林分布專題圖(圖2)。
圖1 各波段組合Fig.1 Band combination map
表2 最大似然分類混淆矩陣Table2 The confusion matrix of maximum likelihood classification
圖2 建甌市竹林分布Fig.2 Bamboo forest distribution in Jian'ou city
3.2.1 樣本描述性統(tǒng)計分析 由建甌市毛竹林土壤屬性描述性分析表(表3)可知,建甌市竹林土壤有機質(zhì)含量變幅為7.22~95.27 g/kg ,平均值為39.47 g/kg,最大值與最小值相差接近13 倍。變異系數(shù)(CV)可反映土壤特性的空間變異程度[41]。從表3 可以看出,有機質(zhì)屬于中等變異水平(CV 介于10%和100%之間)。
3.2.2 土壤有機質(zhì)的趨勢面特征 由于全局趨勢的存在,會對半變異分析過程產(chǎn)生影響,因此在半變異分析過程中要對全局趨勢進行剔除[42]。從建甌市采樣點土壤有機質(zhì)數(shù)據(jù)的趨勢分析(圖3)可以看出,在東西方向上二階趨勢不明顯,在南北方向呈“U”型變化。
3.2.3 土壤有機質(zhì)的空間結(jié)構(gòu)特征 數(shù)據(jù)的正態(tài)性檢驗結(jié)果表明,土壤有機質(zhì)不符合正態(tài)分布,在對數(shù)據(jù)進行對數(shù)轉(zhuǎn)換后,通過正態(tài)性檢驗,符合半方差分析的要求。半方差函數(shù)擬合的結(jié)果表明,土壤有機質(zhì)的最優(yōu)擬合模型是指數(shù)模型,其參數(shù)見表4。
塊金值(C0)反映由實驗誤差和小于實驗取樣尺度引起的變異,如果塊金值較大,則表明較小尺度上的某種過程不可忽視[43]。基臺值(C0+C)表示系統(tǒng)屬性的最大變異,塊金值和基臺值的比值(C0/C0+C)稱為塊基系數(shù),它反映的是變量的空間相關(guān)程度[44]。塊基系數(shù)小于25%說明變量的空間相關(guān)性強,土壤屬性受到結(jié)構(gòu)性因素(如氣候、土壤母質(zhì)、地形等)的影響較大,而受到隨機因素(如施肥、灌溉、耕作制度等)的影響較??;系數(shù)在25%~75%之間說明變量具有中等的空間相關(guān)性,土壤屬性受到隨機因素和結(jié)構(gòu)性因素的共同影響;系數(shù)大于75%,說明空間相關(guān)性很弱,土壤屬性受到結(jié)構(gòu)性因素的影響較大[45-46]。變程(A0)指區(qū)域化變量在空間上的最大相關(guān)距離[47]。由表4 可知研究區(qū)竹林土壤的有機質(zhì)表現(xiàn)出中等的空間相關(guān)性,其空間相關(guān)性范圍為24 990 m。
3.2.4 土壤有機質(zhì)的空間分布特征 在ArcGIS 10.2 軟件平臺下,應用普通克里格法和反距離權(quán)重法進行插值,插值精度檢驗結(jié)果見表5。由表5可知,兩種插值方法的平均誤差和平均絕對誤差雖然相近,但普通克里格法的均方根誤差明顯小于反距離權(quán)重法,表明普通克里格的插值結(jié)果更為精確。
以竹林分布圖為掩膜對插值后的數(shù)據(jù)進行提取和重分類,如圖4 展示了研究區(qū)內(nèi)竹林土壤有機質(zhì)的空間分布特征。插值后的竹林土壤中的有機質(zhì)變幅為22.53~76.37 g/kg,較樣本點數(shù)據(jù)變幅有所縮減。這是由于普通克里格插值的平滑效應以及將插值結(jié)果轉(zhuǎn)換為柵格數(shù)據(jù)所致[48]。其中東北部的龍村鄉(xiāng)和川石鄉(xiāng)以及西南部的南雅鎮(zhèn)竹林土壤有機質(zhì)較高;中部的芝山街道、建安街道和甌寧街道竹林土壤有機質(zhì)含量較低。土壤有機質(zhì)含量介于22.53~32.88 g/kg 的竹林面積最大,占比35.95%,土壤有機質(zhì)含量介于60.97~76.37 g/kg 的竹林面積最小,僅占比7.11%。
表3 土壤有機質(zhì)描述性分析Table 3 Descriptive statistics of SOM
圖3 建甌市竹林土壤有機質(zhì)趨勢分析Fig.3 Trend analysis of SOM in Jian'ou city
表4 半方差理論模型與參數(shù)Table 4 Semivariogram theoretical model andparameters
圖4 建甌市竹林土壤有機質(zhì)空間分布Fig.4 Space distribution map of SOM in Jian'ou city
4.1 本研究采用遙感影像解譯的方式,以及時準確地更新竹林分布信息。遙感影像信息提取的方法眾多,包括平行六面體分類法[49]、最小距離分類法[50]、馬氏距離分類法[51]、神經(jīng)網(wǎng)絡(luò)分類法[52]、最大似然分類法[53]及支持向量機法[54]等。而研究人員[55]基本認同最大似然分類法的效果較好,在Kappa 系數(shù)和總體分類精度上往往優(yōu)于其他分類方法的分類效果,因而應用較為廣泛。如白宇興[56]在提取西安市未央?yún)^(qū)土地利用分類信息時分別使用了最大似然分類法和支持向量機法進行了解譯,前者的整體分類精度(97.80%)和Kappa 系 數(shù)(97.07%) 均 高 于 后 者(91.78%、88.81%)。本文在竹林信息提取時選用了最大似然分類法和支持向量機法進行解譯,其中基于最大似然分類法的信息提取精度達到了89.4%,分類精度也高于支持向量機法,驗證了這一觀點。
就土壤屬性而言,董蔚[16]測得的建甌市毛竹林土壤有機質(zhì)均值為39.8 g/kg,其空間分布以西南及東北區(qū)為高;林振清[17]的文獻中土壤有機質(zhì)含量均值為39.8 g/kg,呈從中部向西南及東北方向逐漸升高的趨勢;Zhang 等[18]得到的有機質(zhì)均值為40.8 g/kg,中部地區(qū)有機質(zhì)含量最低。本研究測得的結(jié)果與他們較為相近,土壤有機質(zhì)變幅為7.2~95.3 g/kg,均值為39.5 g/kg,其中西南部以及東北部的鄉(xiāng)鎮(zhèn)(街道)有機質(zhì)較高,而中部的有機質(zhì)含量較低。中間區(qū)域的海拔相對較低,通常海拔越高,有機質(zhì)含量往往越高[29]。傅平[57]于2008 年的研究中指出建甌市耕地土壤有機質(zhì)含量平均值為28.5 g/kg,屬中等偏上水平,其中最大值為85.9 g/kg,最小值為2 g/kg。本研究結(jié)果要明顯高于該值,這與竹林土壤的特性有關(guān)。竹林土壤之上覆蓋有森林凋落物,在降水、森林土壤動物和微生物等因素的作用下,轉(zhuǎn)移到土壤中,因此林地土壤有機質(zhì)含量普遍要高于農(nóng)田。本研究測定的樣本有機質(zhì)及其變幅與已有文獻基本一致,體現(xiàn)了數(shù)據(jù)的合理性,在此基礎(chǔ)上通過插值得到的空間分布特征也呈現(xiàn)相同的特點,也體現(xiàn)了研究結(jié)果的合理性。
4.2 本文以建甌市為研究區(qū)域,基于遙感影像提取建甌市竹林分布信息并形成研究底圖,在此基礎(chǔ)上以鄉(xiāng)(鎮(zhèn))為單位,進行竹林土壤采樣,通過半方差分析和空間插值來揭示研究區(qū)竹林土壤有機質(zhì)的空間變異特征。最佳指數(shù)法的計算結(jié)果表明,356 波段為竹林提取的最佳組合方式。在此基礎(chǔ)上采用最大似然法得到的竹林解譯信息精度較高。利用遙感影像來提取竹林信息,一方面可以獲取及時更新的竹林分布專題圖,另一方面也令后續(xù)的空間插值結(jié)果更為客觀、準確。插值分析表明,研究區(qū)竹林土壤有機質(zhì)在不同方向的空間分布存在差異。由于土壤有機質(zhì)受到結(jié)構(gòu)性和隨機因素的共同影響,即自然條件和人為因素的綜合影響,因此建甌市竹林的管理和規(guī)劃應關(guān)注人為因素的影響[58]。建議在今后的研究中納入其他土壤屬性并探討其空間分布和變異特征,同時注重遙感數(shù)據(jù)的應用與更新,從而更好的推動當?shù)刂窳之a(chǎn)業(yè)的發(fā)展。