劉長峰,侯 鷹,陳衛(wèi)平,崔昊天

1 中國科學(xué)院生態(tài)環(huán)境研究中心城市與區(qū)域生態(tài)國家重點(diǎn)實(shí)驗(yàn)室, 北京 100085 2 中國科學(xué)院大學(xué), 北京 100049

中國在過去40年經(jīng)歷了快速城市化的過程,在推動(dòng)各個(gè)地區(qū)社會(huì)經(jīng)濟(jì)快速發(fā)展的同時(shí),也對區(qū)域生態(tài)環(huán)境產(chǎn)生了諸多不利影響[1- 3],如人類生產(chǎn)生活產(chǎn)生的大量污染物損害了生態(tài)系統(tǒng)的結(jié)構(gòu)、功能和健康[4],城鎮(zhèn)建設(shè)用地?cái)U(kuò)張導(dǎo)致城市周邊的森林、草地、濕地、農(nóng)田等自然或半自然生態(tài)系統(tǒng)面積下降,景觀斑塊破碎,并造成生態(tài)系統(tǒng)服務(wù)損失[5- 6]。城市化導(dǎo)致的生態(tài)系統(tǒng)的不利變化目前已成為學(xué)界關(guān)注的熱點(diǎn)問題[7]。

城市化的生態(tài)風(fēng)險(xiǎn)是指城市經(jīng)濟(jì)發(fā)展、城區(qū)建設(shè)和城市人類生產(chǎn)生活活動(dòng)所導(dǎo)致的生態(tài)系統(tǒng)等的不利變化或?qū)θ祟惤】翟斐傻牟涣加绊懙某潭群涂赡苄訹8]。過去對于城市化的生態(tài)風(fēng)險(xiǎn)研究主要集中在土地利用變化和環(huán)境污染兩個(gè)方面:(1)以區(qū)域景觀格局作為生態(tài)風(fēng)險(xiǎn)評價(jià)終點(diǎn),基于景觀指數(shù)和土地利用類型得分構(gòu)建生態(tài)風(fēng)險(xiǎn)指數(shù),用以評估土地利用變化帶來的生態(tài)系統(tǒng)受損的風(fēng)險(xiǎn)[9- 12]；(2)城市化過程中大量污染物的排放影響人體健康和區(qū)域生態(tài)系統(tǒng)特定功能的風(fēng)險(xiǎn)[13- 15]。以上研究的風(fēng)險(xiǎn)評價(jià)終點(diǎn)往往只是生態(tài)系統(tǒng)的某一種屬性或過程,不能綜合地反映城市化導(dǎo)致的區(qū)域整體的生態(tài)風(fēng)險(xiǎn)水平。

生態(tài)系統(tǒng)服務(wù)能夠反映生態(tài)系統(tǒng)的整體性,使評價(jià)者能夠明確保護(hù)對象和屬性,其貨幣形式的評價(jià)結(jié)果便于同決策者交流,是理想的生態(tài)風(fēng)險(xiǎn)評價(jià)終點(diǎn)[16-17]。近年來以生態(tài)系統(tǒng)服務(wù)作為生態(tài)風(fēng)險(xiǎn)評價(jià)終點(diǎn)的研究逐漸增多,Thomsen 等[18]基于土壤生態(tài)系統(tǒng)服務(wù)構(gòu)建了土壤污染的生態(tài)風(fēng)險(xiǎn)評估概念模型；Munns等[19]提出了以生態(tài)系統(tǒng)服務(wù)的損失耦合場地生態(tài)風(fēng)險(xiǎn)評估和自然資源破壞評估的概念框架；曹祺文等[20]提出了基于生態(tài)系統(tǒng)服務(wù)的景觀生態(tài)風(fēng)險(xiǎn)評價(jià)框架。實(shí)例研究方面,Kang等[21]耦合區(qū)域生態(tài)系統(tǒng)服務(wù)價(jià)值和生態(tài)系統(tǒng)健康評估了京津冀城市群的生態(tài)風(fēng)險(xiǎn)；Pártl等[22]沿用傳統(tǒng)風(fēng)險(xiǎn)定義從生態(tài)系統(tǒng)服務(wù)的供給脆弱度和風(fēng)險(xiǎn)源的脅迫程度兩個(gè)方面,評估了捷克共和國的綜合生態(tài)風(fēng)險(xiǎn)水平。然而,目前基于生態(tài)系統(tǒng)服務(wù)的生態(tài)風(fēng)險(xiǎn)評價(jià)研究仍處于起步階段,有關(guān)風(fēng)險(xiǎn)表征方法的研究還很缺乏。此外,由于缺乏同生態(tài)系統(tǒng)服務(wù)相關(guān)的風(fēng)險(xiǎn)閾值,使得傳統(tǒng)的風(fēng)險(xiǎn)商表征方法不能適用。因此,大多數(shù)現(xiàn)有研究還停留在評價(jià)框架或概念模型階段,實(shí)例研究較少。

針對已有研究的不足,本研究旨在建立適用于城市化區(qū)域的基于生態(tài)系統(tǒng)服務(wù)的生態(tài)風(fēng)險(xiǎn)表征方法。將多種類型的生態(tài)系統(tǒng)服務(wù)納入到生態(tài)風(fēng)險(xiǎn)評估中,參考TOPSIS (Technique for Order Preference by Similarity to an Ideal Solution) 模型(優(yōu)劣解距離法),以空間上各類生態(tài)系統(tǒng)服務(wù)價(jià)值的最大損失量作為風(fēng)險(xiǎn)參考值,構(gòu)建了生態(tài)風(fēng)險(xiǎn)綜合指數(shù),建立了以生態(tài)系統(tǒng)服務(wù)為評價(jià)終點(diǎn)的生態(tài)風(fēng)險(xiǎn)表征方法,并使用構(gòu)建的方法分析了北京市的生態(tài)風(fēng)險(xiǎn)；最后使用基于蒙特卡洛模擬的不確定性分析和參數(shù)敏感性分析對生態(tài)風(fēng)險(xiǎn)指數(shù)的可靠性進(jìn)行了分析。

1 研究方法

1.1 生態(tài)風(fēng)險(xiǎn)評估框架

生態(tài)風(fēng)險(xiǎn)評估框架主要包括問題形成、風(fēng)險(xiǎn)分析和風(fēng)險(xiǎn)表征三個(gè)過程。本研究關(guān)注的是城市化區(qū)域生態(tài)系統(tǒng)服務(wù)下降的風(fēng)險(xiǎn),城鎮(zhèn)土地?cái)U(kuò)張是該類區(qū)域中城市化過程最顯著的特征。過去很多研究表明土地利用變化是導(dǎo)致生態(tài)系統(tǒng)服務(wù)變化的最主要因素[6,23],因此城市化區(qū)域的風(fēng)險(xiǎn)脅迫主要是城市化帶來的土地利用等人類活動(dòng),相應(yīng)的風(fēng)險(xiǎn)評價(jià)終點(diǎn)是區(qū)域生態(tài)系統(tǒng)服務(wù),服務(wù)下降是風(fēng)險(xiǎn)脅迫產(chǎn)生的不利生態(tài)效應(yīng)(圖1)。

圖1 基于生態(tài)系統(tǒng)服務(wù)的城市化區(qū)域生態(tài)風(fēng)險(xiǎn)評價(jià)基本框架Fig.1 A basic ecological risk assessment framework for urbanized areas based on ecosystem services

1.2 生態(tài)風(fēng)險(xiǎn)表征方法

傳統(tǒng)的評價(jià)環(huán)境污染物導(dǎo)致的生態(tài)風(fēng)險(xiǎn)是通過量化風(fēng)險(xiǎn)脅迫的程度(如污染物濃度)或危害程度(如健康損害指數(shù))來表征生態(tài)風(fēng)險(xiǎn)水平[13- 15]。本研究中風(fēng)險(xiǎn)脅迫(土地利用)與風(fēng)險(xiǎn)評價(jià)終點(diǎn)(生態(tài)系統(tǒng)服務(wù))之間的關(guān)系明確,因此以生態(tài)系統(tǒng)服務(wù)損失量來表征生態(tài)風(fēng)險(xiǎn)(圖1)。由于缺乏關(guān)于生態(tài)系統(tǒng)服務(wù)損失的風(fēng)險(xiǎn)閾值,本研究提出以生態(tài)系統(tǒng)服務(wù)價(jià)值量的相對損失量來表征生態(tài)風(fēng)險(xiǎn),構(gòu)建適用于城市化區(qū)域的生態(tài)風(fēng)險(xiǎn)綜合指數(shù)。

生態(tài)風(fēng)險(xiǎn)綜合指數(shù)參考了TOPSIS模型[24]。首先找出區(qū)域內(nèi)所有空間樣本(柵格像元或矢量多邊形區(qū)域等)中各類生態(tài)系統(tǒng)服務(wù)價(jià)值達(dá)到最大值和最小值的兩個(gè)樣本作為理想樣本點(diǎn)(正理想樣本點(diǎn)和負(fù)理想樣本點(diǎn)),以正理想樣本點(diǎn)的服務(wù)價(jià)值量作為風(fēng)險(xiǎn)閾值,以這兩個(gè)理想樣本點(diǎn)的歐氏距離作為該研究區(qū)域生態(tài)系統(tǒng)服務(wù)損失達(dá)到最大的表征,則綜合風(fēng)險(xiǎn)指數(shù)定義為每一個(gè)真實(shí)樣本與正理想樣本點(diǎn)間的歐氏距離與生態(tài)系統(tǒng)服務(wù)損失達(dá)到最大的歐氏距離之比。具體計(jì)算方法如下:

(1)基于現(xiàn)有的某種生態(tài)系統(tǒng)服務(wù)價(jià)值計(jì)算方法計(jì)算出m個(gè)空間樣本(柵格像元或矢量多邊形區(qū)域等)的n類生態(tài)系統(tǒng)服務(wù)價(jià)值,提取所有空間樣本構(gòu)建生態(tài)系統(tǒng)服務(wù)價(jià)值矩陣ESV=(esvij)mn。

(1)

(2)為了消除量綱和數(shù)量級的影響,對生態(tài)系統(tǒng)服務(wù)價(jià)值矩陣ESV進(jìn)行歸一化處理(式2),之后計(jì)算加權(quán)樣本矩陣Y=(yij)mn(式3)。

(2)

式中,rij:樣本i的第j種生態(tài)系統(tǒng)服務(wù)價(jià)值的歸一化數(shù)值；max(esv·j):第j種生態(tài)系統(tǒng)服務(wù)價(jià)值所有樣本中的最大值；min(esv·j):第j種生態(tài)系統(tǒng)服務(wù)價(jià)值所有樣本中的最小值。

(3)

式中,yij:樣本i的第j種生態(tài)系統(tǒng)服務(wù)價(jià)值的加權(quán)歸一化數(shù)值；ωj:第j種生態(tài)系統(tǒng)服務(wù)價(jià)值的權(quán)重,對不同生態(tài)系統(tǒng)服務(wù)類型的賦權(quán)體現(xiàn)了在生態(tài)風(fēng)險(xiǎn)評價(jià)中對不同生態(tài)系統(tǒng)服務(wù)的側(cè)重不同。

(3)正理想樣本點(diǎn)和負(fù)理想樣本點(diǎn)的識(shí)別

(4)

(5)

式中,y+:正理想樣本,即所有的生態(tài)系統(tǒng)服務(wù)價(jià)值量均達(dá)到最大值的理想樣本點(diǎn)；y-:負(fù)理想樣本,即所有的生態(tài)系統(tǒng)服務(wù)價(jià)值量均達(dá)到最小值的理想樣本點(diǎn)；max(y·j):生態(tài)系統(tǒng)服務(wù)j在所有區(qū)域中加權(quán)歸一化最大值；min(y·j):生態(tài)系統(tǒng)服務(wù)j在所有區(qū)域中加權(quán)歸一化最小值。

(4)生態(tài)系統(tǒng)服務(wù)價(jià)值相對損失量的表征:

(6)

(7)

(5)風(fēng)險(xiǎn)指數(shù)的計(jì)算:

(8)

式中,Ri:樣本i的生態(tài)風(fēng)險(xiǎn)指數(shù),數(shù)值介于0到1之間。越接近于0表明i的生態(tài)系統(tǒng)服務(wù)相對損失量越小,即生態(tài)風(fēng)險(xiǎn)越低；越接近于1表明i的生態(tài)系統(tǒng)服務(wù)相對損失量越大,即生態(tài)風(fēng)險(xiǎn)最高。

(6)風(fēng)險(xiǎn)等級劃分

為了減少風(fēng)險(xiǎn)等級劃分的主觀性,通過對樣本性質(zhì)進(jìn)行探索,找到相應(yīng)的風(fēng)險(xiǎn)等級分界點(diǎn)。由于構(gòu)建的風(fēng)險(xiǎn)指數(shù)只有單一屬性列表,故采用高斯混合模型對數(shù)據(jù)的聚類數(shù)目進(jìn)行探索[25]。如果所有風(fēng)險(xiǎn)指數(shù)樣本呈多峰分布且每一個(gè)峰均服從正態(tài)分布,則可以使用每一個(gè)正態(tài)分布峰的均值或分位點(diǎn)作為風(fēng)險(xiǎn)等級的分界線；如果不是正態(tài)分布,則使用Jenks自然斷點(diǎn)分類法劃分生態(tài)風(fēng)險(xiǎn)等級[26]。

2 案例研究

2.1 研究區(qū)概況

近幾十年來北京市的城市化發(fā)展迅速,城市區(qū)域快速擴(kuò)張,現(xiàn)已形成五條環(huán)城公路。一般認(rèn)為五環(huán)以內(nèi)的區(qū)域?yàn)楸本┑氖袇^(qū),五、六環(huán)之間的區(qū)域?yàn)楸本┑慕紖^(qū)；六環(huán)路以外的區(qū)域?yàn)楸本┑霓r(nóng)村地區(qū)。北京市的人工地表主要集中于六環(huán)以內(nèi),人工地表率隨著離城市中心的距離增加而明顯下降(圖2)。因此,北京市是典型的快速城市化區(qū)域,是研究快速城市化發(fā)展導(dǎo)致的區(qū)域生態(tài)風(fēng)險(xiǎn)的理想?yún)^(qū)域。

圖2 2015年北京市土地利用類型分布圖Fig.2 The land use map of Beijing in 2015

2.2 生態(tài)系統(tǒng)服務(wù)價(jià)值核算方法及權(quán)重的確定

生態(tài)系統(tǒng)服務(wù)價(jià)值核算方法主要有基于生態(tài)系統(tǒng)服務(wù)功能量的價(jià)值定量方法[27- 29]和當(dāng)量因子法[30-34]。本研究使用當(dāng)量因子法進(jìn)行生態(tài)系統(tǒng)服務(wù)價(jià)值核算,當(dāng)量因子法所需參數(shù)和數(shù)據(jù)量較少,易于在區(qū)域尺度核算各類生態(tài)系統(tǒng)服務(wù)價(jià)值[23,35-36]?；诋?dāng)量因子法的生態(tài)系統(tǒng)服務(wù)價(jià)值量計(jì)算如下:

(9)

式中,E:1 個(gè)標(biāo)準(zhǔn)當(dāng)量因子的生態(tài)系統(tǒng)服務(wù)價(jià)值量 (元/hm2),指 1 hm2全國平均產(chǎn)量的農(nóng)田每年自然糧食產(chǎn)量的經(jīng)濟(jì)價(jià)值；Aik:區(qū)域i內(nèi)土地利用類型k的總面積(hm2),本研究用到6種土地利用類型:森林、草地、耕地、濕地、人工地表和裸地；Djk:土地利用類型k所提供的生態(tài)系統(tǒng)服務(wù)j的當(dāng)量因子,本研究參考了謝高地等制定的適用于全國陸地生態(tài)系統(tǒng)服務(wù)的當(dāng)量因子表[34](表1),根據(jù)北京市的生態(tài)系統(tǒng)服務(wù)當(dāng)量因子是全國平均水平的1.04倍,對應(yīng)用于北京市的當(dāng)量因子進(jìn)行修正[31]。

本研究選取了適用于北京市的9類生態(tài)系統(tǒng)服務(wù),使用熵權(quán)法對服務(wù)的權(quán)重進(jìn)行計(jì)算[37- 38](表1)。

表1 生態(tài)系統(tǒng)服務(wù)單位面積當(dāng)量因子和權(quán)重

2.3 不確定分析與參數(shù)敏感分析

由于本研究提出的生態(tài)風(fēng)險(xiǎn)指數(shù)無法使用實(shí)測數(shù)據(jù)驗(yàn)證,因此使用不確定性分析和參數(shù)敏感性分析對風(fēng)險(xiǎn)指數(shù)的可靠性進(jìn)行評價(jià)。如果模型(本研究中為生態(tài)風(fēng)險(xiǎn)指數(shù))計(jì)算結(jié)果的不確定性太大,說明模型的可靠性低[39- 40]；反之,如果模型模擬結(jié)果的不確定性可接受,則說明模型較為可靠。模型參數(shù)的敏感性分析是用于評價(jià)參數(shù)對模型輸出結(jié)果影響程度的方法[41- 42],通過參數(shù)敏感性分析,可以找到對模型計(jì)算結(jié)果影響較大的參數(shù),從而有針對性地對模型進(jìn)行優(yōu)化。

為提高風(fēng)險(xiǎn)指數(shù)可靠性評價(jià)的效率,不確定分析與參數(shù)敏感分析在模擬柵格中進(jìn)行。模擬柵格的土地利用情況為森林、草地、耕地、濕地和人工地表各占20%,即不考慮柵格內(nèi)部土地利用類型比例的不確定性,針對9種生態(tài)系統(tǒng)服務(wù)的權(quán)重和當(dāng)量因子(當(dāng)量因子的不確定性綜合考慮5類土地利用類型當(dāng)量因子之和)進(jìn)行不確定分析與參數(shù)敏感分析。

2.3.1不確定性分析

本研究使用蒙特卡羅模擬的方法進(jìn)行不確定性分析。蒙特卡羅模擬是通過從模型參數(shù)的概率分布中隨機(jī)抽樣,計(jì)算出一系列運(yùn)算結(jié)果樣本的過程[43]。本研究使用蒙特卡洛模擬結(jié)果樣本的變異系數(shù)、最大變異范圍和相對變異指數(shù)來表征風(fēng)險(xiǎn)指數(shù)計(jì)算結(jié)果的不確定性,三個(gè)指標(biāo)值越小,說明風(fēng)險(xiǎn)指數(shù)計(jì)算結(jié)果的不確定性越小。

(1)變異系數(shù)(coefficient of variance,CV)

(10)

式中,μ為樣本均值；s為樣本標(biāo)準(zhǔn)差。變異系數(shù)是衡量一系列數(shù)據(jù)變異性的常用指標(biāo),適用于任何概率分布。

(2)最大變異范圍

(3)相對變異指數(shù)(relative variability index,RVI)

2.3.2參數(shù)敏感性分析

本研究采用局部參數(shù)敏感性分析,檢驗(yàn)每個(gè)參數(shù)值增加50%對風(fēng)險(xiǎn)指數(shù)計(jì)算結(jié)果的影響[42],計(jì)算每個(gè)參數(shù)的敏感性在全部參數(shù)的敏感性中的貢獻(xiàn)度,計(jì)算公式如下:

(11)

式中,S(pi)為參數(shù)pi的敏感性貢獻(xiàn)度,以某個(gè)參數(shù)變化導(dǎo)致風(fēng)險(xiǎn)指數(shù)計(jì)算結(jié)果的變化量絕對值占所有參數(shù)變化導(dǎo)致風(fēng)險(xiǎn)指數(shù)計(jì)算結(jié)果的變化量絕對值的總和的比例來表示；pi為參數(shù)i的實(shí)際取值；f(·)為模型運(yùn)算公式。

2.4 數(shù)據(jù)與分析

本研究主要使用的數(shù)據(jù)為:2015年北京市土地利用分類30 m分辨率柵格數(shù)據(jù)[44],北京市行政邊界數(shù)據(jù)(http://www.resdc.cn/Default.aspx),生態(tài)系統(tǒng)服務(wù)價(jià)值當(dāng)量因子表[34]。 在ArcMap 10.5中將北京市劃分為1 km×1 km的網(wǎng)格,然后使用局域統(tǒng)計(jì)工具提取出每個(gè)網(wǎng)格的土地利用配置；將數(shù)據(jù)導(dǎo)出后,在MATLAB 2017b中執(zhí)行生態(tài)風(fēng)險(xiǎn)指數(shù)的運(yùn)算以及蒙特卡羅模擬；使用R 3.4.3的mixtools工具包進(jìn)行高斯混合模型訓(xùn)練以及數(shù)據(jù)的分布檢驗(yàn)。使用ArcMap 10.5將生態(tài)風(fēng)險(xiǎn)計(jì)算結(jié)果進(jìn)行空間化表達(dá)；使用MATLAB 2017b繪制風(fēng)險(xiǎn)指數(shù)頻率分布直方圖和蒙特卡羅模擬結(jié)果圖；使用Excel 2016繪制參數(shù)敏感性分析結(jié)果圖。

3 結(jié)果與討論

3.1 生態(tài)風(fēng)險(xiǎn)等級劃分結(jié)果

2015年北京市生態(tài)風(fēng)險(xiǎn)指數(shù)空間樣本的頻率分布直方圖顯示風(fēng)險(xiǎn)指數(shù)呈多峰分布(圖3)。經(jīng)過高斯混合模型的訓(xùn)練,將所有樣本分為4個(gè)峰的貝葉斯信息準(zhǔn)則值(Bayesian Information Criterion, BIC)和p值(Kolmogorov-Smirnov test)最高(表2),表明最合適的劃分是分成4個(gè)子分布。進(jìn)一步對每個(gè)子分布進(jìn)行正態(tài)檢驗(yàn)(Shapiro-Wilks test),結(jié)果表明4個(gè)子分布均服從正態(tài)分布(表3)。

圖3 生態(tài)風(fēng)險(xiǎn)指數(shù)的頻率分布直方圖以及正態(tài)分布擬合曲線Fig.3 The frequency histogram and normal distribution fitting curve of ecological risk index

表2 高斯混合模型訓(xùn)練結(jié)果

表3 子分布正態(tài)性檢驗(yàn)結(jié)果

生態(tài)風(fēng)險(xiǎn)指數(shù)的多峰分布反映了空間樣本的分類,每類樣本具有一定的相似性?；谶@種分類可以在沒有現(xiàn)成標(biāo)準(zhǔn)的情況下,比較客觀地劃分生態(tài)風(fēng)險(xiǎn)等級。本研究以4個(gè)子正態(tài)分布的均值為劃分界限將2015年北京市的生態(tài)風(fēng)險(xiǎn)劃分為五個(gè)等級(表4)。由于本研究提出的生態(tài)風(fēng)險(xiǎn)指數(shù)反映的是區(qū)域的相對生態(tài)風(fēng)險(xiǎn),因此,風(fēng)險(xiǎn)等級也反映了風(fēng)險(xiǎn)的相對高低,當(dāng)研究區(qū)域發(fā)生變化時(shí),風(fēng)險(xiǎn)等級劃分界限也會(huì)相應(yīng)變化。

表4 生態(tài)風(fēng)險(xiǎn)等級劃分以及相應(yīng)風(fēng)險(xiǎn)指數(shù)范圍

3.2 生態(tài)風(fēng)險(xiǎn)的統(tǒng)計(jì)描述及空間分布

2015年北京市生態(tài)風(fēng)險(xiǎn)指數(shù)最小值為0.19,最大值為0.99,平均值為0.64,生態(tài)風(fēng)險(xiǎn)總體上處于低風(fēng)險(xiǎn)接近中等風(fēng)險(xiǎn)水平。具體而言,極低生態(tài)風(fēng)險(xiǎn)區(qū)域面積為869.7 km2,僅占北京市總面積的5.3%,主要分布于北京市北部的部分區(qū)域。 低風(fēng)險(xiǎn)區(qū)域占據(jù)了最大的面積,為8041 km2,占北京市總面積的49%,主要分布于北京的西部和北部；中等風(fēng)險(xiǎn)區(qū)域和高風(fēng)險(xiǎn)區(qū)域面積分別為3889 km2和1640.9 km2,分別占北京市總面積的24.7%和10%,主要分布于北京市的東部和南部的郊區(qū)地帶；極高風(fēng)險(xiǎn)水平區(qū)域面積為1641.1 km2,占北京市總面積的10%,主要分布于北京中部的城市區(qū)域(圖4)。從整體上看,2015年北京市生態(tài)風(fēng)險(xiǎn)呈東高西低和南高北低的空間分布特征,并且具有明顯的空間異質(zhì)性和城郊梯度特征(圖4)。已有其他區(qū)域的研究發(fā)現(xiàn)了相似的現(xiàn)象,如張小飛等[9]和楊勇等[11]的研究分別發(fā)現(xiàn)淮北城市綜合生態(tài)風(fēng)險(xiǎn)和西安長安區(qū)土地利用生態(tài)風(fēng)險(xiǎn)具有很高的空間異質(zhì)性；汪翡翠等[10]的研究顯示京津冀區(qū)域中非城市地區(qū)的土地利用生態(tài)風(fēng)險(xiǎn)低于城市地區(qū)。此外,本研究發(fā)現(xiàn)北京市生態(tài)風(fēng)險(xiǎn)從高風(fēng)險(xiǎn)到低風(fēng)險(xiǎn)的空間變化趨勢與該區(qū)域人工地表從密集到稀疏的空間變化趨勢大體上一致(圖2和圖4),這表明城市化過程驅(qū)動(dòng)的城鎮(zhèn)擴(kuò)張會(huì)導(dǎo)致周邊區(qū)域生態(tài)系統(tǒng)服務(wù)下降,造成生態(tài)風(fēng)險(xiǎn)水平跨級增長,該結(jié)果與汪翡翠等[10]和Kang等[21]的研究結(jié)果相似。這兩項(xiàng)研究表明京津冀的快速城市化地區(qū)較高和高生態(tài)風(fēng)險(xiǎn)等級區(qū)域范圍明顯增加。本研究案例的風(fēng)險(xiǎn)評價(jià)結(jié)果與過去多項(xiàng)研究結(jié)果的相似性表明本研究提出的風(fēng)險(xiǎn)表征方法合理可靠。

圖4 2015年北京市不同生態(tài)風(fēng)險(xiǎn)等級的空間分布 Fig.4 Spatial distribution of different ecological risk levels in Beijing in 2015

3.3 生態(tài)風(fēng)險(xiǎn)指數(shù)與生態(tài)系統(tǒng)服務(wù)價(jià)值間的定量關(guān)系

為分析本研究提出的生態(tài)風(fēng)險(xiǎn)指數(shù)能否準(zhǔn)確地反映區(qū)域生態(tài)系統(tǒng)服務(wù)的損失,以研究區(qū)域全部柵格為樣本,進(jìn)行生態(tài)風(fēng)險(xiǎn)指數(shù)與生態(tài)系統(tǒng)服務(wù)價(jià)值當(dāng)量因子總數(shù)間的回歸分析,結(jié)果顯示生態(tài)風(fēng)險(xiǎn)指數(shù)的對數(shù)值與生態(tài)系統(tǒng)服務(wù)當(dāng)量因子總數(shù)之間呈極顯著線性關(guān)系(P<0.001,R2=0.9715)。因此,本研究提出的風(fēng)險(xiǎn)指數(shù)增加量可以反映生態(tài)系統(tǒng)服務(wù)價(jià)值的損失量。例如,當(dāng)風(fēng)險(xiǎn)指數(shù)從0.4上升到0.5時(shí),相應(yīng)的生態(tài)系統(tǒng)服務(wù)價(jià)值當(dāng)量因子數(shù)損失量為562個(gè)/km2。參考謝高地等[45]的研究,相應(yīng)的生態(tài)系統(tǒng)服務(wù)價(jià)值損失量為191.44萬元/km2。

3.4 不確定性分析和參數(shù)敏感性分析結(jié)果

3.4.1不確定性分析結(jié)果

本研究對生態(tài)風(fēng)險(xiǎn)指數(shù)進(jìn)行了10000次蒙特卡羅模擬,正態(tài)性檢驗(yàn)表明模擬結(jié)果呈正態(tài)分布(Shapiro-Wilks test,P>0.05)。模擬結(jié)果的變異系數(shù)僅為0.079,且3倍的標(biāo)準(zhǔn)差僅為0.081,相對變異指數(shù)(RVI)僅為-5%—5%。幾個(gè)表征不確定性的統(tǒng)計(jì)指標(biāo)值均較小,表明本研究提出的生態(tài)風(fēng)險(xiǎn)指數(shù)的不確定性是可以接受的,能夠較好地表征一個(gè)土地利用類型組成確定的研究區(qū)域的生態(tài)風(fēng)險(xiǎn)(表5和圖6)。

圖5 生態(tài)風(fēng)險(xiǎn)指數(shù)與生態(tài)系統(tǒng)服務(wù)價(jià)值間的定量關(guān)系Fig.5 Quantitative relationship between ecological risk index and ecosystem service value

圖6 蒙特卡羅模擬結(jié)果的累計(jì)概率分布 Fig.6 Cumulative probability distribution of Monte Carlo simulation results

表5 蒙特卡羅模擬結(jié)果的統(tǒng)計(jì)描述

3.4.2參數(shù)敏感性分析結(jié)果

基于蒙特卡羅模擬的結(jié)果,對每個(gè)參數(shù)的敏感性進(jìn)行了分析,結(jié)果如圖7所示。敏感性貢獻(xiàn)度最大的三類參數(shù)分別是水文調(diào)節(jié)服務(wù)當(dāng)量因子、文化娛樂服務(wù)當(dāng)量因子和氣體調(diào)節(jié)服務(wù)當(dāng)量因子,貢獻(xiàn)度分別為17.65%、13.06%和12.76%。9個(gè)權(quán)重的敏感性貢獻(xiàn)度之和為17.33%,相比于所有當(dāng)量因子,權(quán)重的敏感性貢獻(xiàn)度較低。

圖7 不同參數(shù)的敏感性貢獻(xiàn)度比較Fig.7 Comparison of the contributions of sensitivity of different parameters食物權(quán)重:食物供給服務(wù)權(quán)重；原材料權(quán)重:原材料供給服務(wù)權(quán)重；氣體權(quán)重:氣體調(diào)節(jié)服務(wù)權(quán)重；氣候權(quán)重:氣候調(diào)節(jié)服務(wù)權(quán)重；水文權(quán)重:水文調(diào)節(jié)服務(wù)權(quán)重；環(huán)境權(quán)重:環(huán)境凈化服務(wù)權(quán)重；土壤權(quán)重:土壤保持服務(wù)權(quán)重；生物權(quán)重:生物多樣性維持服務(wù)權(quán)重；文化權(quán)重:文化娛樂服務(wù)權(quán)重；食物當(dāng)量:食物供給服務(wù)當(dāng)量因子；原材料當(dāng)量:原材料生產(chǎn)服務(wù)當(dāng)量因子；氣體當(dāng)量:氣體調(diào)節(jié)服務(wù)當(dāng)量因子；氣候當(dāng)量:氣候調(diào)節(jié)服務(wù)當(dāng)量因子；水文當(dāng)量:水文調(diào)節(jié)服務(wù)當(dāng)量因子；環(huán)境當(dāng)量:環(huán)境凈化服務(wù)當(dāng)量因子；土壤當(dāng)量:土壤保持服務(wù)當(dāng)量因子；生物當(dāng)量:生物多樣性維持服務(wù)當(dāng)量因子；文化當(dāng)量:文化娛樂服務(wù)當(dāng)量因子

圖8 五類參數(shù)的敏感性貢獻(xiàn)度比較 Fig.8 Comparison of the contributions of sensitivities of 5 types of parameters服務(wù)權(quán)重:各項(xiàng)生態(tài)系統(tǒng)服務(wù)權(quán)重；供給當(dāng)量:供給服務(wù)當(dāng)量因子；調(diào)節(jié)當(dāng)量:調(diào)節(jié)服務(wù)當(dāng)量因子；支持當(dāng)量:支持服務(wù)當(dāng)量因子；文化當(dāng)量:文化服務(wù)當(dāng)量因子

本研究中使用的生態(tài)系統(tǒng)服務(wù)價(jià)值定量方法為當(dāng)量因子法,該方法所需參數(shù)少,但未考慮小尺度上生態(tài)過程的差異,相比于其他生態(tài)系統(tǒng)服務(wù)價(jià)值核算方法,在應(yīng)用于中小尺度或人工生態(tài)系統(tǒng)中的結(jié)果不確定性更大[46]。本研究的不確定性分析結(jié)果表明即使使用不確定性較大的當(dāng)量因子法,研究區(qū)域的生態(tài)風(fēng)險(xiǎn)評價(jià)結(jié)果的不確定性并不高。因此未來的研究可以進(jìn)一步優(yōu)化生態(tài)系統(tǒng)服務(wù)核算方法以獲取更加準(zhǔn)確可靠的風(fēng)險(xiǎn)評價(jià)結(jié)果。為了增加風(fēng)險(xiǎn)評價(jià)結(jié)果的準(zhǔn)確性,需要考慮評價(jià)空間尺度、區(qū)域社會(huì)經(jīng)濟(jì)差異和自然環(huán)境差異等要素選取合適的生態(tài)系統(tǒng)服務(wù)價(jià)值計(jì)算模型。

4 結(jié)論

以多種類型生態(tài)系統(tǒng)服務(wù)作為城市化區(qū)域生態(tài)風(fēng)險(xiǎn)的評價(jià)終點(diǎn),參考TOPSIS模型構(gòu)建了基于服務(wù)價(jià)值量的生態(tài)風(fēng)險(xiǎn)表征方法,以北京市為例進(jìn)行了方法的應(yīng)用和不確定性分析。研究發(fā)現(xiàn)2015年北京市生態(tài)風(fēng)險(xiǎn)總體上處于低風(fēng)險(xiǎn)接近中等風(fēng)險(xiǎn)水平,生態(tài)風(fēng)險(xiǎn)具有明顯的沿鄉(xiāng)村到城市的梯度逐漸上升的特征,表明北京市城市區(qū)域的擴(kuò)張?jiān)斐芍苓厖^(qū)域生態(tài)系統(tǒng)服務(wù)的下降,導(dǎo)致生態(tài)風(fēng)險(xiǎn)水平的升高。提出的風(fēng)險(xiǎn)指數(shù)用以計(jì)算土地利用結(jié)構(gòu)確定的區(qū)域內(nèi)的生態(tài)風(fēng)險(xiǎn)時(shí),計(jì)算結(jié)果變異較小,表明其具有較高的可靠性。此外,生態(tài)風(fēng)險(xiǎn)指數(shù)同生態(tài)系統(tǒng)服務(wù)當(dāng)量因子間的定量關(guān)系明確,因此可用于快速估算生態(tài)系統(tǒng)服務(wù)價(jià)值損失量。為優(yōu)化本研究提出的生態(tài)風(fēng)險(xiǎn)指數(shù),未來的研究需要選擇可靠性更高的生態(tài)系統(tǒng)服務(wù)價(jià)值量計(jì)算方法,以獲得更加準(zhǔn)確的風(fēng)險(xiǎn)評價(jià)結(jié)果。