張玉娟，曲建光，葉猛猛

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

0 引言

流域生態(tài)環(huán)境變化及其引發(fā)的各種負(fù)面影響越來越受到人們關(guān)注[1-2]，改變土地的利用方式是人類開發(fā)利用自然環(huán)境最主要的表現(xiàn)形式，不合理的土地利用方式會威脅區(qū)域的生態(tài)環(huán)境健康[3-4]. 生態(tài)風(fēng)險評價是一種能有效支持生態(tài)系統(tǒng)管理的工具，通過生態(tài)風(fēng)險評價能夠評估環(huán)境污染、 人為脅迫或者自然災(zāi)害等外界干擾對區(qū)域生態(tài)系統(tǒng)及其結(jié)構(gòu)產(chǎn)生負(fù)面效應(yīng)的可能性及危害程度[5-7]. 這有助于找出區(qū)域發(fā)展過程中存在的環(huán)境問題，對區(qū)域土地科學(xué)利用與合理開發(fā)、 人類經(jīng)濟(jì)活動和區(qū)域可持續(xù)發(fā)展具有重要意義[8-9]. 生態(tài)風(fēng)險評價出現(xiàn)由傳統(tǒng)生態(tài)風(fēng)險評價, 到區(qū)域生態(tài)風(fēng)險評價, 再到景觀生態(tài)風(fēng)險評價的發(fā)展趨向. 景觀生態(tài)風(fēng)險評價通過選擇與生態(tài)風(fēng)險相關(guān)的景觀指數(shù)結(jié)合權(quán)重構(gòu)建評價模型，如呂樂婷等[5]采用景觀干擾度和景觀脆弱度結(jié)合景觀類型面積比例構(gòu)建景觀生態(tài)風(fēng)險模型對1985—2015年細(xì)河流域景觀生態(tài)風(fēng)險時空分布特征及空間關(guān)聯(lián)格局進(jìn)行評價；汪翡翠等[10]采用斑塊密度和蔓延度結(jié)合層次分析法構(gòu)建景觀生態(tài)風(fēng)險模型對京津冀城市群土地利用生態(tài)風(fēng)險空間分布特征進(jìn)行分析；胡綿好等[11]通過對景觀類型設(shè)置風(fēng)險等級值結(jié)合景觀類型面積比例對德興市生態(tài)風(fēng)險進(jìn)行了分析. 景觀生態(tài)風(fēng)險評價注重風(fēng)險的時空異質(zhì)性和尺度的效應(yīng)，致力于實(shí)現(xiàn)風(fēng)險的空間表征及其可視化，目前得到了廣泛應(yīng)用.

松花江流域哈爾濱段位于松花江流域中部，地屬東北平原，地域平原分布遼闊、 土壤肥沃，是重要的糧食生產(chǎn)區(qū), 在行政區(qū)劃上屬于 “一帶一路”戰(zhàn)略節(jié)點(diǎn)區(qū)域—哈爾濱市. 因此，這一流域的景觀格局不僅對本區(qū)域生態(tài)安全產(chǎn)生影響，對整個松花江流域生態(tài)格局乃至東北對外開放水平都產(chǎn)生重要的作用. 近20年來，該區(qū)域土地利用方式不斷發(fā)生變化[12]，生態(tài)安全性也隨之發(fā)生變化，目前針對該區(qū)域糧食安全和水環(huán)境評價方面有較多的研究[13-14]，但景觀生態(tài)風(fēng)險未見研究. 鑒于此，本文以松花江流域哈爾濱段為研究對象，利用ENVI5.3、 GS+、 ArcGIS10.3、 Geoda095i等軟件提取研究區(qū)域1995年、 2005年和2015年三期土地利用類型，通過景觀易損度指數(shù)及景觀類型面積比重結(jié)合層次分析權(quán)重賦值法構(gòu)建并計(jì)算景觀生態(tài)風(fēng)險指數(shù)，借助冷熱點(diǎn)、 空間變差函數(shù)法，對松花江流域哈爾濱段景觀生態(tài)風(fēng)險的時空變化特征進(jìn)行研究. 松花江流域哈爾濱段景觀生態(tài)風(fēng)險狀態(tài)研究成果是景觀格局優(yōu)化的基礎(chǔ)，能夠明確景觀格局規(guī)劃的保護(hù)重點(diǎn)，有助于完善松花江流域哈爾濱段的生態(tài)格局，為區(qū)域生態(tài)環(huán)境的進(jìn)一步治理以及土地的合理利用提供科學(xué)依據(jù).

1 數(shù)據(jù)準(zhǔn)備及研究方法

1.1 數(shù)據(jù)來源及處理

本文遙感數(shù)據(jù)均來源于地理空間數(shù)據(jù)云，年份分別為1995年、 2005年和2015年，行帶號分別為116/28、 117/28、 118/28、 117/29、 118/29. 非遙感數(shù)據(jù)主要包括哈爾濱市行政區(qū)劃圖，哈爾濱市部分土地利用現(xiàn)狀數(shù)據(jù)，外業(yè)調(diào)查數(shù)據(jù)，哈爾濱市1996、 2006和2016年鑒統(tǒng)計(jì)中各區(qū)縣人口密度數(shù)據(jù)等. 根據(jù)2017年國家標(biāo)準(zhǔn)《土地利用現(xiàn)狀分類》，結(jié)合本研究區(qū)域的主要土地利用方式，將研究區(qū)劃分為耕地、 林地、 建設(shè)用地、 水域、 草地以及未利用地6種土地利用與覆蓋類型. 利用ENVI 5.3遙感影像處理軟件，對下載獲得的初始影像分別進(jìn)行輻射定標(biāo)、 大氣校正、 裁剪、 鑲嵌等預(yù)處理，根據(jù)外業(yè)調(diào)查數(shù)據(jù)、 土地利用現(xiàn)狀數(shù)據(jù)及Google earth地圖等相關(guān)資料，利用最大似然法進(jìn)行遙感影像監(jiān)督分類，獲得松花江流域哈爾濱段1995年、 2005年和2015年土地利用與覆蓋類型分布圖(見圖1)，三期影像分類精度分別為0.871、 0.862和0.888.

1.2 研究方法

1.2.1 景觀生態(tài)風(fēng)險指數(shù)

為建立土地利用類型與區(qū)域生態(tài)風(fēng)險之間的聯(lián)系，利用景觀易損度結(jié)合各景觀類型面積比重及層次分析權(quán)重賦值法構(gòu)建景觀生態(tài)風(fēng)險指數(shù)(ecological risk index, ERI)，該值的大小反映了區(qū)域景觀生態(tài)風(fēng)險的高低，其計(jì)算模型如下：

(1)

式中：IER為生態(tài)風(fēng)險指數(shù)；n為景觀類型的數(shù)量；Ai為研究區(qū)內(nèi)景觀類型i的面積(km2)；A為研究區(qū)域總面積(km2)；Wi為研究區(qū)內(nèi)景觀類型i所反映的生態(tài)風(fēng)險強(qiáng)度權(quán)重；Ri為景觀類型i的易損度，選擇與景觀類型易損程度相關(guān)的景觀破碎度、 分維數(shù)和分離度組合構(gòu)建，按下式進(jìn)行計(jì)算.

Ri=aCi+bNi+cFi

(2)

式中：Ci為景觀類型i的破碎度；Ni為景觀類型i的分維數(shù)；Fi景觀類型i的分離度[15]，為了消除不同量綱的影響，對計(jì)算所得的各景觀指數(shù)進(jìn)行歸一化處理；a、b、c分別為Ci、Ni和Fi的權(quán)重.

評價指標(biāo)權(quán)重賦值是否合理將直接影響景觀生態(tài)風(fēng)險指數(shù)的合理性，本文采用AHP[16-17]的原理和方法確定各景觀類型生態(tài)風(fēng)險強(qiáng)度權(quán)重和景觀破碎度、 分維數(shù)、 分離度指數(shù)權(quán)重. 首先查閱現(xiàn)有的景觀類型生態(tài)風(fēng)險強(qiáng)度參數(shù)相關(guān)文獻(xiàn)中各種景觀類型權(quán)重值[18-20]，將這些權(quán)重按平均值大小進(jìn)行排序，作為初始專家經(jīng)驗(yàn)構(gòu)造判斷矩陣，通過對該判斷矩陣進(jìn)行一致性檢驗(yàn)來判斷所構(gòu)造的矩陣是否具有滿意的一致性. 本文通過計(jì)算獲得一致性比例的計(jì)算結(jié)果為0.007 8(<0.1)，表明判斷矩陣具有滿意的一致性，獲得6種土地利用類型生態(tài)風(fēng)險的AHP權(quán)重分別為：建設(shè)用地0.366、 林地0.049、 水域0.107、 草地0.075、 耕地0.146、 未利用地0.257. 同理，獲得景觀破碎度、 分維數(shù)、 分離度指數(shù)權(quán)重值分別為0.539、 0.297、 0.164，一致性比例的計(jì)算結(jié)果為0.003 5(<0.1).

1.2.2 景觀生態(tài)風(fēng)險空間統(tǒng)計(jì)方法

1) 空間變異函數(shù). 在地統(tǒng)計(jì)學(xué)上，通過半方差函數(shù)的理論方法分析地理變量空間上的關(guān)系及分布格局[21]，實(shí)現(xiàn)判斷地理變量在單元網(wǎng)格尺度劃分下是否具有空間相關(guān)性，從而獲得地理變量在多尺度網(wǎng)格劃分中，哪個尺度更為合理. 半方差函數(shù)的理論模型計(jì)算公式為：

(3)

式中：γ(h)為變異函數(shù)值；Z(xi)和Z(xi+h) 是在空間單元xi和xi+h上的景觀生態(tài)風(fēng)險值(i= 1, 2, 3, …,N(h))；N(h)是分割距離h的樣本量.

2) 冷熱點(diǎn)分析. 莫蘭指數(shù)被用于判定地理變量空間自相關(guān)關(guān)系，包括全局自相關(guān)莫蘭指數(shù)和局部自相關(guān)莫蘭指數(shù). 其中全局自相關(guān)莫蘭指數(shù)用于分析地理變量空間相似性聚集度，局部自相關(guān)莫蘭指數(shù)用來識別地理變量在不同空間單元的高值簇或熱點(diǎn)區(qū)與低值簇或冷點(diǎn)區(qū)的空間分布[23]. 全局自相關(guān)莫蘭指數(shù)和局部自相關(guān)莫蘭指數(shù)計(jì)算公式為：

(4)

(5)

式中：Pe和Pf分別為地理單元e和地理單元f位置的函數(shù)值(e≠f)；Wef為權(quán)重系數(shù)矩陣；n為地理單元數(shù). 莫蘭指數(shù)的判斷區(qū)間為[-1, 1]，區(qū)間(0, 1]表示正自相關(guān)，值越大，其正自相關(guān)性越高；區(qū)間[-1, 0)表示負(fù)自相關(guān)，值越小，其負(fù)自相關(guān)性越高；值等于0時表明變量值在空間上是相互獨(dú)立分布的.

2 實(shí)證研究

2.1 研究區(qū)域概況

松花江流域哈爾濱段總面積5.31萬km2，區(qū)域位于松嫩地塊-松嫩斷陷東南隆起區(qū)，東南臨張廣才嶺支脈丘陵，北部為小興安嶺山區(qū)，中部有松花江通過，山勢不高，河流縱橫，平原遼闊. 區(qū)域?qū)俦睖貛Ъ撅L(fēng)氣候區(qū)，大陸性氣候特點(diǎn)明顯，多年平均降水量一般在500 mm左右，汛期6-9月的降水量占全年的60%～80%，冬季12月至次年2月的降水量僅為全年的5%左右. 研究區(qū)域內(nèi)的大小河流均屬于松花江水系和牡丹江水系，松花江干流由西向東貫穿研究區(qū)域中部，是研究區(qū)域灌溉量最大的河道. 研究區(qū)域水資源特點(diǎn)是自產(chǎn)水偏少，過境水較豐富，時空分布不均勻，表現(xiàn)為東富西貧.

2.2 景觀生態(tài)風(fēng)險小區(qū)劃分及生態(tài)風(fēng)險指數(shù)計(jì)算

景觀生態(tài)風(fēng)險指數(shù)由景觀格局指數(shù)構(gòu)建，因此具有尺度效應(yīng)，本文分別按照邊長4和6 km正方形單元將研究區(qū)域分割成2 300和1 065個研究小單元，利用ArcGIS 10.3中的漁網(wǎng)工具和分割工具對三期土地利用與覆蓋圖進(jìn)行單元分割，對獲得矢量小單元進(jìn)行批量轉(zhuǎn)換生成為柵格小單元. 利用景觀指數(shù)提取軟件Fragstats 4.2批量提取三期每個網(wǎng)格4個類型水平的景觀指數(shù)：斑塊面積(TA)、 斑塊數(shù)(NP)、 分維數(shù)(FRAC_MN)和分離度指數(shù)(DIVISION). 在Excel中按照式(1)編輯函數(shù)，計(jì)算三期每個小單元的景觀生態(tài)風(fēng)險值，并將該值作為網(wǎng)格中心點(diǎn)的屬性值.

2.3 空間變異函數(shù)理論模型擬合分析

不同步長下的景觀生態(tài)風(fēng)險變異函數(shù)模型能準(zhǔn)確反映景觀生態(tài)風(fēng)險在不同尺度上的變化特征. 確定最佳的采樣尺度是進(jìn)行景觀生態(tài)風(fēng)險空間分異及關(guān)聯(lián)研究的基礎(chǔ)與前提. 本文采用的是規(guī)則格網(wǎng)取樣，可將格網(wǎng)間距4和6 km作為步長大小來分析在這兩種采樣尺度下的空間變異性，從而確定能在一定空間范圍內(nèi)反映景觀生態(tài)風(fēng)險的最佳尺度. 利用地統(tǒng)計(jì)學(xué)GS+軟件進(jìn)行兩種尺度下景觀生態(tài)風(fēng)險的統(tǒng)計(jì)分析，構(gòu)造景觀生態(tài)風(fēng)險空間變異函數(shù)理論模型，完成景觀生態(tài)風(fēng)險數(shù)據(jù)變異函數(shù)模型擬合(見表1和表2).

表1 4 km尺度變異函數(shù)擬合模型參數(shù)

注：Rc為結(jié)構(gòu)方差占總方差的比值C/(C+C0)，下同.

表2 6 km尺度變異函數(shù)擬合模型參數(shù)

續(xù)表2

4 km尺度下，三期景觀生態(tài)風(fēng)險均為線性模型擬合效果最好，三期決定系數(shù)R2分別為0.974、 0.974和0.969，此時殘差RSS均最小. 結(jié)構(gòu)方差占總方差的比值Rc呈現(xiàn)逐漸增加趨勢，表明在4 km尺度劃分下，1995—2015年，研究區(qū)域景觀生態(tài)風(fēng)險受結(jié)構(gòu)性因子的影響有所增加[24-25]. 研究區(qū)三期景觀生態(tài)風(fēng)險空間分異的有效變程A為18.31 km(大于4 km)，因此研究區(qū)域景觀生態(tài)風(fēng)險在4 km尺度劃分下有較高度的空間相關(guān)性.

6 km尺度下，三期景觀生態(tài)風(fēng)險均為指數(shù)模型擬合效果最好，三期決定系數(shù)R2分別為0.985、 0.981和0.990，此時殘差RSS均最小. 結(jié)構(gòu)方差占總方差的比值Rc呈現(xiàn)逐漸增加趨勢，表明在6 km尺度劃分下，1995—2015年，研究區(qū)域景觀生態(tài)風(fēng)險受結(jié)構(gòu)性因子的影響也是有所增加的. 研究區(qū)三期景觀生態(tài)風(fēng)險空間分異的有效變程A為37.08 km(大于6 km)，因此研究區(qū)域景觀生態(tài)風(fēng)險在6 km尺度劃分下具有非常高度的空間相關(guān)性.

4 km和6 km兩種尺度劃分都適合于研究區(qū)域景觀生態(tài)風(fēng)險研究，二者相比較，6 km尺度下決定系數(shù)R2值更大、 殘差RSS值更小、 結(jié)構(gòu)方差占總方差的比值Rc更大些，景觀生態(tài)風(fēng)險空間相關(guān)性更高些，因此，選擇6 km網(wǎng)格劃分尺度進(jìn)行研究區(qū)景觀生態(tài)風(fēng)險評價. 6 km尺度下，區(qū)域景觀生態(tài)風(fēng)險網(wǎng)格分布圖(見圖2).

2.4 生態(tài)風(fēng)險時空變化分析

根據(jù)自然斷點(diǎn)法將景觀生態(tài)風(fēng)險指數(shù)劃分為5個景觀生態(tài)風(fēng)險等級，即低風(fēng)險級、 較低風(fēng)險級、 中等風(fēng)險級、 較高風(fēng)險級、 高風(fēng)險級. 采用ArcGIS 10.3地統(tǒng)計(jì)模塊的克里格插值法對劃分等級后的景觀生態(tài)風(fēng)險離散點(diǎn)進(jìn)行插值，形成三期景觀生態(tài)風(fēng)險等級連續(xù)分布圖(見圖3)，分別對三期景觀生態(tài)風(fēng)險等級進(jìn)行統(tǒng)計(jì)(見表3).

表3 三期景觀生態(tài)風(fēng)險分級統(tǒng)計(jì)

近20年間，各風(fēng)險等級分布面積不斷變化，但分布的主要區(qū)域位置無明顯變化. 高風(fēng)險區(qū)主要分布在主城區(qū)，較高風(fēng)險區(qū)主要分布在高風(fēng)險區(qū)周圍以及各縣域城鎮(zhèn)中心，原因?yàn)橹鞒菂^(qū)分布較大面積的建設(shè)用地景觀和極小面積的林地景觀，縣域城鎮(zhèn)中心主要以建設(shè)用地景觀和耕地景觀為主. 低風(fēng)險區(qū)主要分布在通河縣、 木蘭縣、 方正縣和賓縣境內(nèi)，原因?yàn)檫@些縣域分布著很大面積的林地景觀. 中等風(fēng)險區(qū)面積最大，占研究區(qū)域41%～47.6%，主要分布于大面積的耕地景觀范圍內(nèi). 較低風(fēng)險區(qū)主要分布在中等風(fēng)險區(qū)與低風(fēng)險區(qū)的過渡區(qū)域范圍. 1995—2015年，高風(fēng)險區(qū)和較高風(fēng)險區(qū)面積呈現(xiàn)先增加后減少的變化，高風(fēng)險區(qū)分布面積比例由6.60%到11.00%再到7.90%，較高風(fēng)險區(qū)，分布面積比例由13.80%到19.90%再到19.30%. 低風(fēng)險區(qū)、 較低風(fēng)險區(qū)和中等風(fēng)險區(qū)面積呈現(xiàn)先減少后增加的變化，低風(fēng)險區(qū)分布面積比例由14.50%到14.00%再到14.50%，較低風(fēng)險區(qū)分布面積比例由17.50%到14.10%再到16.10%，中等風(fēng)險區(qū)分布面積比例由47.60%到41.00%再到42.20%. 從整體水平上看，研究區(qū)景觀生態(tài)風(fēng)險呈現(xiàn)先上升后下降的變化.

2.5 區(qū)域生態(tài)風(fēng)險等級的空間關(guān)聯(lián)格局分析

表4 全局莫蘭指數(shù)及驗(yàn)證值

利用ArcGIS 10.3空間統(tǒng)計(jì)工具分析研究區(qū)域景觀生態(tài)風(fēng)險空間關(guān)聯(lián)性. 利用全局莫蘭指數(shù)計(jì)算工具計(jì)算景觀生態(tài)風(fēng)險全局莫蘭指數(shù)值和Z值(見表4). 三期全局莫蘭指數(shù)分別為0.512 5、 0.555 2和0.605 5，均大于0；Z值分別為24.25、 25.69和27.44，說明研究區(qū)域景觀生態(tài)風(fēng)險等級存在很強(qiáng)的正相關(guān)性，呈聚類模式的空間分布.

通過對景觀生態(tài)風(fēng)險等級分布進(jìn)行熱點(diǎn)分析，獲得相鄰區(qū)域之間景觀生態(tài)風(fēng)險的空間關(guān)聯(lián)模式. 利用ArcGIS 10.3中聚類和異常值分析工具進(jìn)行相似景觀生態(tài)風(fēng)險的空間聚類分析，結(jié)果見圖4. 從圖4分析可知，1995—2015年，研究區(qū)域景觀生態(tài)風(fēng)險高-高聚集、 低-低聚集分布區(qū)域并未發(fā)生較大變化，但逐漸表現(xiàn)為較分散的分布. 高-低、 低-高聚集出現(xiàn)的非常少，只在1995年和2015年聚集分布于研究區(qū)域內(nèi).

3 結(jié)語

研究區(qū)域景觀生態(tài)風(fēng)險具有高度的空間相關(guān)性，結(jié)構(gòu)性因素是三期景觀生態(tài)風(fēng)險空間相關(guān)性的主要影響因素，且結(jié)構(gòu)性因素對景觀生態(tài)風(fēng)險的影響程度逐漸加深.

從整體上看，1995—2015年松花江流域哈爾濱段景觀生態(tài)風(fēng)險先升高后下降. 源于1995—2005年間，區(qū)域忽視生態(tài)環(huán)境保護(hù)，轉(zhuǎn)林為耕、 開荒建城，使得生態(tài)用地不斷轉(zhuǎn)換為非生態(tài)用地，導(dǎo)致區(qū)域景觀生態(tài)風(fēng)險增加. 2005—2015年間，隨著全球生態(tài)環(huán)境保護(hù)意識增強(qiáng)，區(qū)域制定了多項(xiàng)自然保護(hù)區(qū)發(fā)展戰(zhàn)略，逐步恢復(fù)松花江沿岸生態(tài)環(huán)境，退耕還林、 退耕還草，使得部分非生態(tài)用地逐漸轉(zhuǎn)換為生態(tài)用地，區(qū)域景觀生態(tài)風(fēng)險降低.

1995—2015年，研究區(qū)域生態(tài)風(fēng)險高-高聚集、 低-低聚集分布區(qū)域并未發(fā)生較大變化，但逐漸表現(xiàn)為較分散的分布. 高-低、 低-高聚集出現(xiàn)的非常少，只在1995年和2015年聚集分布于研究區(qū)域內(nèi).

本文通過景觀易損度指數(shù)及景觀類型面積比重結(jié)合層次分析權(quán)重賦值法構(gòu)建景觀生態(tài)風(fēng)險指數(shù)，具有一定的創(chuàng)新性. 生態(tài)系統(tǒng)服務(wù)是鏈接人類福祉和生態(tài)系統(tǒng)的橋梁，將景觀生態(tài)服務(wù)價值與景觀格局相結(jié)合進(jìn)行生態(tài)風(fēng)險研究將是課題進(jìn)一步研究內(nèi)容.