自从分形的概念由Mandelbrot[1]提出,很快便引起许多学科领域的关注,成为描述复杂、不规则几何特征的主要工具,广泛应用于各个领域研究中。自20世纪80年代初,Batty将分形理论应用于城市形态研究以来[2],国内外展开了大量的基于分形理论的城市形态研究。其中理论分析主要包括:Batty和Longley研究了城市边界和城市土地利用形态分形[3],陈彦光等研究了城市形态熵与分形维数[4,5],Fotheringham、Batty、Longley、White、Engelen、陈彦光、况颐等在城市形态分形模拟方面开展了大量工作[6-10]。实证研究主要包括:Frankhauser测算了包括中国北京和台北在内的世界上许多城市形态的分维值[11],Benguigui等人将特拉维夫整个城市区域的分形维数作为时间的函数研究城市动态演化的发展趋势[12]。Feng等用网格法计算了杭州城市形态的分维特征,揭示了杭州市1949-1996年间城市形态和土地利用结构的演化特征[13]。姜世国等通过计算容量维数,半径维数和信息维数研究了北京市的城市分形形态[14]。

以上城市分形理论研究及典型应用,为深入理解城市分形特性和城市空间的复杂性提供了有价值的参考。但是,由于城市建筑高度分布信息难以获取,在过去的研究中,城市形态分形的测度方法都是基于二维平面。然而在城市的演化中,物质空间形态不仅是用地面积上的扩展,同时建筑高度也不断发生着变化,呈现出一种动态的随机的变化过程[15]。而现代城市建设为提高城市空间利用和空间配置效率,城市形态呈现出越来越强的三维特征,城市三维有机增长成为城市空间发展的重要目标,需要新的基于城市三维空间结构的定量测度方法,来准确全面地探寻城市空间结构演变规律和时空特征。因此,本文将传统的基于二维计盒法的城市空间形态分维测度方法扩展到三维,详细介绍了三维计盒法的计算原理与实现过程;并且针对城市分形研究中人工判定法不能够精确识别无尺度区的问题,借鉴了二阶导数的方法进行分形无尺度区的自动识别;最后使用2003年及2012年扬州市中心城区数据作为案例区,对三维计盒维数方法进行了实证研究。

计盒法的基本思想是利用不同尺度的标尺量测空间。对于特定尺度r,假设量测空间所需的标尺数量为N(r) 。那么,在一个分形体维数估计过程中,测量体 (如量规、盒子) 的线性尺度与测量体的数量之间服从负幂律关系:

式中:D即为分形维数。在一维情况下,标尺为线段。在二维情况下,标尺为正方形或者矩形。在三维情况下,标尺为正方体或者长方体。计盒法的“盒子”指的就是线段、矩形 (正方形) 和立方体 (长方体) 等标尺,覆盖整个空间所需的标尺数量即为“非空盒子数”。在二维情况下,不同矩形标尺必须具有相同的长宽比,以统一尺度确定的标准。理论上,在长宽比恒定的条件下,以长作为尺度和以宽作为尺度的维数计算结果一样。类似的,在三维情况下,不同立方体标尺必须具有相同的长宽高比,如果长宽比不同则无法统一尺度标准,长宽高都可以作为尺度的衡量标准。对(1)式两边取对数得:

式中:C为常数。由此可知,在分形情况下,尺度r与非空盒子数N(r)在双对数坐标系上呈直线关系。那么,可利用最小二乘法对对数变换后的尺度r与非空盒子数N(r)进行线性回归,回归直线的斜率-D的绝对值就是分形维数。

Benguigui等曾对二维计盒法测度城市分形特征的原理进行了详细的阐述[12,13,16]。三维计盒法是二维计盒法向高维度的扩展,算法基本思想与二维算法类似,只是城市空间结构不需要经过任何的信息压缩而被直接投影到三维空间,而且用于分割被测对象的不是小正方形,而是小长方体。也即,采用不同尺度大小的长方体对三维空间进行量测,计算该尺度下覆盖城市建筑空间所需要的长方体数目(N(r)),计算结果用于建立标度r与非空盒子数N(r)的回归曲线。算法的详细描述如下:假设三维城市空间为M,可以用长为L、宽为W、高为H的长方体完全覆盖。首先以L×W×H长方体为标尺,取其尺度r₁= L。那么,量测整个城市空间M所需的立方体数目,也即非空盒子数目必然为N(r₁)= 1。接着将长方体长宽高各二等分,形成大小为L/2×W/2×H/2的标尺,其尺度为r₂= L/2 ,覆盖整个城市空间所需的标尺数目为N(r₂) 。进一步利用大小为L/4×W/4×H/4,尺度r₃= L/4的长方体标尺量测城市三维空间,所需标尺数目为N(r3) 。同样,逐步按比例减小长方体标尺的大小,到第n步的长方体标尺大小为 ,并取尺度 ,覆盖M空间所需非空长方体数目记为N(r_n) 。如果城市三维空间是分形的,则尺度r与其对应的非空盒子数N(r)满足(1)式的幂次关系,分维值可根据(2)计算得出。

以扬州市中心城区为研究案例区域,其范围包括扬溧高速公路以东、廖家沟以西、启扬高速公路以南和仪扬河—吴州路以北的区域,如图1所示。2012年该区域总人口为76.92万,面积为220 km,城市建设用地面积为89.26 km。

计算城市三维分维需要的数据主要有:城市建筑物轮廓,建筑物高度及建筑物属性数据。使用的基础数据包括2012年扬州市中心城区0.5 m分辨率的Geo Eye影像,与200年扬州市中心城区0.61 m分辨率的Quik Bird影像,《扬州市城市总体规划 (20092020) —规划区土地利用现状图》以及《扬州市城市总体规划 (2002-2020) —规划区土地利用现状图》CAD图。数据获取流程如图2所示,首先通过e Cognition软件对两期高分辨率遥感影像自动获取建筑物轮廓数据[17],并对结果进行人工检验和校正;然后通过2003年土地利用现状图,2009年土地利用现状图 (对近两年数据进行人工补充),为建筑物轮廓数据进行属性赋值;最后由“E都市—扬州”三维模型 (http://yangzhou.edushicom/) 获取建筑物楼层数,通过实际调查数据获取不同属性建筑物的楼层平均高度,在Arc GIS中对提取的建筑物轮廓进行高度赋值。由于“E都市—扬州”没有历史数据,所以2003年建筑高度数据通过建筑物阴影及建筑物高度之间的关系估算建筑物的高度[18,19]。并选择两年建筑物高度无变化区域,抽样调查建筑物实际层数进行误差验正,绝大多数的误差能控制在6%以内,结果是完全可以接受的。城市的土地利用可按功能划分成不同的职能类型,参照《城市用地分类与规划 建设用地 标准GB50137-2011》中城市用地分类,将建筑物类型划分为住宅建筑,商业与公共建筑,工业建筑。2003年与2012年不同类型建筑三维模型如图3所示。

图1 扬州中心城区研究范围Fig. 1 The study area in Yangzhou city

图2 扬州市三维数据获取流程图Fig. 2 Yangzhou 3D morphology data acquisition

扬州市中心城区三维城市空间可用16820 m×16820 m×84 m的长方体完全覆盖,其中16820 m×16820 m为扬州市中心城区的外接矩形,84 m为建筑物最高高度。首先,将扬州市中心城区建筑物高度矢量数据进行栅格化。为便于计算,栅格分辨率取0.5133 m,以使栅格总数为32768×32768,也即2×2。这样,长方体标尺的平面维度可从32768×32768依次对半递减,直到1个栅格为止;而其垂直维度则从84 m依次对半递减直到84/2m为止,如表1和表2所示。理想情况下,为最大程度反映被测空间的性质,最小长方体标尺越小越好。但是,栅格数据大小随栅格分辨率的减小成倍增长,太小的栅格分辨率造成数据处理的困难。因此,长方体标尺最小取到0.5133 m×0.5133 m×(84/2) m为止。实验结果表明,0.5133 m的栅格大小能够将扬州市中心城区分形特征的无尺度区包含进来。2003年及2012年不同类型建筑尺度r及其对应的非空盒子数值统计如表1和表2所示。

图3 扬州市2003年与2012年不同类型建筑三维模型Fig. 3 The 3D model of different types of buildings in Yangzhou

对于自然界与人类社会中大量存在的随机分形,不像数学上的分形,具有在无穷尺度上的自相似或自仿射性,其分形特征只是在一定尺度范围出现,这个尺度范围称之为无尺度区[20]。无尺度区是研究随机分形的基础,准确求解无尺度区,不仅有助于抑制各种分形量的发散,也是分形特征值准确、可靠的基础[21]。城市地理系统的分形性质当然也只表现在一定尺度范围之内,而关于无尺度区的判定,在城市分形研究中大部分采用人工判定法或者经验判定,即根据lnr - ln N(r) 曲线,凭人的经验找出其中线性关系较好的一段作为无尺度区[13,14,22,23]。这种判别方法主观性比较强,往往不能够进行精确的识别和计算,并且对于同一组数据,不同观察者得到的无尺度区的也会有一定的差异,导致最后得到的分形维数可能有较大的差别。

为了提高无尺度区识别的客观性与准确性,采用基于lnr - ln N(r) 曲线二阶导数信息的自动识别方法[24]。第i个点 的一阶导数,也即局部斜率,可由下式计算:

式中: ln′N(r) 为第i个点的局部斜率;K为lnr - ln N(r) 曲线上点的数量。理想分形情况下, lnr - ln N(r) 曲线应为一条直线,即所有 (ln r - ln N(r)) 点的局部斜率保持一样。实际情况中, lnr - ln N(r) 曲线仅在某一区域近似为一条直线,这个区域就是无尺度区。在无尺度区内, (ln r - ln N(r)) 点的局部斜率在某一固定值上下微幅波动。因此,若对局部斜率再次求导,即求取lnr - ln N(r) 曲线的二阶导数。那么,无尺度区内点的二阶导数,也即一阶导数 (局部斜率) 的变化率,都应在0附近微幅波动。第i个点 (ln r_i- ln N(r_i)) 的二阶导数可由下式计算:

通过以上分析可知,无尺度区的寻找可转换为寻找lnr - ln″N(r) 曲线中取值在0附近的一段连续区域。并且,应寻找尽量长的连续区域,以使无尺度区包含更多的点。具体搜索方案为:在lnr - ln″N(r) 曲线中选择M(3≤M≤K - 5) 个连续的点,假设起始点为第j(3≤j≤K - M - 1) 个点,那么这M个点对应的坐标则为 (ln r_i,ln N(r_i)) (i = j,j + 1,...,j + M - 1) 。假设这M个点合围x轴的面积为S_(M,j)。那么,在点数M相同的情况下, S(M,j)越小就代表点列在0附近的波动幅度越小。因此,无尺度区的寻找就转化为寻找尽量小的S_(M,j)和尽量大的M。S_(M,j)可由下式计算:

表1 2003年扬州市不同类型建筑各尺度下非空盒子数N(r)、盒子边长与高度数据Tab. 1 Data of N(r), length and height of the box of different types of buildings in Yangzhou in 2003

表1 2003年扬州市不同类型建筑各尺度下非空盒子数N(r)、盒子边长与高度数据Tab. 1 Data of N(r), length and height of the box of different types of buildings in Yangzhou in 2003

表2 2012年扬州市不同类型建筑各尺度下非空盒子数N(r)、盒子边长与高度数据Tab. 2 Data of N(r), length and height of the box of different types of buildings in Yangzhou in 2012

表2 2012年扬州市不同类型建筑各尺度下非空盒子数N(r)、盒子边长与高度数据Tab. 2 Data of N(r), length and height of the box of different types of buildings in Yangzhou in 2012

对于lnr - ln″N(r) 曲线中的K - 4个点,可得到2 + 3 + ... + K - 6 =(K - 4)×(K - 8)/2个这样的连续点序列,即计算这 (K - 4)×(K - 8)/2个点列对应的所有S_(M,j)值。实际选择过程中,需要设定一个S_(M,j)的阈值μ,满足S_(M,j)≤μ的点列的斜率波动幅度被认为满足无尺度区要求,成为无尺度区的候选点。选取μ=0.5,即在μ≤0.5的情况下,选取点列数最多的连续区域作为无尺度区。

依据式 (3)、(4)、(5) 获得城市三维分形无尺度区后,分别对2003年及2012年无尺度区的点列使用最小二乘法进行线性拟合,分维值结果分别列于表3,分维无尺度区点列拟合效果如图4和图5所示,其中直线斜率的绝对值即为扬州市分形盒维数D。为了对比无尺度区分形维数计算结果,同时对全部点列的对数值也进行了线性拟合,所得结果如表3、图4与图5。对比表3中拟合结果的测定系数R可以得出,全部序列线性拟合度与无尺度区的线性拟合度相比较差,全部序列拟合结果的R远低于无尺度区拟合的R,全部点列拟合计算得到的分维值与无尺度区得到的分维值也有较大差异,由此可得无尺度区的选取对得到精确的分维值具有重要的意义。由图4、图5可以看出通过二阶导数的方法获取城市三维空间形态分维无尺度区点列,在双对数图上线性拟合结果较好,能准确地找到不同标度区的分界点 (无尺度区前面的数据点与后面的数据点对无尺度区回归曲线有一定的偏移)。说明本文给出的通过二阶导数来确定无尺度区的方法,可以很好的确定城市形态的分形特征。由表3可得两年无尺度区点列线性拟合的测定系数R的值均在0.996以上,将此值对比Benguigui在测量特拉维夫分维值时制定的判定二维分维的标准R>0.996[12],计算所得的扬州市三维城市形态分维值符合这个标准。为了进一步验证上述判断,对拟合测定系数最低的2003年工业建筑无尺度区点列,采用负指数函数、对数函数和二次多项式等可能匹配的模型进行拟合,得到的测定系数分别为0.3962,0.2245,0.0363,拟合效果远不如幂指数模型。综上,可以得出,扬州市三维形态及不同职能属性建筑均具有分形特征。

对比两年无尺度区三维分维数D可得,扬州市整体三维空间形态分维数由2003年的2.4161增加为2012年的2.4624,这反映了随着时间的推移,扬州市不断协调其内部三维空间的发展,城市三维建筑在向外围延伸的同时,也进行着旧城改造及城市内部空间的填充,从而使三维空间的利用更为有效和紧凑,使城市建设用地的分布更加均匀。各职能类型建筑中,住宅建筑、工业建筑以及商业与公共建筑的三维分维数均有上升,其中住宅建筑分维值增加最大,但各职能建筑类型的分维数值均低于建筑总体分维数,这与二维城市形态分维的包容原理,城市各个职能类的土地利用形态维数不得高于整个城市形态的分维数相一致[25]。住宅建筑维数呈现出绝对优势的增长,维数的增加表明扬州市居住用地已经摆脱了过去过分集中在旧城和古城的局面,大量向新城与郊区发展,从而使居住用地的空间分布变得均匀,并且旧城内部城中村的改造,高层住宅数量的增加也推动了居住建筑维数的增长。居住建筑维数的增长对同期城市整体三维空间形态维数的增长起了重要作用。工业建筑维数的上升与工业郊区化和工业化发展方式有一定的关系。一方面,扬州市工业企业由新城大量外迁直接促使了工业郊区化的发生发展,使得大部分工业建筑集中分布在近郊的工业园区及经济开发区,但是工业郊区化存在着分散搬迁、二次搬迁等现象,整个市区工业发展的遍地开花现象还是比较严重,这使得2012年工业建筑的维数没有因为集中分布而降低,相对于2003年分维数值反而较大幅度的增加。商业与公共建筑维数增加,主要来自于城市商圈地段高层商业与公共建筑数量的增加如广陵新城商务行政中心高层建筑的建设,京华园商圈高层建筑的建设。

表3 扬州市2003年与2012年不同类型建筑分维计算结果Tab. 3 Evolution of the fractal dimension of different types of buildings in Yangzhou

表3 扬州市2003年与2012年不同类型建筑分维计算结果Tab. 3 Evolution of the fractal dimension of different types of buildings in Yangzhou

图4 2003年扬州城市形态三维分维分析图Fig. 4 The fractal dimension of Yangzhou three-dimensional morphology in 2003

图5 2012年扬州城市形态三维分维分析图Fig. 5 The fractal dimension of Yangzhou three-dimensional morphology in 2012

通过对城市形态分形维数的二维盒维数方法进行扩展,提出了一种城市三维形态的三维盒维数计算方法,以扬州市中心城区为例,详细介绍了城市三维分形维数计算过程,并通过二阶导数方法实现了分形无尺度区自动精确的寻找。试验结果显示: 1二阶导数方法可以有效选取城市分形的无尺度区。2003年及2012年取对数后的无尺度区点列线性拟合优度远高于整体拟合优度,并且无尺度区拟合确定系数R最低值均大于0.996,说明与以往城市分形研究中人工经验判定无尺度区相比,使用该方法能够准确客观地选取分形无尺度区。2城市三维形态具有分形特征,城市三维空间结构在不断演化发展。2003年及2012年无尺度区数据拟合的确定系数R均大于0.996,说明点列线性拟合效果均较好,扬州市城市三维分维特征明显;2012年城市三维空间形态分维数较2003年有了较高的增长,不同职能建筑的三维分维数也均有上升,三维分维数的增加表明城市三维空间利用总体上趋于有效与紧凑。

三维盒维数算法为准确全面地测度城市空间结构和城市形态的演化过程提供了新的方法手段,促进了城市研究向高维扩展,对提高城市三维空间集约利用和空间配置效率具有重要的指导意义。但由于目前城市三维空间形态分维特征的相关研究处于起步阶段,参考较少,对此本文仅仅是进行了探索性的研究,随着研究的进展,以下问题有待于进一步的探讨分析: 1仅通过对扬州市中心城区两个时间节点的数据来验证所提出的方法,无论是在时间的尺度上,还是在样本的广泛性上都有所缺乏,下一步的研究中需要通过大量城市数据样本来检验,以对所提出的理论与方法做进一步补充。2依据Batty等的模拟研究,认为D=1.71左右是二维城市分维期望数值[7],并且大量实证都已经证实了,城市形态的二维分维值波动上升到一定程度会向平均状态回归。对于城市三维分维,发现扬州市三维空间无尺度区整体及各职能建筑的分维值均在2.0~2.5之间,因此猜想,在此之间是否也同样存在一个三维分维期望值,这有待于进一步的理论与实证研究来确定。