大理市森林景观形态与热环境的动态研究

张顺敏,黄晓园,刘俊泽,陈 蓉,李 翔

(西南林业大学 地理与生态旅游学院，昆明 650224)

快速城市化进程导致了热岛效应的发生[1]。根据联合国政府间气候变化专门委员会(IPCC)第六次评估报告显示,2011—2020年,全球地表温度比1850—1900年升高了1.1℃,预计在2021—2040年内达到1.5 ℃[2]。全球气温的升高会导致危害多发并发,其中之一就是城市热环境。城市热环境所体现的城市热岛效应作用于城市,造成了空气污染严重、疾病发病率上升和生物多样性降低等一系列负面影响,热环境正严重影响城市的正常运作和居民的生活质量。在此背景下,可持续发展理念受到世界各国的认同,热环境的研究越发受到全球学者的关注[3]。当前对热环境的研究主要是利用Landsat数据来反演地表温度,如单通道算法[4]、多通道算法[5]、劈窗算法[6]。研究内容主要体现在土地利用类型[7-8]、城市空间格局[9]、植被覆盖度[10]、地理单元划分[11]、景观格局[12]、数值模拟[13]和源汇景观[14-15]等方面,而在森林景观方面多采用景观格局指数来衡量[16]。研究所涉及的指标都是从数值角度对城市的森林空间布局进行量化,反映的是森林整体形态特征,缺乏对森林内部空间格局的考量。运用形态学空间格局分析(Morphological Spatial Pattern Analysis,MSPA)研究森林景观已受到国内外学者的广泛关注[17-19]。目前,MSPA主要运用于森林景观破碎化[20]和构建生态网络[21-23]方面。虽然有部分学者在对热环境进行研究时有所运用,但主要运用在热岛效应的景观连通性方面[24-25],在森林景观形态与热环境的研究上极少运用。

大理市作为中国旅游发展最快的地区之一,是国家园林城市和中国优秀旅游城市,但高强度的人类活动和快速城市化已对区域生态环境产生了巨大的影响,部分地区生态系统功能下降,热环境问题日渐严重,已成为经济社会可持续发展的瓶颈。随着《大理市国土空间规划(2021—2035年)》的推进[26],大理市经济社会发展已迎来新的机会,生态环境必将面临新的挑战。森林景观形态在不同季节对热环境有何差异?区域森林景观与热环境之间存在怎样的耦合关系?是推进生态环境优化亟须解决的迫切问题。针对以上问题,本文基于大理市2006年、2016年遥感数据(http://www.gscloud.cn)和森林资源二类调查数据①开展研究:1)利用2006年、2016年遥感影像进行地温反演,并采用MSPA分析森林景观与热环境的相关性:2)运用普通最小二乘(OLS)和地理加权回归(GWR)建立森林景观形态与热环境的耦合关系模型。通过本文的研究,以期为缓解大理市热环境提供参考依据。

1 研究区概况

大理市位于云南省西部,大理白族自治州中部,地处北纬25°25′—25°58′,东经99°58′—100°27′之间,市辖九镇、一乡、三个街道,面积1 815 km2。市域地形复杂,地势西高东低,平均海拔2 090 m。地貌以山地为主,山地占总面积的70.5%;坝区占总面积的15.8%;洱海水域占总面积的13.7%[27]。属亚热带高原季风气候类型,平均气温15.2 ℃。境内苍山是世界纬度带相同或相近地区植被垂直带谱最完整、植被类型最多样、山地植被保存最自然的山脉之一。大理市森林总面积808.9 km2,森林覆盖率大于46%①,主要植被类型为针叶林、常绿阔叶林、落叶阔叶林等。

2 研究数据和方法

2.1 数据来源

本文所涉及的2006年、2016年遥感影像从地理空间数据云(http://www.gscloud.cn)获取,条代号为131/42,分辨率为30 m(表1)。土地利用数据来源于两期森林资源二类调查数据(1)大理市人民政府.大理市省级森林城市建设总体规划(2017—2026年).2017.。在地温反演之前,对遥感影像进行预处理,包括辐射定标,大气校正,影像裁剪等。所有影像及数据统一至WGS_1984_UTM_Zone_47N坐标系。

表1 研究区遥感影像数据源

2.2 研究方法

2.2.1地表温度反演

本文基于大气校正法,其原理是把大气对卫星传感器的影响从观测到的热辐射总量中去除,进而通过相应公式转换为地表温度(LST)[28]。首先,利用 Landsat热红外波段[TM,ETM(第6波段),TIRS(第10波段)]数据进行热红外波段辐射定标,将像元灰度值转换为辐射亮度值,并通过归一化植被指数(NDVI)以及植被覆盖度(FVC)进行地表比辐射率计算[29];其次,估算出同温下黑体辐射亮度,根据普朗克函数对研究区LST进行反演。参考文献[30],具体计算公式如式(1)—式(3)所示。

Lχ=KGain×VDN+CBias

(1)

(2)

(3)

式中:Lχ为云顶辐射亮度[W/(m2·sr·μm)];

KGain,CBias和VDN分别为转换函数的增益系数、偏移系数和像元灰度值,由影像自带的元数据中可获得;Ly为地表辐射亮度[W/(m2·sr·μm)];Lm,Ln分别为上、下行辐射强度[W/(m2·sr·μm)],t为大气透过率,这3个参数均可由NASA网站(https://atmcorr.gsfc.nasa.gov/)查询得到;e为地表发射率,参考覃志豪等[29]的研究成果计算得出;TS为反演的LST(K);K1,K2均为定标常数,对于TM,K1=607.76 W/(m2·sr·μm),K2=1260.56 K;对于ETM,K1=666.09 W/(m2·sr·μm),K2=1282.71 K;对于TIRS,K1=774.89 W/(m2·sr·μm),K2=1321.08 K。为验证反演结果的准确性,将反演平均LST与同期MODIS数据得到的LST进行对比,结果显示,误差绝对值在2 K以内。可知,反演结果能够满足研究需要,数据可用。

2.2.2景观分类

选取2006年11月17日和2016年11月20日遥感影像,采用监督分类法对土地利用类型进行分类,将其分为林地、耕地、水域、建设用地、未利用地等5类(图1)。为检验分类的准确度,根据每个土地的比例,随机选择100个点,并参考现有土地利用数据和谷歌历史影像,利用混淆矩阵对分类结果进行精度评价。结果显示:两期土地利用类型总体精度分别为92%和90%;Kappa系数为88.04%和85.29%(表2),皆大于80%。

图1 土地利用类型

表2 土地利用精度评价

在土地利用类型基础上,依据《云南植被》[31]和我国现行的《土地利用现状分类标准》[32]中的分类系统,结合两期森林资源二类调查数据(2)大理市人民政府.大理市省级森林城市建设总体规划(2017—2026年).2017.对林地植被进行分类,可将植被类型分为针叶林类和阔叶林类(常绿阔叶林和落叶阔叶林)(图2、图3)。

图2 植被类型

图3 阔叶林类型

2.2.3MSPA方法分析

MSPA是基于腐蚀、膨胀、开启、闭合等数学形态学原理对栅格图像的空间格局进行度量、识别和分割的一种图像处理方法。基于边缘宽度和森林像素之间的连通度,通过将森林斑块分为核心(Core)、孔隙(Perforation)、孤岛(Islet)、桥接(Bridge)、边缘(Edge)、支线(Branch)、环岛(Loop)等7个景观类型(表3),实现对景观碎片量化描述[33]。

表3 景观类型分类

1)基于大理市两期土地利用类型矢量数据,将其转换为30 m栅格数据。2)以研究区的林地作为前景数据,赋值为2;其它用地类型划分为非林地,作为背景数据,赋值为1。生成30 m空间分辨率的森林/非森林二值栅格影像。3)采用MSPA方法,以边缘宽度为2个像元,且8个相邻像元连通的结构来表征每个森林像素所属景观类型(图4)。

图4 2006年、2016年大理市森林空间形态

2.2.4回归分析

1)多元线性回归分析。线性回归分析是一种预测性的建模技术,它研究的是因变量和自变量之间的关系。通过构建2006年和2016年的多元线性回归模型,并运用OLS方法进行诊断,以此作为与其它模型比较的基准。

2)空间自相关和GWR分析。空间自相关常用的检验指标有莫兰指数(Moran′sI)、邻接比(Geary′sC)、热点分析(Getis)、聚合函数(Join count)等,其中Moran′sI指数因其普适性而应用最为广泛,Moran′sI指数值介于-1～1之间,数值大于0,则说明要素之间存在空间正相关;数值小于0,则为负相关;等于0,为不相关[34]。本文将研究尺度设定为1km×1km,首先,利用ArcGIS 10.8对8期LST平均值进行格网分析;其次,运用Geoda软件对其进行局部空间自相关分析,局部莫兰指数(Local Moran′sI)值分别为0.760,0.766,0.760,0.809,0.746,0.751,0.734,0.759,显著性水平P<0.05,通过显著性检验。表明大理市LST值在所选尺度下均呈很强的空间正相关。但显著的空间自相关违背了标准非参数检验(如OLS中的相关检验)所要求的样本间独立性的基本假设,违反这一假设会导致对相关系数的过高估计,从而增加变量之间没有关系的无效假设的相似度,因此还需采用GWR进行分析。基于赤池信息准则(AIC)值,选用高斯函数为权重函数,以每个像元的LST平均值作为GWR的因变量,各驱动因子为自变量,进行GWR分析。

3 结果分析

3.1 热环境时空变化特征

基于均值—标准差法对LST进行分级,分为低温、次低温、中温、次高温、高温区(图5)。由图5可知:2006年高温区主要集中在上关镇、双廊镇、挖色镇、太邑乡和下关镇,以人类活动区域为热源中心,呈现出明显的多中心格局。2016年,在经济快速发展,企业数量急剧上升的背景下,伴随着建筑面积增加,绿地面积减少,新兴工业区排放的二氧化碳和热量增多,高温区逐渐向海东镇和凤仪镇扩散。低温区主要在苍山地区,该区域海拔高且垂直落差大,大面积分布生态林斑块,林木茂密且人类活动少,较符合城市热岛效应的基本特征。

注:“①”上关镇;“②”双廊镇;“③”挖色镇;“④”太邑乡;“⑤”下关镇;“⑥”海东镇;“⑦”凤仪镇;“⑧”苍山。

大理市作为一个高原山地城市,地形对温度起着重要的作用。苍山地区由于海拔较其它地区高,始终为低温分布区。2006年、2016年,热环境的分布格局存在季节性差异,高温区分布集中于夏季,春季则相对较低。随着城市化的发展,在2016年,高温区和次高温区范围逐渐扩大,这表明随着人类活动的加强,其影响逐渐超过四季本身的气候影响。

3.2 影响机制分析

3.2.1不同植被类型的热环境特征

从不同植被类型占大理市林地总面积比例来看,大理市林地植被类型以针叶林为主(3)大理市人民政府.大理市省级森林城市建设总体规划(2017—2026年).2017.。2006年:1)针叶林类占区域总面积的29.23%,占森林总面积的51.64%;阔叶林类占区域总面积的27.59%,占森林总面积的48.36%。2)在阔叶林中,落叶阔叶林占阔叶林的66.34%,常绿阔叶林占阔叶林的33.66%。2016年:1)针叶林类占区域总面积的28.45%,占森林总面积的51.44%;阔叶林类占区域总面积的26.64%,占森林总面积的48.56%。2)在阔叶林中,落叶阔叶林占阔叶林的68.94%,常绿阔叶林占阔叶林的31.06%。林地面积的减少与城市发展有较大关系,其中最明显的为阔叶林类,其落叶阔叶林所占比例增加2.60%,而常绿阔叶林面积则减少,这与落叶阔叶林类多为经济林有关①。

按春、夏、秋、冬来进行比较,由图6可知:1)2006年,常绿阔叶林与LST的相关系数分别为0.002,-0.21,-0.31,-0.27;2016年的为-0.16,-0.17,-0.23,-0.27。2)2006年,落叶阔叶林与LST的相关系数为-0.09,-0.28,-0.17,-0.21;2016年的为-0.02,-0.32,-0.36,-0.25。3)2006年,针叶林与LST的相关系数为0.05,-0.07,-0.22,-0.21;2016年的为-0.26,-0.40,-0.38,-0.26。2006年、2016年,整体上各植被类型与LST呈负相关关系,且落叶阔叶林与LST在相关性上强于常绿阔叶林和针叶林。2016年与2006年相比,整体上植被类型与LST的相关系数有所提高,这与城市扩张致使森林面积减少有极大关系。由于植物生长具有季节性变化,NDVI与LST的相关性表现为夏季最强,冬季最弱。

注:“*”表示P<0.05;“**”表示P<0.01;“***”表示P<0.001。

3.2.2地形因子的热环境特征

在地形因子上:1)海拔与LST相关性最强,呈显著负相关。其中,秋季最强,两年相关系数皆为-0.72;春季最弱,为-0.43和-0.39。这与太阳辐射和大气透射率有关。大理市处于高原山区,海拔高,秋季多为无云期,提高了地面对太阳辐射的吸收。2)由于太阳高度角原因,坡向对温度的影响主要表现在春季。3)在坡度方面,整体上LST随坡度升高而递减,秋季坡度对LST的影响最强,春季最弱。这种差异表明,春季坡度对LST的影响并不是一个决定性因素。

3.2.3森林景观形态的热环境特征

通过格网法将研究尺度设为1km×1km,每个网格内为气温平均值和景观形态的斑块面积比例,通过计算得出森林景观形态与LST的相关系数。2006年、2016年,核心类景观形态与温度LST的关系最强,呈显著负相关关系,其次为支线类,环岛类相关性最低。2006年,相关性强度为核心>支线>边缘>桥接>孔隙>孤岛>环岛;2016年,为核心>支线>边缘>孔隙>孤岛>桥接>环岛。这表明,林地景观形态面积占比越大,与LST的相关程度越高。但随着城市化的不断发展,2016年孔隙、孤岛类的相关性强度强于桥接类,其森林景观逐渐趋于破碎化。

3.3 LST回归分析

3.3.1OLS分析

以2006年和2016年各月份的LST为因变量,海拔、坡度、坡向、NDVI、植被类型,以及7类森林景观形态为自变量做回归分析。由于每个自变量的量纲不统一,首先,对各因素值进行标准化处理。采用高斯函数进行拟合,运用AIC值确定带宽,施瓦兹准则(SC)值确定合适的滞后期长度。利用ArcGIS10.8进行OLS分析,得到模型诊断结果(表4)。在该模型中,调整后,R2最大值为0.786,最小值为0.445,表示该模型可解释因变量中78.6%和44.5%的变化,整体上各月份的模型拟合效果较好。且对于LST,自变量NDVI、常绿阔叶林、孤岛和支线类与其呈正相关关系,NDVI对LST的贡献度最大,其次为支线类(表5)。

表4 OLS模型诊断结果

表5 2016年OLS模型变量

3.3.2GWR分析

由于DEM和核心类在进行GWR模型分析时,存在共线性问题,因此需剔除该变量。GWR模型诊断结果(表6)与OLS模型相比,GWR模型的拟合优度R2均大于 OLS 模型,且改善明显。表明:GWR模型对LST变化量的解释能力以及拟合效果优于OLS模型,可以有效地解释森林景观形态对于LST的演变机制及空间非平稳性特征。

表6 GWR模型诊断结果

通过GWR模型分析可知:植被类型对LST呈现以负相关为主且较集中的区域分布特征;坡度对LST的影响以负相关为主,呈现多格局的分布特征;坡向对LST的影响以正相关为主,呈现中心向四周延伸的分布特征;NDVI对LST的影响以正相关为主,呈现“中部—东西”格局;在7类景观形态中,桥接类、环岛类、孔隙类对LST的影响以负相关为主,孤岛类、边缘类、支线类对LST的影响以正相关为主,均表现为空间多极化分布格局。

4 讨论

由于大理市处于高原山地地区,地形的分割作用使热环境出现多格局的分布特征,在空间上呈聚集状态。这体现了地形对山地城市的重要性[35-36]。因此,要改善城市热环境,必须充分利用城市的地形特征。LST则随土地利用类型变化而变化,建设用地与热环境呈正相关[7],林地与热环境呈显著负相关[11]。以往对于热环境的研究主要基于夏季表面温度以及单一年份的分析[37],对森林景观与热环境的研究通常采用景观格局指数[38-39]进行相关性分析,而MSPA主要运用在森林破碎化[20]、生态网络构建[21-23]和景观连通性分析方面[24-25]。本文运用MSPA方法,拟合多因素对热环境进行综合影响分析,相比以往利用景观格局指数得出的森林景观类型,该方法更为直观,所表达的生态学意义更为丰富。但由于本文主要研究自然因素对热环境的影响,人文因素方面涉及较少,且LST数据主要是基于遥感影像数据获取,容易受到云层和阴影的影响,从而降低了数据质量。因此,在今后的研究中:首先,可以采用新的高质量遥感数据来减少数据质量对研究结果的影响。其次,在因素选择中可进一步丰富,如增加旅游效应、交通通达度等人文因素的影响,以及具有季节特征的风、湿度等气象数据的研究,使研究结果更加显著和准确。最后,森林景观对缓解热环境呈现出季节性差异,可加强对现有森林景观的保护,如:对裸露土地景观进行植树处理,以减少裸露土地面积;通过合理地配置森林景观,以缓解热环境。

5 结论

1)基于多时相遥感数据,采用MSPA、回归分析、空间分析等相结合的方法,对大理市森林景观与热环境的耦合关系进行了定量和定性分析,结果表明:2016年与2006年相比,大理市LST呈上升趋势,并表现为多极化,即出现低温区和高温区聚集状态。高温区主要分布于研究区北部、西南地区和洱海周围,这些区域多为城乡居住区。2016年和2006年相比,建筑面积大量增加,温度也随之上升。主要增长点位于海东镇和凤仪镇,温度增幅超过2 ℃。而低温区主要分布于苍山地区,该区域植被茂密,森林对热环境的缓解作用明显,极大地降低了地表温度。

2)森林景观形态与热环境表现出明显的季节性差异,不同季节多重因素综合影响的强度表现为秋季>夏季>春季>冬季;LST与核心类、孔隙类呈显著负相关,与支线、边缘、孤岛、桥接、环岛类均呈显著正相关。2006年,相关性强度为核心>支线>边缘>桥接>孔隙>孤岛>环岛;2016年的为核心>支线>边缘>孔隙>孤岛>桥接>环岛。森林景观形态面积占比越大,对温度的降温效应越明显。

3)GWR 模型在1km×1km空间尺度上对热环境演变影响因素的解释优于OLS模型,各时期调整后的R2值均大于OLS的,且最大值可达0.815,表明其能较好地反映森林景观与地表温度间的定量关系。