張玉娟，曲建光，王 雷

(黑龍江工程學(xué)院測繪工程學(xué)院，黑龍江 哈爾濱 150050)

0 引言

生態(tài)系統(tǒng)服務(wù)是生態(tài)系統(tǒng)形成和變化過程中提供給人類生存與發(fā)展所需要的各種效用[1-2]，是鏈接生態(tài)系統(tǒng)與社會經(jīng)濟系統(tǒng)的橋梁[3]. 生態(tài)系統(tǒng)服務(wù)能力是綜合評價生態(tài)環(huán)境優(yōu)劣和生態(tài)經(jīng)濟價值的重要指標(biāo). 對生態(tài)系統(tǒng)服務(wù)價值(ecosystem services value，ESV)進(jìn)行定量評估已成為生態(tài)學(xué)、 地理學(xué)、 環(huán)境經(jīng)濟學(xué)的熱點和前沿問題[4]. ESV分布均勻度是指不同等級的ESV在數(shù)量和空間上分布的均勻程度. 通過分析不同等級的ESV數(shù)量與空間的分布均勻度可以為研究區(qū)域景觀格局調(diào)整提供新的參考，以實現(xiàn)小代價換取區(qū)域ESV增值的目的.

1997年，Costanza在《Nature》中第一次提出了ESV估算原理和方法[5]. 近30年各國學(xué)者研究了多種ESV估算方法，其中國內(nèi)參考量最高的是謝高地團(tuán)隊[6]依據(jù)Costanza等研究的部分成果，對中國近200位生態(tài)學(xué)者進(jìn)行問卷調(diào)查，制定出中國區(qū)域ESV當(dāng)量因子表. 謝高地團(tuán)隊于2006年又對中國500位生態(tài)學(xué)者進(jìn)行問卷調(diào)查，得出新的生態(tài)系統(tǒng)服務(wù)評估體系，并于2015年進(jìn)行修訂[7]. 很多國內(nèi)學(xué)者以謝高地團(tuán)隊的研究為基礎(chǔ)，根據(jù)各自研究區(qū)域情況調(diào)整生態(tài)服務(wù)當(dāng)量因子表，估算研究區(qū)域ESV，為區(qū)域生態(tài)環(huán)境建設(shè)提供依據(jù). 其中較多研究集中在景觀格局變化對生態(tài)系統(tǒng)服務(wù)價值的影響上[8-9]. 也有部分研究集中在生態(tài)系統(tǒng)服務(wù)價值權(quán)衡協(xié)同關(guān)系上[10-11]. 近年來，在生態(tài)系統(tǒng)服務(wù)的研究上注重景觀格局變化對ESV的影響和ESV權(quán)衡協(xié)同關(guān)系，對ESV空間分布特征研究較少，對ESV的分布均勻度更是未見研究.

近20年，松花江流域哈爾濱段草地減少，鹽堿地增加，自然植被用地轉(zhuǎn)化頻繁，土地利用結(jié)構(gòu)趨于不穩(wěn)定狀態(tài)[12-13]，不斷變化的土地利用方式必然導(dǎo)致ESV時空分布的變化. 目前對該區(qū)域在景觀格局、 糧食安全和水環(huán)境評價方面有較多的研究[14-17]. 對該區(qū)域綜合ESV定量評估研究文獻(xiàn)只有1篇：龔文峰等[18]針對研究區(qū)域2006年和2005年ESV總量進(jìn)行計算，但學(xué)者并未實現(xiàn)區(qū)域小單元劃分，沒有進(jìn)行ESV空間分布變化研究，且所研究的年限較為久遠(yuǎn).

鑒于此，本文以松花江流域哈爾濱段為研究區(qū)域，以謝高地ESV模型[7]為基礎(chǔ)，修正本研究區(qū)域的ESV模型，對區(qū)域進(jìn)行網(wǎng)格劃分，通過地統(tǒng)計學(xué)及空間相關(guān)分析方法對研究區(qū)域ESV時空變化特征及空間相關(guān)性進(jìn)行研究. 引入信息熵均勻度理論和獨占圓被含均勻度理論對研究區(qū)域不同等級的ESV在數(shù)量分布和空間分布的均勻度進(jìn)行計算與變化分析. 研究成果有助于完善松花江流域哈爾濱段的生態(tài)格局，并為其它研究區(qū)域ESV分布均勻度的研究提供范例.

1 數(shù)據(jù)來源及研究方法

本文數(shù)據(jù)來源及處理過程、 獲得的松花江流域哈爾濱段1995年、 2005年和2015年景觀類型分布圖，參見參考文獻(xiàn)[16]. 考慮數(shù)據(jù)可獲取性，在計算ESV時，本文將9個轄區(qū)合并為主城區(qū).

1.1 ESV計算

1.1.1ESV計算公式

根據(jù)謝高地團(tuán)隊給出的ESV計算方法[7]，按下式計算區(qū)域ESV總價值.

(1)

式中：VEStotal為ESV總量，VCij為第i類景觀對應(yīng)第j類型服務(wù)的單位面積ESV(元·hm-2)；Si為研究區(qū)第i類景觀類型面積(hm2)；n為景觀類型個數(shù).

1.1.2單位面積ESV當(dāng)量表修訂

謝高地團(tuán)隊在2015年研究所得的單位面積ESV當(dāng)量表[7]是基于二級生態(tài)系統(tǒng)分類的，本研究所確定的生態(tài)系統(tǒng)分類為一級分類，且與謝高地團(tuán)隊的一級分類結(jié)果存在差異，需要對單位面積ESV當(dāng)量表進(jìn)行修訂. 根據(jù)本研究區(qū)域景觀特征，將其中的草地、 未利用地與謝高地團(tuán)隊二級分類結(jié)果中的草甸、 裸地相對應(yīng)，本研究區(qū)域冬季11月至次年2月為霜期，3—10月為無霜期，將謝高地團(tuán)隊當(dāng)量表中二級分類的水系與冰川積雪按3∶1加權(quán)計算本研究區(qū)域水體當(dāng)量. 通過查閱哈爾濱市1996年、 2006年和2016年統(tǒng)計年鑒，2015年林業(yè)志等參考資料，統(tǒng)計研究區(qū)域、 縣域及主城區(qū)林地主要二級類型面積比例，按二級分類面積比例加權(quán)計算各區(qū)域林地的服務(wù)價值當(dāng)量. 研究區(qū)同一區(qū)域不同年份耕地中水田與旱田面積比例變化較大，通過查閱哈爾濱市統(tǒng)計年鑒，獲得不同年份區(qū)域、 縣域及主城區(qū)耕地主要二級類型面積比例，按二級分類面積比例加權(quán)計算各區(qū)域耕地服務(wù)價值當(dāng)量. 依據(jù)哈爾濱市2013年城市綠化普查數(shù)據(jù)及Google Earth地圖數(shù)據(jù)，計算建設(shè)用地中含主要生態(tài)景觀面積比例和非生態(tài)景觀面積比例，按各類景觀對應(yīng)面積比例加權(quán)計算建設(shè)用地服務(wù)價值當(dāng)量.

1.1.3ESV當(dāng)量因子修正

ESV的變化具有時間和空間效應(yīng)，謝高地團(tuán)隊將2010年全國單位面積農(nóng)田生態(tài)系統(tǒng)糧食生產(chǎn)的凈利潤當(dāng)作1個標(biāo)準(zhǔn)當(dāng)量因子，確定中國2010年1個生態(tài)系統(tǒng)價值當(dāng)量的經(jīng)濟價值為3 406.50元·hm-2. 根據(jù)《2010年哈爾濱市統(tǒng)計年鑒》，哈爾濱市耕地地均糧食產(chǎn)量為6.059 1 t·hm-2，同期全國地均糧食產(chǎn)量為5.152 4 t·hm-2，據(jù)此對本研究區(qū)域ESV當(dāng)量因子進(jìn)行修正，修正系數(shù)為2010年哈爾濱市耕地地均糧食產(chǎn)量與全國地均糧食產(chǎn)量比值1.176，得到區(qū)域內(nèi)1個生態(tài)系統(tǒng)價值當(dāng)量的經(jīng)濟價值為4 006.04元·hm-2，為了剔出通貨膨脹引起的不同年份ESV差異，本研究所有年份ESV的計算均參照2010年當(dāng)量因子計算結(jié)果.

1.2 基于半方差函數(shù)的空間變異結(jié)構(gòu)擬合

空間變異結(jié)構(gòu)可通過半方差函數(shù)的理論模型模擬，計算公式為[19]：

(2)

式中：γ(h)為變異函數(shù)值；Z(xi)和Z(xi+h)是Z(x)在空間單元xi和xi+h上的ESV(i=1, 2, 3, …,N(h))；N(h)是分割距離h的樣本量.

最優(yōu)變異函數(shù)擬合模型的選擇主要考慮決定系數(shù)R2取得最大值，殘差值RSS取得較小值，并結(jié)合結(jié)構(gòu)方差與基臺值的比值Rc(該值反映了空間相關(guān)性的強弱，將它定義為空間相關(guān)系數(shù)，即該系數(shù)越大空間相關(guān)性越強，說明樣本間的變異更多是由結(jié)構(gòu)性引起的)以及有效變程A(變程值表示某一特征在空間上自相關(guān)的空間幅度，在大于變程的空間尺度上該變量沒有自相關(guān)性).

1.3 ESV均勻度理論

1) 基于信息熵理論的均勻度. 信息熵最早由克勞德﹒香農(nóng)提出，用于判斷不同通信信號出現(xiàn)概率的均勻性，當(dāng)不同信號出現(xiàn)概率呈現(xiàn)比較均勻分布時所能呈現(xiàn)的信息量較大，當(dāng)不同的信號出現(xiàn)概率呈現(xiàn)不均勻分布時所能呈現(xiàn)的信息量較小[20]. 其計算公式為：

(3)

式中：H為信息熵；n為不同信號個數(shù)；pi為信號i出現(xiàn)概率.

在景觀生態(tài)學(xué)中基于信息熵理論的均勻度用來描述景觀里不同生態(tài)系統(tǒng)的分配均勻程度，其計算公式為：

(4)

式中：ISHE為基于信息熵理論的均勻度指數(shù)；H為信息熵；Hmax為理論最大均勻信息.

本文引入信息熵理論，將計算獲得的ESV按照大小進(jìn)行相等間隔法等級劃分，將不同等級的ESV看成不同信號，統(tǒng)計網(wǎng)格不同等級ESV出現(xiàn)的概率，通過式(3)～(4)計算研究區(qū)域三期ESV基于信息熵理論的均勻度，該均勻度可衡量研究區(qū)域不同等級ESV在數(shù)量上分布的均勻程度，不能衡量空間分布的均勻程度.

2) 獨占圓與被含均勻度理論. 獨占圓與被含均勻度理論最早由羅傳文教授[21]提出，用于判斷林分中植物格局的均勻程度，定義為：二維空間中具有相同屬性的一組離散點Pi(i=1, 2, …,n)，假設(shè)A∈Pi(i=1, 2, …,n)，將距離點A的最近鄰體點記為B，B∈Pi(i=1, 2, …,n)，AB兩點間距離稱為點A的緊鄰距離，記為D(A)； 以A為圓心，以D(A)/2為半徑的圓稱為點A的獨占圓，獨占圓的任意外切正方形稱為點A的獨占體，則點A獨占體面積記為S(A).

被含均勻度U被定義為離散點Pi(i=1, 2, …,n)獨占體總面積與研究區(qū)域面積之比，計算公式為：

(5)

通過式(5)計算研究區(qū)域三期不同等級的ESV被含均勻度，該均勻度可衡量不同等級的ESV空間分布的均勻程度.

2 實驗結(jié)果與分析

2.1 ESV計算

2.1.1ESV當(dāng)量表

根據(jù)節(jié)1.1中單位面積ESV當(dāng)量表修訂方法，分別對主城區(qū)和6個縣進(jìn)行不同年份單位面積ESV當(dāng)量表計算. 其中林地、 草地、 建設(shè)用地、 水域和未利用地對各種生態(tài)服務(wù)類型的單位面積ESV當(dāng)量表見表1，耕地不同區(qū)縣不同年份對各種生態(tài)服務(wù)類型的單位面積ESV當(dāng)量表見表2.

表1 單位面積ESV當(dāng)量表Tab.1 Equivalent scale of ESV per unit area

表2 耕地單位面積ESV當(dāng)量表Tab.2 Equivalent scale of ESV per unit area of cultivated land

2.1.2網(wǎng)格劃分及ESV計算

按照5 km×5 km尺度將研究區(qū)域劃分成1 499個正方形單元，分別對三期按景觀類型融合后的景觀類型分布矢量、 網(wǎng)格矢量和區(qū)縣邊界矢量進(jìn)行相交運算，獲取每個斑塊所屬的網(wǎng)格、 景觀類型和區(qū)縣屬性，根據(jù)表1～2和式(1)計算斑塊ESV，再按所屬網(wǎng)格進(jìn)行合并，計算出三期網(wǎng)格ESV，結(jié)果見圖1.

圖1 5 km尺度ESV網(wǎng)格分布圖Fig.1 Distribution of ESV at 5 km scale

2.2 基于半方差函數(shù)的空間變異結(jié)構(gòu)擬合分析

在地統(tǒng)計學(xué)軟件GS+的支持下，構(gòu)造三期ESV變異函數(shù)模型，完成數(shù)據(jù)變異函數(shù)理論模型的擬合(見表3). 1995年、 2005年和2015年ESV分別為線性、 球狀和指數(shù)模型擬合效果最好，決定系數(shù)R2分別為0.990、 0.993和0.826，殘差值RSS均為最小值，結(jié)構(gòu)方差與基臺值的比值Rc均大于0.5，且呈現(xiàn)遞增的變化，說明在5 km×5 km單元網(wǎng)格劃分下，1995—2015年，結(jié)構(gòu)性因素對ESV影響程度加深[22]. 1995年、 2005年和2015年ESV空間分異的變程分別為182.624、 229.100和53.100 km，表明5 km×5 km單元網(wǎng)格劃分下研究區(qū)ESV具有高度的空間相關(guān)性.

表3 5 km尺度ESV變異函數(shù)擬合模型參數(shù)Tab.3 Parameters of ESV variation function fitting model at 5km scale

2.3 ESV時空變化分析

按照等間距法對ESV劃分為最低值ESV、 低值ESV、 較低ESV、 中等ESV、 較高ESV、 高值ESV、 最高值ESV. 運用普通克里格插值法對劃分ESV等級后的5 km尺度離散點進(jìn)行插值生成三期ESV連續(xù)空間分布圖(見圖2).

圖2 ESV連續(xù)空間分布圖Fig.2 Continuous spatial distribution of ESV

各級別ESV空間分布存在明顯的異質(zhì)性，最低值、 低值ESV主要分布在以耕地和建設(shè)用地為主的主城區(qū)、 巴彥縣、 賓縣和依蘭縣； 高值、 最高值ESV在2005年均沒有分布，1995年分布在通河縣與方正縣交界以水域景觀為主的網(wǎng)格分布區(qū)域，2015年分布在松花江位于主城區(qū)西部水域較寬區(qū)域； 較高值ESV主要分布在水域景觀與其他景觀交界處或林地景觀中心區(qū)域，中等ESV主要分布在林地耕地交界處靠近林地一側(cè)，較低ESV主要分布在林地耕地交界處靠近耕地一側(cè)或狹窄水域景觀與耕地景觀交界處. 不同年份各級別ESV分布區(qū)域差異不大，不同級別ESV分布面積受景觀類型轉(zhuǎn)換影響，最高值、 高值ESV分布較大因素是受不同年份降雨量影響，2005年研究區(qū)域降水量小，以水域為主的景觀網(wǎng)格沒有形成聚集分布，導(dǎo)致高值、 最高值ESV均沒有分布.

2.4 ESV分布均勻度變化

在進(jìn)行ESV分布均勻度計算前，根據(jù)信息熵理論和獨占圓與被含均勻度理論，按照相等間隔法對ESV進(jìn)行5個不同等級的劃分，即低值ESV、 較低ESV、 中等ESV、 較高ESV、 高值ESV.

2.4.1基于信息熵理論的均勻度變化與分析

分別對三期各等級ESV網(wǎng)格數(shù)量進(jìn)行統(tǒng)計，計算不同等級ESV在數(shù)量上分布的均勻度，結(jié)果見表4，1995年、 2005年和2015年不同等級ESV數(shù)量分布均勻度分別為0.390 0、 0.366 3和0.389 6. ESV數(shù)量分布均勻度均比較低，源于研究區(qū)域內(nèi)各級別ESV網(wǎng)格分布數(shù)量極不均衡，近50%的網(wǎng)格ESV都處于低值區(qū)域，近30%的網(wǎng)格ESV處于較低值區(qū)域，中等、 較高和高值ESV總網(wǎng)格數(shù)不足20%.

表4 三期各等級ESV網(wǎng)格數(shù)及均勻度Tab.4 Number of ESV grids at different levels and uniformity

由表1、 2可知，不同的景觀類型單位面積產(chǎn)生的ESV差異較大，水域景觀在單位面積上產(chǎn)生的服務(wù)價值最大，分別是林地、 耕地、 草地、 建設(shè)用地、 未利用地景觀單位面積產(chǎn)生的服務(wù)價值的3.850 6倍、 22.066 6倍、 7.625 8倍、 41.704 8倍和435.815 0倍. 因此，當(dāng)網(wǎng)格內(nèi)水域面積比例較大時所能產(chǎn)生的ESV也非常大，研究區(qū)域內(nèi)水域面積比例大的網(wǎng)格個數(shù)非常少，導(dǎo)致較高和高值ESV網(wǎng)格個數(shù)也非常少； 未利用地單位面積產(chǎn)生的服務(wù)價值非常小，但未利用地斑塊面積也非常小，一般在網(wǎng)格內(nèi)部所占面積非常小，對網(wǎng)格ESV影響不足以表現(xiàn)為使ESV出現(xiàn)極小值的情況； 耕地和林地為在研究區(qū)域優(yōu)勢景觀，尤其是耕地景觀作為基質(zhì)占有最大比例的面積，當(dāng)網(wǎng)格內(nèi)主要為耕地景觀時呈現(xiàn)低值ESV，網(wǎng)格內(nèi)主要為林地景觀時呈現(xiàn)較低ESV，因此低值ESV和較低ESV網(wǎng)格比例非常大. 三期不同等級ESV在數(shù)量分布均勻度總體上呈現(xiàn)先減小后增加的變化，主要由于相對于1995年和2015年，2005年研究區(qū)域內(nèi)作為廊道的水域景觀寬度較窄，沒有網(wǎng)格范圍內(nèi)水域景觀面積比例大到能使網(wǎng)格ESV值達(dá)到較高ESV和高值ESV的級別，在1995年和2015年零星分布的較高ESV和高值ESV網(wǎng)格在2005年都呈現(xiàn)為低值ESV和較低ESV，因此2005年各等級ESV分布較1995年和2015年更不均勻.

2.4.2獨占圓與被含均勻度理論變化與分析

根據(jù)獨占圓定義，以網(wǎng)格中心點坐標(biāo)為獨占圓圓心，將三期網(wǎng)格ESV等級值賦給中心點，在ArcGIS10.3中編寫代碼繪制每期不同等級ESV的獨占圓(見圖3)，計算獨占體面積及被含均勻度，結(jié)果見表5.

圖3 不同等級ESV獨占圓Fig.3 ESV exclusive circles of different levels

表5 ESV空間分布均勻度Tab.5 Spatial distribution uniformity of ESV

從表5可以看出，各等級ESV空間分布不均勻. 低值ESV空間分布均勻度均大于0.7，在研究區(qū)域內(nèi)空間分布均勻，其余各等級ESV空間分布均勻性均小于0.5，都處于不均勻分布狀態(tài)，尤其是較高、 高值ESV及1995年和2015年的中等ESV空間分布均勻性極低，源于這些等級ESV網(wǎng)格個數(shù)非常少，且較集中分布于局部范圍內(nèi). 三期ESV平均空間分布均勻度均小于0.3，表現(xiàn)為不均勻空間分布，呈現(xiàn)增加后減小的變化，主要源于2005年中等ESV雖然在數(shù)量上分布非常少，但中等ESV網(wǎng)格較離散分布于研究區(qū)域內(nèi)，因此2005年中等ESV空間分布均勻度較1995年和2005年高，致使2005年研究區(qū)域ESV空間分布均勻度高于1995年和2005年.

3 結(jié)語

研究區(qū)域ESV存在很強的空間相關(guān)性，各級別ESV空間分布存在明顯的異質(zhì)性，低值ESV主要分布在主城區(qū)、 巴彥縣、 賓縣和依蘭縣，高值ESV主要分布在以水域景觀為主的網(wǎng)格分布區(qū)域，較高值ESV主要分布在水域景觀與其他景觀交界處或林地景觀中心區(qū)域，中等ESV主要分布在林地耕地交界處靠近林地一側(cè)，較低ESV主要分布在林地耕地交界處靠近耕地一側(cè)或狹窄水域景觀與耕地景觀交界處，不同年份各級別ESV分布差異不大.

不同等級ESV在數(shù)量分布均勻度上呈現(xiàn)先減小后增加的變化，均勻度都比較低. 各等級ESV空間分布不均勻，低值ESV空間分布均勻度均大于0.7，在研究區(qū)域內(nèi)空間分布均勻，其余各等級ESV空間分布均勻性均小于0.5，都處于不均勻分布狀態(tài)，三期ESV平均空間分布均勻度呈現(xiàn)增加后減小的變化，分別為0.239 6、 0.296 7和0.231 8，均小于0.3，表現(xiàn)為不均勻空間分布.