本文为自本人珍藏 版权所有 仅供参考
都市区城市扩展模型 — 以北京城市扩展模拟为例
1 引言
大都市区是指一定规模的中心城市和与中心城市具有紧密社会经济联系的外围地域。当城市规模较小时，城市功能的对应空间主要集中在中心区和建成区；而当核心城市的焦聚超过一定的门槛规模时，核心与周围地域的空间相互作用将产生具有一体化特征的紧密联系区，即大都市区[1]。事实上，从空间上看，中国近20多年的城市化过程主要可以概括为2大类，一是新兴城市的发展兴起，一是原有大城市在改革开放背景下的快速发展，并在部分地区形成大都市区和都市连绵区[2]。后者与前者相比，在中国的城市化进程中也是普遍存在的，而且无论是空间形态，变化过程，驱动机制还是后果影响都比前者更为复杂[1]，同时已有的研究表明，随着中国城市化进程的进一步深入发展，大都市将成为中国城市化进程中最引人注目的地区，同时也是中国未来城市化最具有活力的地区[3]。因此，开展大都市区城市化过程中的土地利用/覆盖变化研究，认识其格局特征和演变过程，并在此基础上进行有效的模拟、预测，从而最大限度地降低城市化过程的风险水平，促使城市化与社会经济的协调发展，已经成为了当前区域土地利用/覆盖变化研究的重要内容[4]。 CA是一种时间、空间、状态都离散， (空间上的) 相互作用和 (时间上的) 因果关系皆局部的格网动力学模型，它“自下而上”的研究思路，强大的复杂计算功能、固有的并行计算能力、高度动态特征以及具有空间概念等特征，使得它在模拟空间复杂系统的时空演变方面具有很强的能力，在地理学研究中具有天然优势[5]。其中，CA在城市增长、扩散和土地利用演化的模拟方面研究最早，最为深入，同时也是当前CA应用的热点[5]。Tobler在20世纪70年代认识到CA在模拟地理复杂现象方面的优势，首次正式采用CA的概念来模拟当时美国5大湖区底特律城市的迅速扩展[5]。随后，
Coucleis的工作使人们看到了应用CA进行城市模拟的巨大潜力[6,7]，引起了人们应用CA开展城市模拟的极大兴趣。Batty、Clarke、White、Wu等人先后开展了相关的研究，并取得了积极的进展[8-18]。如Clarke等根据城市发展的历史数据对美国西海岸的San Francisco和东部的Washingto-Baltimore都市区城市发展的模拟和长期预测[11, 12]；White等应用CA模型对美国Cincinnati市的城市增长、全球气候变化对Caribbean岛的土地利用构成变化影响进行的系列研究[13-16]；Wu.集成CA模型和多因子评价模型对中国广州市城市扩展的模拟研究等[17]，均是CA模型在城市模拟中较为成功的应用案例。此外，受国际研究的推动，国内地理学界近年来也开始了类似的研究，周成虎等于1999年出版了“地理元胞自动机”一书，对相关工作进行了介绍和总结，并提出了地理元胞自动机 (GeoCA) 的概念[5]。Li、Yeh等在对广东东莞土地利用变化系统研究的基础上，利用约束CA模型对广东东莞的土地利用变化进行了成功模拟[19,20]。但从研究内容上看，目前国内已有的工作主要还是集中在对城市发展的初期阶段，尤其是对一些新兴城市的模拟上，针对大都市区土地利用/覆盖变化特征的CA模拟工作还比较少见。因此，在理解大都市区土地利用/覆盖变化机制的基础上，建立和发展基于CA的大都市区城市发展动态模型，无疑具有重要的理论和实践意义。
鉴于以上研究背景，在有关研究的基础上[4, 21]，我们发展了一个基于CA和经济学Tietenberg模型的大都市区城市扩展动态模型 (City Expanding Model in Metropolitan Area; CEM)，对北京地区改革开放以来的城市扩展过程进行系统的模拟和预测。研究的基本目的在于： (1) 利用CA模型模拟中国改革开放背景下大都市区城市发展过程； (2) 重建北京地区改革开放以来的城市发展过程 (1975~1997年)，进一步认识和理解城市扩展的基本过程和机制； (3) 模拟预测北京地区未来不同情景下的城市发展格局 (1998~2015年)，从而为区域可持续发展提供决策支持。
2 大都市区城市扩展的CEM模型
2.1 CEM模型的基本特征
2.1.1 CEM的城市扩展模拟是一个在宏观外部约束因素和局部城市单元自身扩展能力变化因素共同影响作用下的变化过程。 大都市区的城市发展过程必然受到不同尺度的社会、经济和自然因素的综合影响和共同制约，很明显城市的整体格局不可能仅仅用局部小尺度城市单元的相互作用结果来反映[22]。此外，城市发展一般包括向心型城市化和离心型城市化2种过程，前者指城市中的商业服务业设施以及政府部门、企事业公司的总部、银行、报社等脑力劳动机关，不断向城市中心集聚，促使城市中心土地利用密度升高，向立体发展，形成中心商业事务区；后者指部分城市设施和部门，如大型企业、煤气厂等则自城市中心向外缘移动扩散，导致城市外围农村地域变质，城市平面扩大[23]。在城市发展的初期，可利用土地资源一般比较丰富，离心型城市化过程往往占有优势，城市发展以二维平面扩展为主，但随着城市的发展，在有限用地资源限制的影响下，向心型城市化过程逐渐加强，原有城市单元上往往向三维方向发展，常常表现为空间容积率的提高和自身平面扩展能力的下降。这种现象，在城市发展的中后期阶段，尤其是大都市阶段，显然是比较
明显的[23]。然而目前已有的CA城市模拟工作主要把城市发展看成是一定外部约束条件控制下，非城市单元在二维平面上向城市单元的转化[4, 20]。尽管已有的部分工作在利用CA进行虚拟城市研究时对城市单元本身的发展演变能力也进行了定义和考虑[10]，但在实际城市的模拟中，却较少考虑城市单元本身扩展能力的变化[4, 20]。
因此，CEM首先对传统的CA模型进行了进一步的修改，将大都市区的城市发展演变看成是一个宏观外部约束因素和局部城市单元自身扩展能力变化因素共同影响作用下的变化过程，认为大都市区的城市发展演变模拟不仅要考虑各种外部约束因素的影响，而且还要考虑局部城市单元自身平面扩展能力变化因素的影响 (图1a)。对于影响城市发展演变的各种外部约束性因素，根据其对城市发展过程的作用效果，CEM将它们分为一般约束性因素和强制约束性因素2大类。前者一般包括交通状况，到城市中心的耗费距离，坡度坡向等因素，它们对城市的发展只起到一般性的限制作用，后者一般包括湖泊，水库，规划保护用地，河流洪泛区等因素，它们对城市的发展起到绝对性的排斥作用。对于城市单元自身扩展能力的变化，CEM则主要定义了一个随时间变化的指数衰减函数来进行反映。同时，为了把大尺度因素和局部因素 (邻域影响) 综合反映在CA模型的转移规则中，很有必要对这些因素进行统一量化并合理决定它们对城市化过程的相对贡献 (权重)。在这种情况下，各种外在因素的定量化和相应因素权重确定就成为了CEM模型十分关注的问题。鉴于各种影响因素在城市扩展过程中的作用一般随时间和空间的变化而不断改变，而目前的许多研究在考虑这个问题时主要根据经验和专家的建议来确定影响因素的权重，不免显得武断和主观，因此CEM模型在利用模糊关系函数对各种因素进行严格的标准化的基础上发展了一种自适应Monte-Carlo随机方法来确定影响因素的权重并进行城市扩展的模拟。
2.1.2 CEM的城市扩展预测是用地总量最优和位置最佳有机统一，从而达到用地效益的最优的过程。 在中国，土地资源十分有限，城市扩展占用耕地与未来粮食安全的矛盾一直十分突出，因而大都市区的城市扩展实际上是受到政策、经济等多种因素严格控制[2]。鉴于此，在预测思路上，CEM认为大都市区未来城市的发展格局应该是一个用地总量最优和位置最佳的有机统一，从而达到用地效益最优的过程。由于目前自下而上的CA模型对未来城市发展的总量难以有效控制，因此，CEM把微观CA过程和宏观经济学模型结合起来进行未来的城市扩展预测。即把土地资源视为一种不可再生资源，首先使用Tietenberg资源分配模型来完成未来各个时段内城市用地数量上的最优分配[19]，然后再利用CA模型完成位置上的最佳分配 (图1b)。此外，考虑到大都市区的城市扩展过程还存在着各种随机的、不完备信息和不确定性的因素，CEM还在模型中引入Monte-Carlo随机方法来确定有关参数，力图在模型中对这些不确定因素也有所反映。
2.2 CEM模型的基本内容
CEM模型的基本框架可以分成2个部分 (图2)：一是在城市增长历史数据的基础上采用自适应Monte-Carlo方法确定各种影响因素的最佳权重并重建城市扩展过程，二是在CA模型和Tietenberg模型的基础上预测未来大都市区城市发展演变格局。
2.2.1 城市单元自身平面扩展能力变化的影响。 Batty等人进行城市模拟工作时，引入城市土地单元活力值的概念，把城市单元划分为青年、中年和老年，认为城市土地单元有一个从产生到死亡的完整的生命周期过程，在不同的阶段具有不同的演化特征[10]。不过城市单元死亡后变成其它用地单元的情况尽管具有理论上的意义，但在实际的城市发展中，尤其是中国过去20多年的快速城市化过程中，还很少出现[2]。在CEM中，参考他们的工作，借用相关概念来反映城市单元平面扩展能力自身的衰减性。我们将城市土地单元划分为青年城市单元和非青年城市单元。认为城市单元的平面扩展能力随着年龄的变化而衰减，青年城市单元年龄的较小，平面扩展能力强，可以对周围的非城市单元产生影响。非青年的城市单元年龄较大，平面扩展能力弱，对周围的非城市单元则不产生影响。
设一个城市土地单元j产生的时刻为tj，则在时刻ti时，该单元j的扩展衰减值 由下式给出，
Pji (t) = Mexp-?姿 (ti - tj) (1)
式中：M为一标准化常数， 表示土地单元平面扩展能力的的衰减速率。
在得到该城市单元扩展衰减值后，就可以利用Monte-Carlo随机方法来确定该单元所处的城市发展阶段。具体的，设P1代表该单元的扩展衰减值，随机在[0，1]内产生一个实数a, 如果a∈[0, P1]，则认为该单元属于青年，可以对周围的城市单元产生影响，如果a∈[P1, 1]，则认为该单元属于非青年，对周围的非城市单元不产生影响。
2.2.2 基于模糊关系函数的因素标准化。 对一般约束性因素，根据其对城市单元的作用特点，我们首先设计了各种模糊关系函数，将其因素标准化到 [0，100]，然后再采用线性权重组合法得到最终的影响值。对强制性约束因素，则用0、1二值数据表示[4]。
2.2.3 CA邻域的定义。 在城市模拟中，将有更多的邻域单元对中心单元发生影响，因此，有必要对一般的4单元诺依曼 (Neumann) 邻域和8单元 (Moore) 邻域进行扩展。参考有关研究[20]，本模型采用一个以7为半径的标准圆形
邻域。具体的，对中心单元j的邻域影响Nj可以表示为：N = A ×I (2)
式中：drj 为邻域范围到中心单元距离，Ik为一个二值变量，如果邻域单元为青年城市单元，则值取1，否则取0，A是标准化因子。N也采用模糊关系函数进行标准化处理表示在[0，100]内。
2.2.4 确定影响因素权重的自适应Monte-Carlo方法。 一般的，模拟结果总是希望尽可能的与实际的城市扩展相一致，因此各种影响因素权重的确定可以看成是一个在已知城市增长结果的情况下确定权重影响参数的优化问题。假设权重总和为100并且所有的权重是正数，则这个问题变成了约束性的求极值问题，可以表示如下：
约束条件： Wk = 100 (3)
目标函数： MaxF (w1, w2, ..., wm) (4) 式中：Wk为因素 k的权重，F 为表示模拟结果和实际结果之间拟合程度的函数值。我们的目标是寻找最佳的权重以使模拟结果和实际结果之间的拟合函数值F达到最大。这个函数可以利用Monte-Carlo方法进行求解，与AHP方法和Delphi方法相比，Monte-Carlo方法更客观并且可以避免寻找专家打分。基于自适应Monte-Carlo方法的城市扩展模拟具体流程如下：
(1) 根据约束条件，在1到100 + m ?C 1之间随机产生产生正整数 (L1, L2,„,Lm-1)，然后对 L1, L2,„, Lm-1 进行升序排列 (L(1)
W1 = L (1)
W2 = L (2) - L (1)
„„„„„„.
Wm = 100 - L (m-1) (5)
(2) 由式(6)计算出t时刻非城市单元ij向青年城市单元转化的概率。
P = W ×S + W ×NC (6)
式中：W ×S 表示各种一般约束性因素的影响，S 代表对于一般约束因素K 的标准化值，Wk 是该因素的权重；N 代表邻域影响，Wm 是其权重；C 是一系列二值变量的乘积，表示对城市扩展的各种强制性约束因素，如果C = 0，则该单元可能是海、河、湖或其它被保护土地，它们在模型中不能被用作城镇用地。
(3) 在得到转移概率 P 后，就可以根据概率的高低对土地单元ij进行初次分配。
(4) 在完成所有非城市单元的初次转化后，则采用上述定义和方法，计算出各种城市单元的平面扩展衰减值，完成城市单元的再次分配 (5) 重复第 (2) 到第 (4) 步，直到完成该模拟时段的土地利用变化总量为止。
(6) 把第5步得到的模拟结果和检验数据进行比较，计算目标函数值F。在北京的模拟中，检验数据主要来自于由1984，1991和1997年的遥感数据得到的土地利用/覆盖图[21]，F值的定义同遥感精度评价中常用的Kappa系数，如式(7)。 F = (7) 式中：xii表示错误矩阵主对角线上的元素，xi+表示错误矩阵i行的和，x+i示错误矩阵i列的和。
第1到第6步之间需要重复多次以使模拟结果和实际结果尽可能接近，最高F值对应的权重就是各因素的最佳权重，而此时的模拟结果就是城市扩展的最佳模拟结果。这里，我们需要确定重复的次数。根据Miyatake and Wakimoto的工作，利用Monte-Carlo方法当重复次数达到500次时，模拟极值逼近实际极值的可能性可以达到0.99[24]。因此，我们在第1步到第6步之间重复500次以确保可以得到可靠的权重和模拟结果。
2.2.5 未来城市最佳用地总量分配 对未来城市格局的预测，主要应该满足城市用地数量最优和位置最佳2个条件。参考有关工作，我们采用Tietenberg模型来确定未来各个时段的最佳城市用地土总量[19]。Tietenberg模型是一个不可再生资源的动态时间分配模型，该模型把涉及时间因素的贴现率放进了模型中，认为对于一定量的不可再生资源，问题是如何在时间上安排它的使用，以获得最大的收益。在n年内最有效地分配Q总量的资源应该满足如下的最大值条件：
(aq - bq - cq)(1 + r)+ ?姿 (Q - q) (8)
式中：Q是所提供的资源总量。a是边际收益曲线的截距，即边际收益曲线的最大理论值。b 是边际收益曲线的斜率，可以选为1，c为边际费用的常数，其值比a小，可以选为c = a/2，r是贴现率，t是时间，?姿是极值公式的常数。 Yeh等采用该模型来进行土地利用总量的分配，建立了下面的方程[19]：
a - bqt / Pta - c ) / (1 + r)t-1 - ?姿 = 0
t = 1, ..., n
Q - qt = 0 (9)
式中：Pta是t时期的增加人口，qt是所对应的用地量，Q是分配的土地总量，其它参数的意义同式(8)。
3 北京地区的城市扩展模拟案例
3.1 模拟区域和使用的数据
具体模拟区域主要包括北京城区的东城、西城、宣武、崇文；近郊区的石景山、海淀、朝阳、丰台和远郊区的昌平、顺义、通县等11个区县级行政单元，范围为115o50′E~116059′E， 39036′N~40023′N。该区1998年人口919.4×104人，面积4649.9 km2，分别占北京市的74 %和28 %；地形上西北高东南低，由西北向东南呈现出低地-丘陵-山前洪积-平原区的有序排列；经济上具有从城市核心区、城乡过渡区到远郊区县的明显过渡，映射出人类活动由强到弱的梯度变化；空间上则呈现从中央大区、城市边缘区到外围地域的明显圈层变化并且整体上联系紧密，表现出大都市区的基本特征[1]。模型中使用的土地利用/覆盖数据主要是来自于编号为123/32的4期Landsat TM/MSS (1975年5月6日获取的MSS以及1984年10月2日，1991年5月6日和1997年5月16日分别获取的TM) 数据[21]，由于遥感影响覆盖能力限制的原因,实际模拟区域总面积为4499.57 km2，像元大小为150m×150m。
3.2 1975~1997城市发展过程重建和模拟
结合北京地区实际情况，主要考虑了9个对区域城市发展演变过程起作用的外部约束因素，其中包括6个一般性约束因素 (对铁路、高速公路、一级公路、首都机场、城市中心的耗费距离和坡度) 和3个强制性约束因素 (河流洪泛区、城市规划中的保护绿地、各种面状水体 (水库、人工湖泊))。首先利用CEM模型对北京地区1975~1984年，1984~1991年，1991~1997年3个时段的城市发展过程进行了模拟重建(图3)，并得到各个时期各种约束因素的最佳模拟权重(表1)。
据表1可见：(1) 各种限制因素的影响能力随着城市的发展在不断发生变化的，现有的许多CA城市模型使用单一影响权重来进行城市模拟显然是不符合实际情况的。(2) 城市单元邻域作用影响的绝对值远远大于其它因素，分别达到了69，60和53，由此可见模拟时间内大都市区旧有城市格局对城市发展的影响是非常巨大的。(3) 交通状况的影响总体上表现出上升趋势，从16上升到了27，其中尤以高速公路的影响上升明显，而同期邻域的影响作用则表现出比较大的下降趋势，可见随着经济的发展和城市单元自身扩展能力的下降，大都市区的城市格局也在一定程度上发生着变化。
同时，1975~1984年，1984~1991年，1991~1997年模拟结果和实际遥感测量结果之间的kappa系数分别达到了0.59、0.65、0.67，这说明CEM可以在一定程度上反映大都市区城市发展演变的基本特征和规律。
模拟结果和实际结果的差异主要表现在一下几个方面：一是模拟结果的中央大区表现出放射状的星型趋势，而实际的中央大区主要还是一种“饼状”的圆形；二是模拟结果中次级中心如昌平，顺义和通县的发展明显高于实际的城市发展；三是模拟结果中对中心大区和次级中心之间的许多中小城镇没有很好的反映。其原因可能主要在于模拟时夸大了区域线状交通状况对旧有城市格局的影响，而CA自身模拟的局限性可能导致了模拟结果对许多点状的中小城镇没有很好的反映。实际上，基于局部自组织规则的CA模型对于许多宏观政策性因素难于有效的反映，如政府20世纪90年代对区域中小城镇发展的鼓励措施等。
3.3 1998~2015年城市发展格局动态预测
由于2-9式分配未来的城市用地总量需要未来的人口数量，考虑到历史资料的可靠性和有关政策背景，采用1980~1998年的区域人口统计资料，依据Logistic模型，从式(10)用最小二乘法拟合出式(11)，并据此完成区域2001~2015年增加人口的预测 (表2)。
X (t) = (10) 式中: X (t) 表示t时间的人口数量，Xm代表人口增长率为0时的人口总量，即人口承载力，r是当人口增量为0时的增长率，称为固有增长率，X0表示预测开始时的人口总量。
y (t) = 1220.05419/(1 + 0.9019*exp(-0.0418*t)) (t = 0, 1, 2, 3 „) (11)
式中：y (t) 表示总人口，它是常住人口和流动人口50 %折算后的总和[19]，t = 0 代表1980年。
以1997年遥感监测数据为基础，可知区域总面积为4499.57 km2，1997年城镇用地总量为1023.55 km2，占研究区域总面积的22.5 %[21]。由于区域城镇用地的增加一方面将导致耕地、园地，尤其是平原区耕地的大量减少，从而影响区域的粮食安全，另一方面将影响区域的生态调节功能，改变区域的生态安全格局，产生一系列生态环境问题，所以会受到政府严格控制。假定政府实行不同的城市化政策，到2015年，区域城镇用地可以分别占到区域总面积的25 %，30 %和35 %，那么区域城市用地增加的总量将分别为101.34 km2, 326.32 km2和551.30 km2。以1997年城镇用地为基础，则城镇用地的年增加率将分别达到0.5 %，1.7 %和2.9 %。考虑1975~1984年，1984~1991年，1991~1997年区域城镇用地的年增加率曾分别达到13.6 %，4.4 %和5.0 %，可以认为这3种速度事实上都是一种严格限制下的城市发展，这也从另一个方面说明了由于1975~1997年的快速城市化过程，区域城市通过进一步占用耕地
平面扩展的潜力事实上也是有限的。
进一步采用Tietenberg模型，依据式2-9，可以首先得到未来各个时期的城镇用地最佳使用总量 (表2)。最后以表2为基础，采用CEM模型，从空间上进行未来城市用地的分配，从而完成对北京未来城市动态变化过程的预测 (图4)。 从2005~2015年的模拟结果可以发现： (1) 中心大区放射性星状格局基本形成，石景山、城区至通县一线，城区、首都机场和顺义一线基本上已经形成了城市廊道，城区至昌平一线也显示出了形成城市廊道的趋势。 (2) 通县，昌平和顺义等区域次级中心都得到了进一步的充分发展，尤其是昌平，发展十分迅速，在中心大区扩展能力减弱的情况下，与周围城镇连成一片，成为了区域的另一个发展中心。
4 结论和讨论
在分析大都市区城市扩展特征的基础上，从宏观外部约束性因素和局部城市单元自身扩展能力变化共同作用影响城市发展演变的角度，构建了一个模拟和预测大都市区城市发展演变过程的城市扩展模型 (CEM)。首先对北京地区1975~1997年的城市发展过程进行了模拟重建，然后在此基础上从城镇用地数量最优和位置最佳两个方面对其2005~2015年的城市发展格局进行了预测。结果表明，在对各种影响因素进行严格标准化并利用自适应Monte-Carlo方法多次模拟确定最佳影响权重的基础上，该模型可以在一定程度上反映城市发展，尤其是大都市区城市发展演变的特征和规律。
通过对北京城市发展变化过程进行重建，我们发现各种限制因素的影响能力随着城市的发展在不断发生变化，目前许多的城市发展CA模型对各种影响因素仅仅分配单一权重的作法明显是不符合实际情况的。同时我们也发现，模拟时间内大都市区旧有城市格局对城市的发展的影响非常巨大，表现出明显不同于新兴城市的发展变化特征。
模型对北京地区未来的城市发展预测，实际上提供了一种相对比较理想的城市发展格局。从模拟结果上看，加强中心大区城市内部改造，鼓励城市立体发展，提高城市容积率；加强区域线状交通体系建设，鼓励远郊区域次级中心城市发展，有可能在一定程度上改变目前北京地区人口与产业过分集中于市区，中心城区膨胀，摊大饼外延发展的城市发展局面。
本文为自本人珍藏 版权所有 仅供参考
都市区城市扩展模型 — 以北京城市扩展模拟为例
1 引言
大都市区是指一定规模的中心城市和与中心城市具有紧密社会经济联系的外围地域。当城市规模较小时，城市功能的对应空间主要集中在中心区和建成区；而当核心城市的焦聚超过一定的门槛规模时，核心与周围地域的空间相互作用将产生具有一体化特征的紧密联系区，即大都市区[1]。事实上，从空间上看，中国近20多年的城市化过程主要可以概括为2大类，一是新兴城市的发展兴起，一是原有大城市在改革开放背景下的快速发展，并在部分地区形成大都市区和都市连绵区[2]。后者与前者相比，在中国的城市化进程中也是普遍存在的，而且无论是空间形态，变化过程，驱动机制还是后果影响都比前者更为复杂[1]，同时已有的研究表明，随着中国城市化进程的进一步深入发展，大都市将成为中国城市化进程中最引人注目的地区，同时也是中国未来城市化最具有活力的地区[3]。因此，开展大都市区城市化过程中的土地利用/覆盖变化研究，认识其格局特征和演变过程，并在此基础上进行有效的模拟、预测，从而最大限度地降低城市化过程的风险水平，促使城市化与社会经济的协调发展，已经成为了当前区域土地利用/覆盖变化研究的重要内容[4]。 CA是一种时间、空间、状态都离散， (空间上的) 相互作用和 (时间上的) 因果关系皆局部的格网动力学模型，它“自下而上”的研究思路，强大的复杂计算功能、固有的并行计算能力、高度动态特征以及具有空间概念等特征，使得它在模拟空间复杂系统的时空演变方面具有很强的能力，在地理学研究中具有天然优势[5]。其中，CA在城市增长、扩散和土地利用演化的模拟方面研究最早，最为深入，同时也是当前CA应用的热点[5]。Tobler在20世纪70年代认识到CA在模拟地理复杂现象方面的优势，首次正式采用CA的概念来模拟当时美国5大湖区底特律城市的迅速扩展[5]。随后，
Coucleis的工作使人们看到了应用CA进行城市模拟的巨大潜力[6,7]，引起了人们应用CA开展城市模拟的极大兴趣。Batty、Clarke、White、Wu等人先后开展了相关的研究，并取得了积极的进展[8-18]。如Clarke等根据城市发展的历史数据对美国西海岸的San Francisco和东部的Washingto-Baltimore都市区城市发展的模拟和长期预测[11, 12]；White等应用CA模型对美国Cincinnati市的城市增长、全球气候变化对Caribbean岛的土地利用构成变化影响进行的系列研究[13-16]；Wu.集成CA模型和多因子评价模型对中国广州市城市扩展的模拟研究等[17]，均是CA模型在城市模拟中较为成功的应用案例。此外，受国际研究的推动，国内地理学界近年来也开始了类似的研究，周成虎等于1999年出版了“地理元胞自动机”一书，对相关工作进行了介绍和总结，并提出了地理元胞自动机 (GeoCA) 的概念[5]。Li、Yeh等在对广东东莞土地利用变化系统研究的基础上，利用约束CA模型对广东东莞的土地利用变化进行了成功模拟[19,20]。但从研究内容上看，目前国内已有的工作主要还是集中在对城市发展的初期阶段，尤其是对一些新兴城市的模拟上，针对大都市区土地利用/覆盖变化特征的CA模拟工作还比较少见。因此，在理解大都市区土地利用/覆盖变化机制的基础上，建立和发展基于CA的大都市区城市发展动态模型，无疑具有重要的理论和实践意义。
鉴于以上研究背景，在有关研究的基础上[4, 21]，我们发展了一个基于CA和经济学Tietenberg模型的大都市区城市扩展动态模型 (City Expanding Model in Metropolitan Area; CEM)，对北京地区改革开放以来的城市扩展过程进行系统的模拟和预测。研究的基本目的在于： (1) 利用CA模型模拟中国改革开放背景下大都市区城市发展过程； (2) 重建北京地区改革开放以来的城市发展过程 (1975~1997年)，进一步认识和理解城市扩展的基本过程和机制； (3) 模拟预测北京地区未来不同情景下的城市发展格局 (1998~2015年)，从而为区域可持续发展提供决策支持。
2 大都市区城市扩展的CEM模型
2.1 CEM模型的基本特征
2.1.1 CEM的城市扩展模拟是一个在宏观外部约束因素和局部城市单元自身扩展能力变化因素共同影响作用下的变化过程。 大都市区的城市发展过程必然受到不同尺度的社会、经济和自然因素的综合影响和共同制约，很明显城市的整体格局不可能仅仅用局部小尺度城市单元的相互作用结果来反映[22]。此外，城市发展一般包括向心型城市化和离心型城市化2种过程，前者指城市中的商业服务业设施以及政府部门、企事业公司的总部、银行、报社等脑力劳动机关，不断向城市中心集聚，促使城市中心土地利用密度升高，向立体发展，形成中心商业事务区；后者指部分城市设施和部门，如大型企业、煤气厂等则自城市中心向外缘移动扩散，导致城市外围农村地域变质，城市平面扩大[23]。在城市发展的初期，可利用土地资源一般比较丰富，离心型城市化过程往往占有优势，城市发展以二维平面扩展为主，但随着城市的发展，在有限用地资源限制的影响下，向心型城市化过程逐渐加强，原有城市单元上往往向三维方向发展，常常表现为空间容积率的提高和自身平面扩展能力的下降。这种现象，在城市发展的中后期阶段，尤其是大都市阶段，显然是比较
明显的[23]。然而目前已有的CA城市模拟工作主要把城市发展看成是一定外部约束条件控制下，非城市单元在二维平面上向城市单元的转化[4, 20]。尽管已有的部分工作在利用CA进行虚拟城市研究时对城市单元本身的发展演变能力也进行了定义和考虑[10]，但在实际城市的模拟中，却较少考虑城市单元本身扩展能力的变化[4, 20]。
因此，CEM首先对传统的CA模型进行了进一步的修改，将大都市区的城市发展演变看成是一个宏观外部约束因素和局部城市单元自身扩展能力变化因素共同影响作用下的变化过程，认为大都市区的城市发展演变模拟不仅要考虑各种外部约束因素的影响，而且还要考虑局部城市单元自身平面扩展能力变化因素的影响 (图1a)。对于影响城市发展演变的各种外部约束性因素，根据其对城市发展过程的作用效果，CEM将它们分为一般约束性因素和强制约束性因素2大类。前者一般包括交通状况，到城市中心的耗费距离，坡度坡向等因素，它们对城市的发展只起到一般性的限制作用，后者一般包括湖泊，水库，规划保护用地，河流洪泛区等因素，它们对城市的发展起到绝对性的排斥作用。对于城市单元自身扩展能力的变化，CEM则主要定义了一个随时间变化的指数衰减函数来进行反映。同时，为了把大尺度因素和局部因素 (邻域影响) 综合反映在CA模型的转移规则中，很有必要对这些因素进行统一量化并合理决定它们对城市化过程的相对贡献 (权重)。在这种情况下，各种外在因素的定量化和相应因素权重确定就成为了CEM模型十分关注的问题。鉴于各种影响因素在城市扩展过程中的作用一般随时间和空间的变化而不断改变，而目前的许多研究在考虑这个问题时主要根据经验和专家的建议来确定影响因素的权重，不免显得武断和主观，因此CEM模型在利用模糊关系函数对各种因素进行严格的标准化的基础上发展了一种自适应Monte-Carlo随机方法来确定影响因素的权重并进行城市扩展的模拟。
2.1.2 CEM的城市扩展预测是用地总量最优和位置最佳有机统一，从而达到用地效益的最优的过程。 在中国，土地资源十分有限，城市扩展占用耕地与未来粮食安全的矛盾一直十分突出，因而大都市区的城市扩展实际上是受到政策、经济等多种因素严格控制[2]。鉴于此，在预测思路上，CEM认为大都市区未来城市的发展格局应该是一个用地总量最优和位置最佳的有机统一，从而达到用地效益最优的过程。由于目前自下而上的CA模型对未来城市发展的总量难以有效控制，因此，CEM把微观CA过程和宏观经济学模型结合起来进行未来的城市扩展预测。即把土地资源视为一种不可再生资源，首先使用Tietenberg资源分配模型来完成未来各个时段内城市用地数量上的最优分配[19]，然后再利用CA模型完成位置上的最佳分配 (图1b)。此外，考虑到大都市区的城市扩展过程还存在着各种随机的、不完备信息和不确定性的因素，CEM还在模型中引入Monte-Carlo随机方法来确定有关参数，力图在模型中对这些不确定因素也有所反映。
2.2 CEM模型的基本内容
CEM模型的基本框架可以分成2个部分 (图2)：一是在城市增长历史数据的基础上采用自适应Monte-Carlo方法确定各种影响因素的最佳权重并重建城市扩展过程，二是在CA模型和Tietenberg模型的基础上预测未来大都市区城市发展演变格局。
2.2.1 城市单元自身平面扩展能力变化的影响。 Batty等人进行城市模拟工作时，引入城市土地单元活力值的概念，把城市单元划分为青年、中年和老年，认为城市土地单元有一个从产生到死亡的完整的生命周期过程，在不同的阶段具有不同的演化特征[10]。不过城市单元死亡后变成其它用地单元的情况尽管具有理论上的意义，但在实际的城市发展中，尤其是中国过去20多年的快速城市化过程中，还很少出现[2]。在CEM中，参考他们的工作，借用相关概念来反映城市单元平面扩展能力自身的衰减性。我们将城市土地单元划分为青年城市单元和非青年城市单元。认为城市单元的平面扩展能力随着年龄的变化而衰减，青年城市单元年龄的较小，平面扩展能力强，可以对周围的非城市单元产生影响。非青年的城市单元年龄较大，平面扩展能力弱，对周围的非城市单元则不产生影响。
设一个城市土地单元j产生的时刻为tj，则在时刻ti时，该单元j的扩展衰减值 由下式给出，
Pji (t) = Mexp-?姿 (ti - tj) (1)
式中：M为一标准化常数， 表示土地单元平面扩展能力的的衰减速率。
在得到该城市单元扩展衰减值后，就可以利用Monte-Carlo随机方法来确定该单元所处的城市发展阶段。具体的，设P1代表该单元的扩展衰减值，随机在[0，1]内产生一个实数a, 如果a∈[0, P1]，则认为该单元属于青年，可以对周围的城市单元产生影响，如果a∈[P1, 1]，则认为该单元属于非青年，对周围的非城市单元不产生影响。
2.2.2 基于模糊关系函数的因素标准化。 对一般约束性因素，根据其对城市单元的作用特点，我们首先设计了各种模糊关系函数，将其因素标准化到 [0，100]，然后再采用线性权重组合法得到最终的影响值。对强制性约束因素，则用0、1二值数据表示[4]。
2.2.3 CA邻域的定义。 在城市模拟中，将有更多的邻域单元对中心单元发生影响，因此，有必要对一般的4单元诺依曼 (Neumann) 邻域和8单元 (Moore) 邻域进行扩展。参考有关研究[20]，本模型采用一个以7为半径的标准圆形
邻域。具体的，对中心单元j的邻域影响Nj可以表示为：N = A ×I (2)
式中：drj 为邻域范围到中心单元距离，Ik为一个二值变量，如果邻域单元为青年城市单元，则值取1，否则取0，A是标准化因子。N也采用模糊关系函数进行标准化处理表示在[0，100]内。
2.2.4 确定影响因素权重的自适应Monte-Carlo方法。 一般的，模拟结果总是希望尽可能的与实际的城市扩展相一致，因此各种影响因素权重的确定可以看成是一个在已知城市增长结果的情况下确定权重影响参数的优化问题。假设权重总和为100并且所有的权重是正数，则这个问题变成了约束性的求极值问题，可以表示如下：
约束条件： Wk = 100 (3)
目标函数： MaxF (w1, w2, ..., wm) (4) 式中：Wk为因素 k的权重，F 为表示模拟结果和实际结果之间拟合程度的函数值。我们的目标是寻找最佳的权重以使模拟结果和实际结果之间的拟合函数值F达到最大。这个函数可以利用Monte-Carlo方法进行求解，与AHP方法和Delphi方法相比，Monte-Carlo方法更客观并且可以避免寻找专家打分。基于自适应Monte-Carlo方法的城市扩展模拟具体流程如下：
(1) 根据约束条件，在1到100 + m ?C 1之间随机产生产生正整数 (L1, L2,„,Lm-1)，然后对 L1, L2,„, Lm-1 进行升序排列 (L(1)
W1 = L (1)
W2 = L (2) - L (1)
„„„„„„.
Wm = 100 - L (m-1) (5)
(2) 由式(6)计算出t时刻非城市单元ij向青年城市单元转化的概率。
P = W ×S + W ×NC (6)
式中：W ×S 表示各种一般约束性因素的影响，S 代表对于一般约束因素K 的标准化值，Wk 是该因素的权重；N 代表邻域影响，Wm 是其权重；C 是一系列二值变量的乘积，表示对城市扩展的各种强制性约束因素，如果C = 0，则该单元可能是海、河、湖或其它被保护土地，它们在模型中不能被用作城镇用地。
(3) 在得到转移概率 P 后，就可以根据概率的高低对土地单元ij进行初次分配。
(4) 在完成所有非城市单元的初次转化后，则采用上述定义和方法，计算出各种城市单元的平面扩展衰减值，完成城市单元的再次分配 (5) 重复第 (2) 到第 (4) 步，直到完成该模拟时段的土地利用变化总量为止。
(6) 把第5步得到的模拟结果和检验数据进行比较，计算目标函数值F。在北京的模拟中，检验数据主要来自于由1984，1991和1997年的遥感数据得到的土地利用/覆盖图[21]，F值的定义同遥感精度评价中常用的Kappa系数，如式(7)。 F = (7) 式中：xii表示错误矩阵主对角线上的元素，xi+表示错误矩阵i行的和，x+i示错误矩阵i列的和。
第1到第6步之间需要重复多次以使模拟结果和实际结果尽可能接近，最高F值对应的权重就是各因素的最佳权重，而此时的模拟结果就是城市扩展的最佳模拟结果。这里，我们需要确定重复的次数。根据Miyatake and Wakimoto的工作，利用Monte-Carlo方法当重复次数达到500次时，模拟极值逼近实际极值的可能性可以达到0.99[24]。因此，我们在第1步到第6步之间重复500次以确保可以得到可靠的权重和模拟结果。
2.2.5 未来城市最佳用地总量分配 对未来城市格局的预测，主要应该满足城市用地数量最优和位置最佳2个条件。参考有关工作，我们采用Tietenberg模型来确定未来各个时段的最佳城市用地土总量[19]。Tietenberg模型是一个不可再生资源的动态时间分配模型，该模型把涉及时间因素的贴现率放进了模型中，认为对于一定量的不可再生资源，问题是如何在时间上安排它的使用，以获得最大的收益。在n年内最有效地分配Q总量的资源应该满足如下的最大值条件：
(aq - bq - cq)(1 + r)+ ?姿 (Q - q) (8)
式中：Q是所提供的资源总量。a是边际收益曲线的截距，即边际收益曲线的最大理论值。b 是边际收益曲线的斜率，可以选为1，c为边际费用的常数，其值比a小，可以选为c = a/2，r是贴现率，t是时间，?姿是极值公式的常数。 Yeh等采用该模型来进行土地利用总量的分配，建立了下面的方程[19]：
a - bqt / Pta - c ) / (1 + r)t-1 - ?姿 = 0
t = 1, ..., n
Q - qt = 0 (9)
式中：Pta是t时期的增加人口，qt是所对应的用地量，Q是分配的土地总量，其它参数的意义同式(8)。
3 北京地区的城市扩展模拟案例
3.1 模拟区域和使用的数据
具体模拟区域主要包括北京城区的东城、西城、宣武、崇文；近郊区的石景山、海淀、朝阳、丰台和远郊区的昌平、顺义、通县等11个区县级行政单元，范围为115o50′E~116059′E， 39036′N~40023′N。该区1998年人口919.4×104人，面积4649.9 km2，分别占北京市的74 %和28 %；地形上西北高东南低，由西北向东南呈现出低地-丘陵-山前洪积-平原区的有序排列；经济上具有从城市核心区、城乡过渡区到远郊区县的明显过渡，映射出人类活动由强到弱的梯度变化；空间上则呈现从中央大区、城市边缘区到外围地域的明显圈层变化并且整体上联系紧密，表现出大都市区的基本特征[1]。模型中使用的土地利用/覆盖数据主要是来自于编号为123/32的4期Landsat TM/MSS (1975年5月6日获取的MSS以及1984年10月2日，1991年5月6日和1997年5月16日分别获取的TM) 数据[21]，由于遥感影响覆盖能力限制的原因,实际模拟区域总面积为4499.57 km2，像元大小为150m×150m。
3.2 1975~1997城市发展过程重建和模拟
结合北京地区实际情况，主要考虑了9个对区域城市发展演变过程起作用的外部约束因素，其中包括6个一般性约束因素 (对铁路、高速公路、一级公路、首都机场、城市中心的耗费距离和坡度) 和3个强制性约束因素 (河流洪泛区、城市规划中的保护绿地、各种面状水体 (水库、人工湖泊))。首先利用CEM模型对北京地区1975~1984年，1984~1991年，1991~1997年3个时段的城市发展过程进行了模拟重建(图3)，并得到各个时期各种约束因素的最佳模拟权重(表1)。
据表1可见：(1) 各种限制因素的影响能力随着城市的发展在不断发生变化的，现有的许多CA城市模型使用单一影响权重来进行城市模拟显然是不符合实际情况的。(2) 城市单元邻域作用影响的绝对值远远大于其它因素，分别达到了69，60和53，由此可见模拟时间内大都市区旧有城市格局对城市发展的影响是非常巨大的。(3) 交通状况的影响总体上表现出上升趋势，从16上升到了27，其中尤以高速公路的影响上升明显，而同期邻域的影响作用则表现出比较大的下降趋势，可见随着经济的发展和城市单元自身扩展能力的下降，大都市区的城市格局也在一定程度上发生着变化。
同时，1975~1984年，1984~1991年，1991~1997年模拟结果和实际遥感测量结果之间的kappa系数分别达到了0.59、0.65、0.67，这说明CEM可以在一定程度上反映大都市区城市发展演变的基本特征和规律。
模拟结果和实际结果的差异主要表现在一下几个方面：一是模拟结果的中央大区表现出放射状的星型趋势，而实际的中央大区主要还是一种“饼状”的圆形；二是模拟结果中次级中心如昌平，顺义和通县的发展明显高于实际的城市发展；三是模拟结果中对中心大区和次级中心之间的许多中小城镇没有很好的反映。其原因可能主要在于模拟时夸大了区域线状交通状况对旧有城市格局的影响，而CA自身模拟的局限性可能导致了模拟结果对许多点状的中小城镇没有很好的反映。实际上，基于局部自组织规则的CA模型对于许多宏观政策性因素难于有效的反映，如政府20世纪90年代对区域中小城镇发展的鼓励措施等。
3.3 1998~2015年城市发展格局动态预测
由于2-9式分配未来的城市用地总量需要未来的人口数量，考虑到历史资料的可靠性和有关政策背景，采用1980~1998年的区域人口统计资料，依据Logistic模型，从式(10)用最小二乘法拟合出式(11)，并据此完成区域2001~2015年增加人口的预测 (表2)。
X (t) = (10) 式中: X (t) 表示t时间的人口数量，Xm代表人口增长率为0时的人口总量，即人口承载力，r是当人口增量为0时的增长率，称为固有增长率，X0表示预测开始时的人口总量。
y (t) = 1220.05419/(1 + 0.9019*exp(-0.0418*t)) (t = 0, 1, 2, 3 „) (11)
式中：y (t) 表示总人口，它是常住人口和流动人口50 %折算后的总和[19]，t = 0 代表1980年。
以1997年遥感监测数据为基础，可知区域总面积为4499.57 km2，1997年城镇用地总量为1023.55 km2，占研究区域总面积的22.5 %[21]。由于区域城镇用地的增加一方面将导致耕地、园地，尤其是平原区耕地的大量减少，从而影响区域的粮食安全，另一方面将影响区域的生态调节功能，改变区域的生态安全格局，产生一系列生态环境问题，所以会受到政府严格控制。假定政府实行不同的城市化政策，到2015年，区域城镇用地可以分别占到区域总面积的25 %，30 %和35 %，那么区域城市用地增加的总量将分别为101.34 km2, 326.32 km2和551.30 km2。以1997年城镇用地为基础，则城镇用地的年增加率将分别达到0.5 %，1.7 %和2.9 %。考虑1975~1984年，1984~1991年，1991~1997年区域城镇用地的年增加率曾分别达到13.6 %，4.4 %和5.0 %，可以认为这3种速度事实上都是一种严格限制下的城市发展，这也从另一个方面说明了由于1975~1997年的快速城市化过程，区域城市通过进一步占用耕地
平面扩展的潜力事实上也是有限的。
进一步采用Tietenberg模型，依据式2-9，可以首先得到未来各个时期的城镇用地最佳使用总量 (表2)。最后以表2为基础，采用CEM模型，从空间上进行未来城市用地的分配，从而完成对北京未来城市动态变化过程的预测 (图4)。 从2005~2015年的模拟结果可以发现： (1) 中心大区放射性星状格局基本形成，石景山、城区至通县一线，城区、首都机场和顺义一线基本上已经形成了城市廊道，城区至昌平一线也显示出了形成城市廊道的趋势。 (2) 通县，昌平和顺义等区域次级中心都得到了进一步的充分发展，尤其是昌平，发展十分迅速，在中心大区扩展能力减弱的情况下，与周围城镇连成一片，成为了区域的另一个发展中心。
4 结论和讨论
在分析大都市区城市扩展特征的基础上，从宏观外部约束性因素和局部城市单元自身扩展能力变化共同作用影响城市发展演变的角度，构建了一个模拟和预测大都市区城市发展演变过程的城市扩展模型 (CEM)。首先对北京地区1975~1997年的城市发展过程进行了模拟重建，然后在此基础上从城镇用地数量最优和位置最佳两个方面对其2005~2015年的城市发展格局进行了预测。结果表明，在对各种影响因素进行严格标准化并利用自适应Monte-Carlo方法多次模拟确定最佳影响权重的基础上，该模型可以在一定程度上反映城市发展，尤其是大都市区城市发展演变的特征和规律。
通过对北京城市发展变化过程进行重建，我们发现各种限制因素的影响能力随着城市的发展在不断发生变化，目前许多的城市发展CA模型对各种影响因素仅仅分配单一权重的作法明显是不符合实际情况的。同时我们也发现，模拟时间内大都市区旧有城市格局对城市的发展的影响非常巨大，表现出明显不同于新兴城市的发展变化特征。
模型对北京地区未来的城市发展预测，实际上提供了一种相对比较理想的城市发展格局。从模拟结果上看，加强中心大区城市内部改造，鼓励城市立体发展，提高城市容积率；加强区域线状交通体系建设，鼓励远郊区域次级中心城市发展，有可能在一定程度上改变目前北京地区人口与产业过分集中于市区，中心城区膨胀，摊大饼外延发展的城市发展局面。