1 引 言

元胞自动机(CA)由Ulam在20世纪40年代提出，具有强大的空间运算能力^[1]。文献^[2]中利用CA模拟了一些地理现象后，许多学者不断扩展CA在地理学的应用领域，如城市扩张^{[3, 4, 5, 6, 7]}、山火蔓延^[8]、种群变化^[9]与荒漠化等^[10]，并取得了一系列有意义的研究成果。

转换规则的定义是CA模型的核心。在模拟过程中，这些规则决定了元胞状态的变化。传统CA模型的转换规则中的参数是固定的，它们不随模拟过程中的时间和空间改变而变化。在模拟时间尺度跨度较大的复杂地理过程时，模拟的结果难以反映实际情况^[8]。例如，在城市发展过程中，由于资源环境约束条件变化和政策等的调整，使得城市空间格局的演变规律在不同经济发展阶段有所变化。利用CA进行模拟预测时，静态的参数组合往往难以反映这种非线性地理过程。

动态调整模型参数的方法被认为能够较好地适应非线性地理模拟，目前已有部分学者开展相关研究。文献^[3]提出利用肉眼判断的方法来调整SLEUTH-CA模型参数值。此方法固定其他参数，反复调整某一参数寻求较优模拟结果^[3]。当参数较多时，无法使用该方法，而且该方法有一定的主观性。文献^[3]中又提出了计算搜索算法。利用计算机计算各个参数组合的模拟结果与实际结果的吻合度，确定最优的参数组合^[12]。但这种试验需要高性能工作站运算几百小时。文献^[6]提出案例推理(CBR)的CA模型，模型能够根据新的案例动态更新转换规则。但此方法忽略了模拟结果的误差积累。文献^[12]提出基于数据同化的CA模型。模型考虑到模型误差的积累，但是没有动态调整转换规则中的参数。总体而言，已有的研究难以在模拟过程中动态调整参数和控制误差积累，这是影响CA模型在复杂非线性地理模拟中应用的主要原因。

近年来兴起的数据同化提供了一种能够利用观测数据自动调整并优化模型状态和参数的方法。它能够融合不同来源、不同精度的直接或间接观测数据，并且集成模型和各种观测因子，自动调整模型模拟轨迹^{[13, 14]}。在数据同化的方法中，非线性滤波方法占有极其重要的地位，其中，集合卡尔曼滤波(EnKF)方法使用最为广泛，在参数优化估计中也有较多的应用^{[15, 16]}。文献^[17]对一个二维非线性海风模型进行参数优化估计；文献^[18]提出了双集合卡尔曼滤波方法优化水文模型参数；文献^[19]在三维的云模型中对模型参数进行优化；文献^[20]把联合状态矩阵的方法应用到核污染扩散模型中研究模型参数变化。运用集合卡尔曼滤波方法能够很好地在模拟过程中动态修正模型参数，并且能够在初始参数不准确的情况下快速捕捉参数特征使得被优化的参数收敛于真实值，这为CA模型动态调整模型参数与状态提供了可能。

本文提出了一种基于集合卡尔曼滤波(EnKF)动态优化CA模型参数的方法。该方法先构建由参数和状态组成的联合状态矩阵，再利用观测数据更新联合矩阵，最后把更新后的状态和参数重新输入到模型中进行下一次模拟迭代。文中以广东省东莞市为试验区，模拟其城市扩张，分析了单参数优化情况下参数的有效性，并讨论了在多参数情况下参数收敛情况和城市模拟情景。

在模拟复杂的城市演变过程中，元胞自动机以其“自下而上”、自组织性等优点，被认为是一种非常有效的模拟复杂城市系统的模型^{[3, 4, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]}。与数学模型相比，CA模型具有简洁的形式并能很好地结合栅格GIS数据。通过定义合适的转换规则，便可模拟出城市的发展情景^{[3, 4, 5, 6, 7]}。研究表明，影响城市扩展的变量主要是空间距离变量，如离商业中心的最短距离、离道路的最短距离等，而这些变量对应的参数值(权重)反映了对模型的“贡献”度^{[4, 5, 6, 7, 21]}。文献^[21]运用这些空间变量并提出了多准则的判断(MCE)来表达其他土地类型转变为城市用地的概率

点击浏览公式

点击浏览公式

全局性开发概率Pg中考虑到各种空间距离变量对元胞转化的影响，在微观层次上，邻域影响是一个非常重要的因素。它通过一个3×3的核来计算元胞在空间上的相互影响，其计算为

点击浏览公式

在实际中，还应考虑单元约束条件。有些土地类型，如陡峭的山地等不适宜城市发展，而平缓的地区则很适宜城市发展，这里用连续的条件约束来表达

点击浏览公式

综合考虑了全局发展概率，邻域范围和约束条件后，任意元胞在t时刻发展为城市用地的概率可由下式来表达

点击浏览公式

求出元胞发展概率后，还需要判断此元胞是否发展为城市用地。文中采用阈值法，表达如下

点击浏览公式

MCE-CA模型是模拟城市发展的经典模型。在模型模拟预测时，其假设条件是城市发展模式不随时间变化。实际上，由于政治和经济等因素，城市发展的模式会改变。例如，某城市前期是镇中心扩张发展为主；后期则以沿公路发展为主。在这种情况下，若模型用前期的转换规则来模拟整个时期的城市扩张，模型误差不仅不会消除，还会随着模拟迭代的进行而传递下去。因此，考虑引入集合卡尔曼滤波算法到CA模型中，通过同化观测数据动态更新转换规则中的参数，调整模型运行轨迹，以减少误差积累。

集合卡尔曼滤波(EnKF)是文献^[16]根据随机动力预报理论而提出的一种顺序同化算法。该算法首先用Monte Carlo方法生成初始集合，然后用集合预报的思想对初始集合进行预测得出预测集合，再把预测集合和观测集合输入到集合卡尔曼滤波更新方程中得到更新后预测集合即分析集合，最后用分析集合作初始集合预测到下一个观测时刻，如此循环。其实现步骤如下：

(1) 初始阶段，模型状态集合初始化。

点击浏览公式

2) 预测阶段，将分析集合(或者初始集合)输入到模型中，并预测下一个观测时刻k+1得到的预测状态集合。

点击浏览公式

(3) 分析阶段，融合预测状态集合和观测集合，并根据集合卡尔曼滤波更新公式，求得分析集合X_k+1^a

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

(4) 判断是否到结束时刻，若未到，返回步骤(2)，否则结束。

集合卡尔曼滤波运用集合预报的思想，把误差统计量的预报隐含在扰动的集合变量中，避免了传统卡尔曼滤波等需要专门对误差协方差进行计算的问题，也解决了卡尔曼滤波应用在非线性系统中的切线性近似的问题。集合卡尔曼滤波不仅简化了矩阵计算，还能利用随机样本的预报值估计协方差得到较高的准确度。

数据同化方法能够充分利用观测值等信息提高模型模拟精度。引入数据同化方法到CA模型中，笔者进行了初步探索并提出了基于数据同化的CA模型^[12]。在该模型的基础上，进一步考虑参数组合对模拟结果的影响，故引入集合卡尔曼滤波的联合状态矩阵方法控制由参数引起的模型误差。该方法主要包括如下几个步骤，其流程如图 1。

图 1 动态优化CA模型参数流程图 Fig. 1 Flowchart of parameter optimization in CA model

图选项

首先，根据2.1节的描述构建起CA模型并确定需要进行估计的参数。确定的参数包括各个空间距离变量的系数，如公式(2)所示。本文中包括4个，分别是b_railways、b_expressways、b_towns和b_roads。运用蒙特卡罗方法生成N(N为EnKF中的集合个数)组不同的参数组合。生成随机参数组合包含两个必要部分，参数初始估计值与参数随机扰动噪声。随机参数集合表示如下

点击浏览公式

其次，把生成的N组不同的参数值输入到CA模型中进行城市模拟，可得出N组不同的模拟结果。此时CA模型可以用如下形式描述

点击浏览公式

再次，整合N组参数的模拟结果到状态矩阵中，以计算状态误差协方差。在计算误差协方差前，考虑到CA模型模拟结果只有两种形式(发展与不发展)，不能直接参与计算。所以，本文考虑把研究区分成若干个规则的方格，每个方格包含m×m个像元(元胞)，用城市发展密度参与数据同化运算。每个方格的城市发展密度为

点击浏览公式

最后，采用文献^[17]中的方法把参数等同于状态，整合N组CA模型参数组合及其模拟结果的发展密度到同一个联合状态矩阵X_j里，即把参数b_irailways、b_i^expressways、b_i^towns和b_i^roads作为状态增加到原有的状态矩阵X中，构成一个包含参数和状态的联合状态矩阵X_j。再将此联合状态矩阵X_j代替X输入集合卡尔曼滤波方法的更新方程(9)中。通过同化观测数据可得到优化后的参数值，并将其输入CA模型，以进行下一时刻的模拟迭代，如此循环。

通过重复上述几个步骤便可实现在模拟过程中动态调整模型参数，进而改变CA模型的转换规则。此外，在集合卡尔曼滤波的分析阶段，也可根据更新方程得出的最优分析场并按照文献^[12]的方法对模拟结果进行修正。

东莞市是珠江三角洲城市化发展中速度最快的城市之一。期间，该区域土地利用变化剧烈、城市扩张速度快，而且还经历多次产业布局调整，城市发展表现出复杂的非线性过程。选择该区域具有典型的研究意义。

选用模拟起始数据为1993年东莞市遥感TM分类数据。城市发展的距离空间变量有到铁路的距离、到高速公路的距离、到道路的距离和到镇中心的距离(图 2)。运用SPSS软件中的逻辑回归方法对1993—1995年训练样本进行运算得出的各个距离变量的权重(公式(2)中的b_k)作为初始参数值。为了观察参数调整后模型误差情况以及方便对所提出的方法进行对比验证，选择观测资料基于如下考虑：若引入多景遥感观测资料来观察参数收敛情况，往往会产生这种情况，纵然获取的参数准确，仍难模拟出“飞地”、“开发区”等人为干预的城市发展用地。而且此时模拟误差主要是由CA根本规则——“邻域发展规则”引起的，难以剥离参数误差与邻域规则误差进行比较验证。若使用较短模拟时间，期间的观测资料又非常难获取。考虑到CA模型在模拟短时段城市发展有较高的精度(本文为87.7%)，故文中用模拟的数据代替无法获取的历史观测数据，同时，模拟的观测数据也给参数收敛结果的验证提供了理想的真值。获取观测数据的具体过程如下，首先用SPSS软件中的逻辑回归方法对1993—1997年训练样本进行运算，得出的各个距离变量的权重并作为参数的真值(见表 1)。然后把这些权重输入CA模型中，进行模拟并得出模拟结果，最后把模拟结果转换成发展密度的形式，并选择如图 3所示的30个观测点数据作为观测数据。

图 2 空间距离变量 Fig. 2 Spatial variables of CA model

图选项

图 3 观察点空间分布 Fig. 3 Position of observation

图选项

集合(ensemble size)过大，运算量会相当巨大且对结果改善不明显；集合过小，体现不出统计特征并且误差较大^{[13, 22]}，本文选择集合大小为30。由于传统CA模型的模拟精度约为80％，故设定模型误差为0.2。观测误差取决于观测数据的精度(87.7％)，故设定观测误差为0.13。

除了设置起始模拟数据、初始参数等相关参数外，还需注意生成扰动参数的方差和参数演化控制。扰动参数的方差大小对模型模拟结果的误差统计有极大影响。文中用参数的标准差等于初始误差，初始误差等于初始估计值与真值之差的绝对值。在早期的滤波同化应用中，同化时刻的参数演化常常是将随机扰动加入原参数的集合里生成新的参数集合。此方法往往会由于参数本身的特性导致滤波发散和信息丢失^[18]。文献^[23]提出了核平滑进行参数更新以避免滤波发散与信息丢失。其参数核平滑表示如下

点击浏览公式

为了检验联合状态矩阵方法在CA模拟中的适用性，进行了单参数优化试验。在真实的模型参数优化中，往往是多参数同时优化，所以进一步实现了多参数优化，并分析了参数收敛效果与参数优化前后的城市模拟情景。

在单参数优化过程中，仅仅被优化的参数按公式(15)进行高斯扰动，其他参数保持真值不变。单参数试验中，模型各个参数优化结果如图 4和图 5所示。除了记录参数优化过程，图中还给出了模拟结果的均方根误差(RMS)。在图例中，“未同化”与“同化后”分别表示未进行参数优化与经过参数优化后的结果。虽然各个参数逼近真值的速度有所差别，但是都能够在较少的同化次数后逼近真值，并在真值的较小范围内波动。各个参数在收敛于真实参数值的时间上显示了一定的相似性，都是在10次同化操作前后收敛于各自的真实参数值。由于各个参数对观测数据的敏感程度不同，所以收敛速度有一所差别。从RMS的变化情况可以注意到，同化后的误差均比没有进行同化的要低。参数“普通公路”与“镇中心”的误差积累约为未同化的1/3，而参数“高速公路”与“铁路”误差只有未同化的1/6。通过动态优化参数和状态，模型能够较好地控制模型误差积累。

图 4 优化的参数值与RMS变化(1) Fig. 4 Parameter value and RMS in parameter optimization (1)

图选项

图 5 估计的参数值与RMS变化(2) Fig. 5 Parameter value and RMS in parameter optimization (2)

图选项

在单参数优化试验可以发现，本方法能够在初始参数不准确的情况下快速获得参数特征，并使得被优化的参数收敛于真实值。

为了完整地检验集合卡尔曼滤波在优化模型参数的效果，分别进行了两组双参数和一组四参数同时优化的试验。在双参数优化中，被优化的参数进行高斯扰动，其他参数保持真值。两组双参数优化后的模拟精度见表 2，参数收敛情况如图 6所示。结果表明，模拟精度有一定程度的提高，参数收敛情况与单参数优化相似。经过较少次的更新操作后，各个优化的参数均能较快地逼近真值，并在相对较小的误差范围内波动。同时，模型误差RMS均有较大程度的降低。

图 6 双参数联合优化的结果 Fig. 6 Parameter values and RMS for two spatial variables

图选项

当把优化的参数扩展到所有的参数集合时，各个参数都进行高斯扰动，而且其初始值均偏离真值。通过同化观测资料，笔者得到的结果见图 7与表 2。总体上，各个参数的收敛速度迅速，模型误差较低，模拟精度有所改善。在优化结果中，“高速公路”、“铁路”与“镇中心”3个参数均能较快地收敛，然而参数“普通公路”却一直波动。这主要是由于两个原因：① 与其他参数相比，参数“普通公路”比其他值要大，这就意味着“贡献”大，参数的敏感性强，观测值对参数的变化影响较大，故在一定时间段内出现振荡，需要较多的更新次数来捕捉参数特征；② 在进行数据同化前，将CA模拟结果转换成发展密度，是用均值代替各个元胞状态值，这在一定程度上降低了参数的敏感性。

图 7 四参数联合优化的结果 Fig. 7 Parameter values and RMS for four spatial variables

图选项

对比模拟结果图(图 8)可以发现，未进行参数优化的模拟结果在市区中南部的大岭山、市区东南部寮步以及研究区的东南部塘厦均有模拟过度的情况。这是由于转换规则的误差引起的，东莞市城市发展前期是围绕城镇中心进行扩张发展，后期沿公路进行发展。未进行优化的CA模型因不能在模拟过程中改变转换规则，会一直以围绕城市中心发展的模式模拟城市演化过程，因此，在某些城镇中心会产生模拟过度的情况。引入集合卡尔曼滤波后，模型能够不断修正转换规则中的参数。修正后的参数能够更好地反映真实情况与各个空间距离变量的关系，因此，参数优化后的模拟结果基本上不存在过度发展的情况。模拟结果与真实参数的模拟结果更接近，也与实际情况很吻合。对比RMS与城市模拟结果，优化后的模拟结果均有一定程度的改善。而且各个参数收敛效果良好，在经过较少次的更新优化操作后，均能够在真值相对较小的误差范围内波动。

图 8 四参数优化的模拟结果 Fig. 8 Simulation results of four model parameters optimization

图选项

随着参与优化的参数越多，与单参数优化相比，部分参数需要的观测数据会越多。在图 6和图 7中，某些时刻的参数值会突然偏离真值，这主要是由于观测误差引起的，一些观测数据在某些时刻会由于观测误差而偏离真实值，从而导致同化结果也会有一定的误差。虽然部分参数需要较多的同化次数来收敛于真值，但是从总体上，观察数据的选择和观测误差的估计是合理的，从而获取到准确的信息，较好地反映了参数自动优化的整个过程。

本文采用集合卡尔曼滤波方法更新CA模型参数与模型状态，并以东莞市为试验区进行模拟试验，得到如下结论：

(1) 本文方法通过构建由模型参数和状态组成的联合状态矩阵，并利用观测数据更新此联合状态矩阵，得出更新后的参数与状态，并将其重新输入到模型中进行下一时刻的模拟。该方法有效地实现了模型参数和状态的动态优化，提高了模型在复杂非线性地理模拟中适应能力。

(2) 单参数和多参数联合优化试验表明，该方法能够很好地优化CA模型参数值并快速收敛于真实参数值，而且还能够有效地减少模型误差。在多参数优化试验中，模拟结果的模拟精度与Kappa系数均有一定程度的提高。

(3) 由于城市发展模式的改变，各个空间距离变量在模拟过程中对模拟结果的“贡献”会变化，即，空间变量对CA模型转换概率的影响会改变。同时空间变量对应权重的不同会引起参数敏感度不同，其参数收敛速度会有所差别。

(4) 空间距离变量对模拟结果有重要的影响，在如何选用合适方法以找寻最优空间变量组合和发展密度网格大小对模型参数的收敛等方面需要作进一步研究。