本文开展的城市热岛时间尺度特征分离,主要面向逐月尺度进行。因此,其前提在于获取研究区的逐月地表温度时间序列。然而,MODIS提供的逐月地表温度产品MOD/MYD11C3空间分辨率较低（0.05°）,无法满足区域尺度上城市热岛分析的需求,同时,MODIS地表温度产品合成周期越长,引起白天、夏季、不同地表覆盖类型下地表城市热岛的偏差越大^{[ 17]}。在另一方面,Terra和Aqua卫星在研究区每日过境时间存在变化,时间纠正难度很大,但该变化具有的周期相似性,使我们可以对具有1 km空间分辨率的MOD/MYD11A1产品进行月合成。然而,我们的前期研究发现,该产品在冬季对于城市区域具有大量的缺失像元^{[ 19]},而MOD/MYD11C3则不存在该问题。因此,采用以下方法构建空间分辨率为1 km的逐月地表温度时间序列。

首先,基于MOD/MYD11A1产品统计各月各像元晴空观测次数。当月内晴空观测天数超过3天时,对MOD/MYD11A1产品提供的地表温度进行平均,作为当月的地表温度月合成值。在该过程中发现,每年12月至次年3月,城市区域有较多缺失像元。

其次,根据尺度转换原理^{[ 24]},对具有较多缺失像元的月份,在Kustas等提出的地表温度采样模型^{[ 25]}基础上进行改进,进而对MOD/MYD11C3地表温度进行空间分辨率提升。除常用的归一化差值植被指数NDVI外,还从地表辐射与能量平衡的角度,选用地表反射率、地表反照率等作为地表温度空间变异的描述参量：

(1)

式中, T _s6km为6 km分辨率上地表温度估计值; NDVI_6km表示对1 km分辨率16天合成NDVI产品进行月合成并重采样至6 km分辨率的植被指数; ρ_6km表示7波段6 km分辨率反射率向量, ρ_{6km, i}表示对MODSIS第 i波段1 km分辨率反射率产品做筛选、月合成,及重采样后,所得6 km分辨率反射率数据; α_6km表示对MODIS 1 km反照率数据重采样后得到的6 km分辨率反照率月合成数据; ξ_6km表示6 km分辨率上误差补偿项; a,b,c,d分别为各项对应系数。

其中,在假定地表为朗伯体的前提下计算地表反照率^{[ 26]}：

(2)

式中, ρ _i（ i=1, 2, …, 7）为MODIS第 i波段的地表反射率。

根据尺度推演的原理,将在6 km尺度上建立的函数 f应用于1 km尺度并进行误差补偿,从而获得1 km尺度上地表温度的估计值：

(3)

式中,下标1km表示参数的空间分辨率为1 km;其余参数的物理意义与式(1)类似。

通过上述方法,实现MOD/MYD11A1和MOD/MYD11C3地表温度产品的集成,构建了研究区的地表温度时间序列。其中,对于地方太阳时10:30的地表温度数据,其时间跨度为2001—2012年,对于地方太阳时13:30则为2003—2012年。

由于温度随海拔高度具有垂直递减的效应,故城市热岛量化需保证城市与作为参考的乡村在同一海拔高度水平上。因此,本文将研究区限定在高程小于100 m的范围。城市选取以空间连续性为标准,提取自研究区地表覆盖类型数据。乡村作为城市热岛研究的基准,选取标准为：与城市隔离、地表覆盖类型以植被为主、温度较为稳定,受城乡距离影响较小。因此,以城市为中心,间隔5 km做缓冲带,选取平均温度趋于稳定的缓冲带作为乡村范围,并确定剔除水体和城镇像元的45~50km缓冲带作为乡村。2001—2012年共12年间乡村的分布情况显示,乡村缓冲带有逐年向外延伸的特征,且在2001—2007年间较为明显,自2007年后,乡村缓冲带位置趋于稳定。为消除乡村位置变化给城市热岛效应研究带来的不确定性,选取2007年乡村缓冲带代表乡村（图1）。

图1

Fig. 1

	Figure Option ViewDownloadNew Window
	图1 各年城乡区域示意图Fig. 1 Urban and rural area in each year.

根据各年城市与乡村提取方案,在乡村平均温度基准上,以城乡温度差作为城市热岛强度的度量标准。

首先,计算逐像元城乡温度差（ Ts _d）,用于衡量城乡地表温度的基本差异：

(4)

式中, Ts _{u_pixel}表示城市中每个像元地表温度值, Ts _r表示乡村平均温度。

其次,以各年提取出的城市与乡村区域平均温度差和城市中发生热岛区域的城市热岛强度作为衡量城市热岛的2项指标,分别记为UHII₁和UHII₂。其中,把城市中 T s _d大于乡村平均温度标准差（ σ _RST）的区域记为发生热岛的区域。

时间序列模型是依据变量自身的变化规律,利用外推机制描述时间序列的变化。当时间序列在年度尺度上具备近线性平稳趋势、在年内尺度上具备明显的季节循环特征时,即可表述成4个不同的成分：季节成分 S _t（ t指某一特定时间点）、趋势项成分 T _t、循环项成分 C _t与随机误差或不规则变动成分 I _t。其中, S _t表示以季节为固定周期的循环特征, C _t是时间序列表现出的周期长度不固定且长于季节周期的循环特性, T _t表示时间序列的长期发展方向。由于测算趋势在实际工作中有一定难度,因此通常趋势项 T _t与循环项 C _t被归并为循环—趋势项 TC _t,不再进行区分^{[ 27]}。

4个成分的函数关系通常可通过加法模型（Additive Model）或乘法模型（Multiplicative Model）进行描述,形式如下^{[ 28]}：

(5)

(6)

其中,式(5)为加法模型,式(6)为乘法模型; X _t表示 t时刻的观测值。

如何选取2种模型,取决于时间序列自身的特点及各组分之间的内在关系。加法模型适用于各组分相互独立的情形,同时季节变动幅度基本保持不变。乘法模型适用于各组分相关的情形,同时季节变动的幅度随趋势上升而增大。

Makridakis^{[ 22]}等给出了经典季节分解方法,基本思路是通过对时间序列进行12项简单移动平均后,对原序列与移动平均后序列作差（乘法模型对原序列与移动平均后序列作商）,在差值（商值）序列中取各年相应月份平均值,即为各月季节项。3×3项移动平均用于减少偶然趋势对序列的影响,以获取趋势—循环项（T-C）。原始序列去除季节项及趋势—循环项后,即得到月度高频项。

加拿大统计局Dagum等^{[ 23]}开发了X-11-ARIMA季节分离模型,它是用于提取给定时间序列组分的季节调整方法。X-11-ARIMA模型的特点是在季节调整之前,通过建立ARIMA模型,对给定时间序列做前向及后向预测,从而补充数据,以保证季节调整过程的数据完整性。具体方法如下：

通过中心化12项移动平均进行趋势—循环项初始估计。移动平均在消除12阶不变季节性时,复合了简单12项移动平均及用于中心化的简单2项移动平均,能够在保留线性趋势的同时减小不规则项的方差。公式如下：

(7)

式中, 表示 t时刻趋势—循环项, Y _t±i（ i=0,1, … ,6）表示 t± i时刻的原始时间序列值。

为保留线性趋势,3×3移动平均用于初步估计季节项,同时,为了使连续12个月内季节因子之和接近0,还需对初步估计的季节项进行标准化。公式如下：

(8)

(9)

式中, 表示 t时刻季节项, S+I表示原始序列 Y与趋势—循环项 T_C⁽¹⁾作差得到的季节及不规则组分, 表示 t时刻标准化后的季节项。

对原始序列做初步季节调整,通过对调整后序列做13项Henderson移动平均^{[ 29, 30]},再次估计趋势—循环项。

(10)

(11)

(12)

式中, A⁽¹⁾为初步季节调整后序列, T_C⁽²⁾为二次估计趋势—循环项。

最终季节项估计通过3×5移动平均给出,标准化后的季节项用于与原始序列求差,估计最终季节调整后序列,公式如下：

(13)

(14)

(15)

式中, S⁽²⁾表示最终季节项, S ^s(2)表示标准化后最终季节项, A⁽²⁾为最终季节调整后序列。

最后,2 H+1项Henderson移动平均用于从最终季节调整后序列中提取出最终趋势项及不规则项。

(16)

(17)

式中, T_C⁽³⁾为最终趋势项, I⁽³⁾为最终不规则项。

地表温度时间序列是量化城市热岛效应的基础。故在使用前文建立的地表温度时间序列前,需验证其有效性和精度。由于逐月地表温度数据采用两种方法集成得到,故分别对12月至次年3月、整个时间序列分别进行了检验。

首先,评价地表温度空间分辨率提升方法。选取每年12月至次年3月这四个月份的地表温度月合成数据,从中提取对提升方法具有验证意义的研究区内城乡有效像元,与对应区域的地表温度空间分辨率提升结果做对比。在验证之前重新筛选回归方程,剔除系数大于100的第五波段反射率及反照率,剔除由Terra卫星第五波段反射率存在条带引起的10:30地表温度数据条带问题。同时剔除显著性大于0.01的回归因子,以提升少量13:30地表温度数据回归效果不佳的问题。验证样本量及验证结果如表2所示。由于本文以城乡平均温差作为城市热岛效应的基本衡量指标,因此,城、乡区域各自的温度误差范围在5K以内是可接受的。

表2

Table 2

表2(Table 2)

表2 提升结果与月合成有效值的误差分析 Table 2 Error analysis between upscaling results and valid data of monthly-integrated MODIS land surface temperature data 地方太阳时 样本量 乡村区域 城市区域 R2 RMSE R2 RMSE 10:30 48 0.9886 0.71 0.9642 1.64 13:30 40 0.9911 0.58 0.9311 1.70

表2 提升结果与月合成有效值的误差分析 Table 2 Error analysis between upscaling results and valid data of monthly-integrated MODIS land surface temperature data

第二,对完整时间序列进行验证时,由于乡村平均温度作为在城市热岛强度计算的基准,需要验证其可靠性。由于乡村地表覆盖类型多为植被,遥感数据得到的地表温度实为冠层温度,理论上,应与气象站获取的离地1.5 m处空气温度较为接近^{[ 31]}。此外,由于每天正午日照最强,午后2小时左右通常是一天中气温最高点,该时刻与Aqua卫星的过境时间较为接近。因此,一天中的最高气温应与Aqua MODIS提供的白天乡村地表温度具有良好的趋势性,并与Terra MODIS提供的白天乡村地表温度表现出相似的变化规律。因此,研究选取中国地面气候资料月值数据集中北京站的月值平均最高气温序列作为参考,分析其与地表温度时间序列中的乡村地表平均温度与之间的相关性,13:30和10:30样本量分别为120和144,结果如图2所示。

图2

Fig. 2

	Figure Option ViewDownloadNew Window
	图2 乡村地表温度与北京站月最高气温对比Fig. 2 Comparison between RST time series andmothly maximum air temperature gotfrom in situ meteorological datasets

北京市平均城乡温差（UHII₁）上午10:30对应的12年时间序列及13:30对应的10年时间序列如图3a所示。从UHII₁时间序列可以看出,北京市主城区在一年中的多数月份出现了热岛效应,每年8月热岛强度最大,12月至次年3月出现冷岛效应。以北京市10年间13:30热岛强度为例,每年1月和8月平均UHII₁分别为-0.022和7.226 K,标准差分别为1.317和0.524。热岛发生区域城市热岛强度（UHII₂）上午10:30对应的12年时间序列及13:30对应的10年时间序列如图3b所示,相较UHII₁序列,UHII₂序列整体趋势与其较为接近,冬季、春季波动较大。

图3

Fig. 3

Figure Option ViewDownloadNew Window

图3 2种量化方法对应的城市热岛强度时间序列 (a) 10:30(2003—2012年) 和13:30(2001—2012年)城乡平均温差时间序列,(b)10:30(2003—2012年) 和13:30(2001—2012年)热岛发生区域城市热岛强度时间序列Fig. 3 Time series of UHI intensities calculated based on two methods (a) Mean temperature difference between urban and rural at 10:30(year2003-2012) and 13:30 (year 2001-2012), (b) Urban heat island intensity where urban heat island happened at 10:30(year2003-2012) and 13:30 (2001-2012)

年内周期性变化表现出明显的季节特征,且在年度上呈近线性平稳趋势。该特征与前文所述时间序列加法模型特征相吻合,适合通过时间序列建模,进行结构性组分分离,挖掘序列各组分的变化规律,从而揭示北京市城市热岛的深层特征。

分别采用经典季节分离模型和X-11-ARIMA模型对UHII₁时间序列进行成分分离,得到其季节项（图4）。由于建模方法不同,经典季节分离模型中各年份对应季节的季节项是以恒定量出现,而X-11-ARIMA模型通过季节异常值调整,得出季节因子的整体趋势与经典季节分离模型相似的同时,还表现出更符合客观规律的年度季节差异。从图4可以看出,季节调整主要集中在冬季。这可以很好地反映两点客观事实：a）X-11-ARIMA模型对每年同一月份进行分析,由于季节成分的年变化不应出现陡变,由调整结果可见,这部分调整对于冬季地表温度空间分辨率提升产生的误差有修正作用;b）由于地表温度及各地表产量数据集的反演误差、研究区域云污染等主要误差源存在,X-11-ARIMA模型的季节调整能够抑制误差对季节的影响,起到去噪作用。UHII₂时间序列进行成分分离结果与UHII₁相似,此处不做单独说明。

图4

Fig. 4

	Figure Option ViewDownloadNew Window
	图4 2种季节分离方法从城乡平均温差时间序列中提取季节项 (a)13:30(2003—2012年),(b)10:30 (2001—2012年)Fig. 4 Seasonal factors of mean temperature difference between urban and rural got from two seasonal decomposition methods (a)13:30(2003-2012), (b) 10:30 (2001-2012)

时间序列的趋势—循环项综合了时间序列的长期演变方向及其持续的、周期性的波动,往往比季节周期长,且每个周期的长度可以不同。分别采用经典季节分离模型和X-11-ARIMA模型对剔除季节因子的UHII₁时间序列进行进一步成分分离,得到趋势—循环项如图5(a)~(b)所示。对比经典季节模型与X-11-ARIMA模型分离出的趋势—循环项可以看出,经典季节模型分离出的趋势—循环项仍具有较强的波动性,X-11-ARIMA模型的趋势—循环项相对平滑,跨度较长,且不在数据所处的周期范围内重复,更符合时间序列趋势—循环项的客观规律。此外,针对发生城市热岛的区域,以地方太阳时13:30数据为例,对UHII₂时间序列做X-11-ARIMA成分分离,如图5c所示,分离出的趋势—循环项时间序列具有很强的平稳性,说明12年间北京市日间出现城市热岛的区域,热岛强度稳定。

图5

Fig. 5

Figure Option ViewDownloadNew Window

图5 趋势—循环项时间序列 (a) 2种季节分离方法从10:30城乡平均温差时间序列中提取(2001—2012年),(b) 2种季节分离方法从13:30城乡平均温差时间序列中提取(2003—2012年),(c)X-11-ARIMA模型提取出13:30热岛发生区域城市热岛强度时间序列趋势—循环项及不规则项(2003—2012年)Fig. 5 Trend Cycle factors (a) get from mean temperature difference between urban and rural at 10:30 based on two seasonal decomposition methods (2001-2012), (b) get from mean temperature difference between urban and rural at 13:30based on two seasonal decomposition methods (2003-2012), (c) get from urban heat island intensity time series where urban heat island happened at 13:30 based on X-11-ARIMA methods (2001-2012)

综合分析2种模型分离出的季节项与趋势—循环项,由于二者在表现UHII₁的季节特征时采用的方法不同,从而对序列趋势—循环项的影响较为显著。城市热岛效应是多种因素共同作用的产物,这些因素随着年份不同也会发生变化。使用各年同一月份季节因子取平均值得到的季节成分,仅仅消除了重要的季节峰值,虽然能在整体趋势上符合客观规律,却会造成趋势—循环项波动难以解释。因此,X-11-ARIMA模型更加适用于分析城市热岛效应的时间尺度特征。

X-11-ARIMA模型分离出的UHII₁季节成分时间序列,通过3×3移动平均结合季节调整,减小了原始UHII₁序列的波动性,放大了其季节特征。如图4所示,北京市热岛强度的季节性最高值和最低值分别出现在每年8月和1月。每年4~9月,北京市季节性热岛强度为正,其中4~8月热岛处于上升期,8~9月热岛强度出现快速下降趋势。同时,每年10月至次年3月,北京市季节性热岛强度为负,其中10月至次年1月热岛强度季节性特征整体呈下降趋势,而1~3月呈上升趋势。另外,11月季节性热岛强度相较10月出现1K左右的提升。

UHII₁城市热岛表现出的季节变化特征与土地利用状况、季节性地表覆盖、地物热特性,及气候因子等有密切联系^{[ 20, 32, 33]}。①由于北京市夏季高温持久稳定,昼夜温差小,城市区域非渗透性地表使得城市区域地表面热容量比郊区大,夏季白天城市地表增温率明显高于乡村地区。同时,夏季农田植被覆盖度较高且土壤含水量丰富,潜热蒸发作用较强,大量土壤和植被冠层中的热量通过潜热蒸发释放,也在一定程度上扩大了城乡地表温度差异,使得城市热岛强度较大。②北京冬季寒冷干燥,降水少,阳光充足,平均云量最少,这些因素提升了城市地表热释放效率。与此同时,秋季农作物收割,乡村地表由大面积植被覆盖变为裸露,造成9~10月城市热岛出现陡降。北方冬春季节植被覆盖度低,乡村区域地表裸露,土壤干燥,潜热蒸发减少,升温速率加大,使得城乡温差减小,热岛强度减弱或出现冷岛。③11月热岛强度有一小幅提升,及次年5月后热岛强度小幅下降,与北方冬小麦播种情况较为吻合。冬小麦通常在9月底或10月初播种,因此11月乡村地表植被覆盖度较10月高,热岛强度上升。次年5月,冬小麦成熟,收割后出现热岛强度下降。

X-11-ARIMA模型分离出的UHII₁趋势—循环项时间序列如图5(a)~(b)显示,2005—2009年间,城市热岛处于较高值,与快速城市扩张联系密切。2009年后,城市选区面积趋于稳定,由于面积较大,以城乡平均温差衡量的热岛强度趋势相对减弱。模型分离出北京城市热岛UHII₁的趋势—循环时间序列,拉伸了UHII₁序列中未能突出表达的年际特征,揭示了以城乡平均温差为衡量指标的城市热岛长期演变规律。此外,针对发生城市热岛的区域,以地方太阳时13:30数据为例,对UHII₂时间序列做X-11-ARIMA成分分离,如图5c所示,分离出的趋势—循环项时间序列具有很强的平稳性,说明12年间北京市日间出现城市热岛的区域,热岛强度稳定。

X-11-ARIMA模型同时还可分离出不规则变化成分,如图5c所示。周纪在研究最大夜间热岛强度时发现,相对湿度和降水等因素也会对城市热岛产生影响^{[ 19]}。在城市地表热岛的研究中,除了考虑地表特征,还需考虑气候、气象因子的驱动作用^{[ 17]}。因此,时间序列建模分析方法从原始序列中提取出的不规则变化成分,为定量研究气候、气象等偶然因素对城市热岛的影响提供了可能。