城市化过程是不连续的、具有梯度性的，其必然会导致区域地表覆被格局的改变^[1]，土地覆被时空演化在LUCC^{[2, 3, 4]}、城镇扩展^{[5, 6, 7, 8]}、区域景观格局^{[9, 10, 11, 12]}及城市热环境^{[13, 14]}、水环境^[15]等方面具有广泛应用。在一定范围内生态系统也随着地表覆被的改变而改变，这种改变在生态不敏感地区引起的效果甚微^[16]，但是在生态脆弱地区，比如石羊河流域，人类活动的干扰会导致整个流域生态系统的紊乱，引发一系列的环境为题^[17]。为了更好的指导城市建设和流域生态保护，找出近十多年的地表覆被变化规律，选择流域地表覆被空间格局演化作为研究生态环境效应的典型代表，将流域土地覆被的时间与空间研究范式结合起来，并以此立题，作为本研究的出发点。

流域生态环境问题是目前城市化与流域可持续发展普遍关注的问题^[18]，本文以Ridd于1995年提出的植被-不透水面-土壤3种基本要素构成的城市地表覆被系统^[19]为基础，结合流域1998—2010年NDVI、夜间灯光数据和土壤含量，运用改进的计算人居环境指数(Human_Settlements_Index，HSI)算法^[20]计算土地人工覆被指数(LCI)。然后对其进行时间和空间尺度上的研究，利用小波分析、空间变差函数和间隙度指数等方法构建测度模型，分析了石羊河流域地区土地覆被空间演化的过程和格局。最后尝试利用回归拟合方程找出影响地区覆被分异和变化的驱动机制，以此反映近十多年来石羊河流域地表演化过程。

石羊河流域位于河西走廊东段，祁连山北麓，介于东经101°22′—104°16′、北纬36°29′—39°27′之间。以高山、中高山、低山丘陵、沙漠、冲积平原构成全流域地貌形态，地势南高北低，自西南向东北倾斜。面积约41600km²，流域行政区划包括武威市凉州区、古浪县、民勤县全部及天祝县乌鞘岭以北部分地区，金昌市永昌县和金川区全部以及张掖市肃南县部分地区，共3市7县(区)。流域生态环境脆弱、沙漠化严重，十多年来“西部大开发”快速的城镇化建设使这些问题更加突出。

图 1 研究区概况 Fig. 1 Sketch map of study area

图选项

本研究1998年、2002年、2006年、2010年流域NDVI、夜间灯光数据和土壤信息分别来源于SPOTVGT数据、美国国家地球物理数据中心(NGDC)的DMSP-OLS和中国科学院计算机网络信息中心国际科学数据服务平台的Landsat TM影像，流域基础数据主要来源于国家自然科学基金委员会“中国西部环境与生态科学数据中心”(http://westdc.westgis.ac.cn)。

统一把流域基础数据集的坐标系转换为西安1980坐标系，由于转换后的NDVI与DMSP-OLS图像均为1km×1km空间分辨率，所以无需进行重采样，用流域边界裁切图像。在ENVI中将4个年份7—8月SPOT Vegetation的6期数据转换为真实NDVI值，并提取最大值；再提取OLS最大值和最小值，根据公式对其进行归一化处理：

点击浏览公式

式中，OLS_nor为归一化灯光指数，OLS是灯光指数实际值，OLS_max和OLS_min分别为灯光指数最大值和最小值。结合Lu等提出的提取人居环境指数(HSI)空间信息算法^[20]：

点击浏览公式

式中，NDVI_max为7—8月份NDVI最大值，解释为若一个像元中居住区所占比重大，那么导致NDVI_max较低、DMSP-OLS较高，所以HSI指数偏高。改进的计算土地覆被指数(LCI)，以土地人工覆被综合表征，修改后公式如下：

点击浏览公式

式中，SOIL为土壤比例指数，为解译的4期16景TM影像(均在7月或8月份)土地利用类型估算土壤比例，然后采用极差标准化方法，对其进行标准化处理。同式(2)，若一个像元中人居环境、即不透水地表所占比重大，那么导致SOIL较低、HSI较高，所以LCI指数偏高。

小波变换是时间(空间)频率的局部化分析，它通过伸缩平移运算对信号(函数)逐步进行多尺度细化，最终达到高频处时间细分，低频处频率细分，能自动适应时频信号分析的要求，从而可聚焦到信号的任意细节。反映景观要素在不同尺度上的特征及其相互关系，揭示空间格局的多尺度和等级结构特征，计算公式为：

点击浏览公式

式中，W(a,x_j)为景观要素在空间尺度为a时的波值，x_j为窗函数的中值，g(x_i-x_j)为窗函数，对地形因子以及土地覆被指数的空间采样序列进行小波方差和小波一致性分析。通过误差分析和相关分析检验模型的模拟效果，统计参数包括均方根偏差(RMSE)、平均偏差(MBE)、相对偏差(RE)^[21]。

空间变差函数也叫半变异函数，是地统计分析的特有函数。区域化变量Z(x)在点x和x+h出的值Z(x)和Z(x+h)差得方差的一半称为区域化变量Z(x)的半变异函数，记为r(h)，2r(h)称为变异函数。根据计算有：

点击浏览公式

半变异函数曲线图反映了一个采样点与其相邻采样点的空间关系(图 2)。理论上，当采样点间的距离为0时，半变异函数值应为0；但由于存在测量误差和空间变异，使得两采样点非常接近时，它们的半变异函数值不为0，即存在块金值(C₀)。当采样点间的距离h增大时，半变异函数r(h)从初始的块金值达到一个相对稳定的常数时，该常数称为基台值(C₀+C)。当半变异函数的取值由初始的块金值达到基台值时，采样点的间隔距离称为变程(a)。在变程范围内，样点间的距离越小，其相似性，即空间相关性越大。当h>R时，区域化变量Z(x)的空间相关性不存在^[22]。

图 2 变差函数理论方差图 Fig. 2 Model variogram on theory

图选项

间隙度是Mandelbrot于1983年在《自然界分形几何》一书中最早提出，被广泛地应用于景观质地分析和景观的分形特征分析中，间隙度指数(Lacunarity Index)用于分析地表景观空间形态演化的过程。间隙度分析方法主要包括移动窗口法和查分盒子法，前者是对土地覆被总体特征及分布的一种度量，后者则用于识别土地覆被空间分布的局部异质性^[23]。

全局间隙度指数计算公式为：

点击浏览公式

式中，M为边长为r的计数窗口中土地覆被栅格变化的数量，Q(M,r)为计数窗口中栅格变化的比例。

局部间隙度指数计算公式为：

点击浏览公式

式中，M_w=∑_i,jn_r(i,j)，n_r(i,j)=v-μ+1，M_w为局部区域(w×w)内边长为r滑动窗口中土地覆被变化的数量，n_r(i,j)是局部区域i和j的列相对高度，v为具有最大象元值的的立方体盒子，μ为具有最小象元值的立方体盒子^[23]。

本文综合以上方法，首先从直线样带的基础上，用小波分析找出影响流域土地覆被演化的因素，然后在面状范围内用变差函数分析LCI在时间点上的变化格局分布，最后从时空动态角度用间隙度指数研究LCI在时间段上的变化过程。

由于所选数据年份的限制，以每隔4a为1个时间点，利用流域3种因子合成1998年、2002年、2006年、2010年土地人工覆被指数(LCI)，综合表征土地覆被。用县区级行政区划裁切生成的栅格，得出流域各县区LCI(表 1)。植被部分主要是利用NDVI体现的，由于采用1km×1km分辨率，干旱内陆河流域城市规模较小，2010年武威市建成区为16.75km²，金昌市建成区为55.51km²，占流域比重分别为0.04%、0.13%，其中城市植被更小，所以本文没有识别，植被部分仅指大区域耕地、林地、草地等。

以土地覆被指数LCI显示石羊河流域地表盖度的空间格局，如表 1，图 3所示。民勤LCI均值在流域中最高，这是由于LCI的计算是以NDVI为基数，其被巴丹吉林和腾格里沙漠包围，NDVI较低；天柱地处祁连山区，LCI均值最低。1998—2010年间流域LCI呈现先增后减的趋势，1998—2002年均值从0.3991增长到0.4282，平均年增长率为0.0073，增幅较大；2002—2010年的两个阶段均值从0.4178持续减小到0.4098，平均每年减少0.54%。除金川外各县区LCI演变过程和流域基本相似，由于金川为工矿型城市，建设集中，在最近十几年间建设用地持续扩张，所以其均值不断增加。从其余6个县区的均值变化中看出，在1998—2002年间增幅最大的为天柱，年均增长2.81%；增幅最小的为肃南，年均增长0.45%。2002—2010年间减幅最大的是天柱，年均减少1.96%；减幅最小的是凉州，年均减少0.41%。

图 3 1998—2010年流域土地覆被空间分布图 Fig. 3 Land-cover in Shiyang River Basin(1998—2010)

图选项

为表现出地表覆被组分在空间演变上的地区差异及幅度，采用剖面分析方法进行组分的空间分异分析。布设2条沿主要交通和水域演变方向的剖面线(图 1)，比较LCI沿2个方向的梯度变化。图 4中N-S、E-W两个方向为典型样带，囊括了所有地表类型，以1km×1km为单元，分析流域LCI在相同位置上不同时间内的变化特征。

图 4 4个时期不同方向上的LCI空间分异坡面图 Fig. 4 The profiles of fractions of LCI in different directions(1998—2010)

图选项

可以看出样带1(N-S方向)的变化不如样带2(E-W方向)明显，样带2大体沿G30连霍高速横贯整个流域，西部处于金川和永昌的城市交界处、中部位于凉州城市建设中心区，由于西部大开发的政策，城市增长迅速，所以LCI变化巨大，并且随着时间的推移明显呈递增趋势，其中凉州城市值中心最高，2010年为1.1349，2002—2006年增幅最大达209%；东部位于古浪和天柱直抵乌鞘岭，地表盖度较低，由于国家实行退耕还林草等，使得2002—2010年间LCI呈递减趋势，说明地区地表植被覆盖度越来越高，环境越来越好。样带1沿石羊大河走势纵贯流域，北部处于两大沙漠的分界线上，植被覆盖度低，所以LCI较高，并且在时间上几乎没有变化；中部位于211省道沿线凉州至民勤段，由于城市建设的原因，1998—2010年间民勤及周边LCI逐渐增大，形成的几个低谷是民勤和凉州耕地；南部位于肃南和天柱边界上的水源涵养区，LCI时间上变化幅度较小，空间上起伏较大，值较高的地区说明有冰川和永久性积雪，较低说明是植被分布，12年间LCI总体呈下降趋势，可以看出山区植被增多，流域荒漠化程度改善。

为找出影响土地覆被空间分布的因素，以流域2条样带的DEM和2010年LCI为例，进行小波一致性分析(图 5)。可以看出，小尺度(32km以下)和中尺度(32—64km)上地形变化对LCI的作用很小，只有零散的圆圈、箭头分布；当达到阀值(64km)以后，地形对LCI的作用逐渐明显，存在空间位置上的一致性，将样带推演到全流域，反映出大尺度地形因子是形成LCI格局分异的主要控制因子，即石羊河流域由南向北“山区-平原-荒漠”大范围上的依次分布，是形成LCI空间格局的一大原因。

图 5 不同方向上LCI与DEM的小波一致性分析 Fig. 5 Multi-scale wavelet coherence analysis between LCI and DEM in different directions

图选项

定义采样步长为5km，将LCI赋予重采样的格网中心点计算实验变差函数，得出石羊河流域4a的方差图及变差拟合曲线(表 2，图 6)。可以看出，通过最小二乘法拟合的模型在前3a相同，都是指数模型，2010年为球状模型，决定系数都达到0.9以上，说明拟合精度很好。

图 6 流域土地覆被梯度变差演化 Fig. 6 Variogram and its evolution of land-cover

图选项

变程(a)的不断增大，特别是在2010年激增为22.5km，说明流域地表覆被在区域内部相互影响，即城市建设和退耕还林等人为影响因素较高。将基台值(C₀+C)定义为系统方差，基台值和块金值(C₀)的差值为结构方差，块金值不断变大，说明表征LCI变化主导因素的空间相关性在增大。结构方差与系统方差的比值即块金系数(C/[C₀+C])，在4个典型年份中不断减小，在采样误差等同的情况下，说明象元内部变异性减小、外部差异变大，即LCI区域间变化增大，由于人为因素对流域地表覆被空间差异的影响越来越大，相反自然因素的相对重要性有所减弱。

LCI的异向方差3D拟合模型显示出变差函数曲线在各个方向上都具有规律性，反映了流域LCI随距离变化的空间自相关性，总体上拟合模型呈现四周高中间低的“凹型”漏斗结构，并且随着时间的推移，这种结构越来越明显，说明石羊河流域LCI高低值在区域上发生了集聚，而且高值区的值越来越高，低值区的值越来越低。在0°(东—西)和90°(南—北)方向上、45°(石羊河沿线)和135°(312国道武威、金昌段)方向上的拟合效果均较好，变差函数曲线的变化趋势基本相同，决定系数都大于0.9，也说明了LCI在区域间各个方向上的空间自相关作用较好。

计算流域LCI 1998—2002年、2002—2006年、2006—2010年的变化量，提取大于0的组分，利用差分盒子法分别测算3个阶段LCI变化量的局部间隙度指数(图 7)。1998—2002年局部间隙度指数高值区分布零散，变化量在全流域没有形成集聚簇；低值区分布在古浪、金川中心城区的周边地区，说明这部分地区LCI空间变化量很小，区域协同性较强。2002—2006年间隙度指数变化值较小，高值区主要分布在祁连山沿线，分散化程度减弱，说明这段时间这片区域内部分地区的地表覆被变化较平凡，并且LCI高低分化明显，这主要是由于地形复杂和土地利用类型较多，寒漠和林草交错分布；另外在金川和凉州城区形成低值集聚区。2006—2010年在民勤的外围两大沙漠和地区祁连山区形成了高值区，表明沙漠环境也在变化；在流域内5大县(区)城市建设中心区形成了低值集聚区，说明建成区面积不断扩大。综合来看，1998—2010年间流域城镇空间扩展速度加剧，同时部分地区地表覆被也在发生变化。

图 7 不同时期流域LCI变化的局部间隙度分析 Fig. 7 Local lacunarity analysis of LCI change in different stages of evolution(1998—2010)

图选项

运用移动窗口法测算3个时间段内LCI变化量的全局间隙度指数(滑窗尺度变化范围为2—39象元)，生成双自然对数曲线(图 8)。可以看出，2002—2006年、2006—2010年ln(全局间隙度指数)在初始阶段都小于1998—2002年，但是随着ln(滑窗尺度)的增大，逐渐超过了1998—2002年，并且所用滑窗尺度越来越小，说明局部间隙度指数高值和低值零散分布区在减少，向着集聚区发展，即小区域内部性的LCI变化区逐渐消失，取而代之出现了一些大范围的LCI变化区，区域异质化程度增强。

图 8 不同时期全局间隙度指数曲线 Fig. 8 The curves of global lacunarity index in different stages of evolution(1998—2010)

图选项

从上文分析可知，影响流域地表覆被变化的是人为和自然双重因素。为了找出空间上导致流域地表覆被变化的因素，人为影响因素的数据主要来自于1998—2010年石羊河流域各县区和周边县区的人口和经济数据，来源于“甘肃、青海、内蒙古统计年鉴”；自然影响因素主要来自于1998—2010年石羊河流域及周边气象站点7、8月份的气温和降水数据(包括山丹、永昌、武威、民勤、马鬃山、松山、景泰、阿拉善右旗、门源、巴音毛道)，来源于中国气象科学数据共享服务网。

对驱动因子进行Kriging插值，即数据空间化，再根据上文公式计算这12a的LCI，与同期人为(人口、GDP)、自然(气温、降水)因素分别进行回归分析，公式如下：

点击浏览公式

式中，R_LCI表示相关系数，x_i为其中一个影响因子， 为因子4a平均值，n=4，导入ArcGIS计算，结果如图 9所示(通过P=0.001的显著性水平检验)。对整个流域的4组相关系数进行统计，其中相关系数平均值最大为0.83(图 9)，57.4%的区域相关性在0.8以上；最小为0.66(图 9)，26.8%的区域相关性在0.8以上，表明研究区LCI与自然和人为因素的相关性都较好。

图 9 1998—2010年LCI与自然因素相关系数分布图 Fig. 9 Linear correlations of LCI with natural factors for the period 1998—2010

图选项

从自然因素与流域LCI的相关性分布图(图 9)可以看出，相关性较高的区域主要分布在祁连山水源涵养区，相关性平均值最大，为0.74(气温)和0.76(降水)，这部分地区为林地和草地；另外，高值区还分布在北部荒漠与绿洲交互区，表明沙漠和戈壁演化的主要驱动力来自于自然因素，综合说明了在人迹罕至的高山和荒漠地区，自然条件是影响流域地表覆被的主要因素。

图 10 1998—2010年LCI与人为因素相关系数分布图 Fig. 10 Linear correlations of LCI with human factors for the period 1998—2010

图选项

本文以NDVI、DMSP-OLS和Landsat TM为基础数据，研究了石羊河流域地表覆被的时空演化分布格局，发现整个流域LCI随时间的改变而不断演化。1998—2010年间，LCI均值经历了由小变大、再由大变小的过程，其中2002年达到峰值，之后均值逐渐减小。各县区LCI的演化规律也各不相同，凉州和金川的城市中心区由于城市不断扩张，所以LCI一直在增大，肃南、天祝高植被覆盖区地表盖度逐渐变好，民勤、古浪等荒漠区LCI变化不明显。接着利用典型样带和小波分析，印证了这一情况。然后利用实验变差函数和间隙度指数模拟了流域地表覆被正向演化过程，揭示了区域性聚集演化规律逐渐增强，零散性演化分布逐渐减少，即LCI高值和低值分布区逐渐形成集聚簇，而零散分布区逐渐减少。

最后根据流域气象和人口经济数据找出驱使LCI演化的规律分布，可以看出，流域城镇集中分布的区域及其周边地区LCI变化主要受人口和经济因素影响，在远离城镇人烟稀少的地区，LCI变化受气温和降水(即自然因素)影响较强，研究结果为石羊河流域生态保护建设提供了参考。

地表覆被时空变化是区域城镇空间演化的主要因素，研究城镇演化和地表演化的方法由一定的相通性，所以本文将研究城镇空间扩展的方法引入到土地覆被变化过程中，由于所选年份比较接近，土地覆被变化不是很明显；并且流域及周边的统计数据有限，驱动力分析不够充分全面。下一步将利用长时间序列，结合高精度遥感影像，选择流域典型城市进一步研究。