張亦漢，喬紀(jì)綱，艾 彬

1.廣東商學(xué)院 資源與環(huán)境學(xué)院，廣東 廣州510230；2.中山大學(xué) 地理科學(xué)與規(guī)劃學(xué)院，廣東 廣州510275

1 引 言

元胞自動機（CA）由Ulam在20世紀(jì)40年代提出，具有強大的空間運算能力［1］。文獻［2］中利用CA模擬了一些地理現(xiàn)象后，許多學(xué)者不斷擴展CA在地理學(xué)的應(yīng)用領(lǐng)域，如城市擴張［3－7］、山火蔓延［8］、種群變化［9］與荒漠化等［10］，并取得了一系列有意義的研究成果。

轉(zhuǎn)換規(guī)則的定義是CA模型的核心。在模擬過程中，這些規(guī)則決定了元胞狀態(tài)的變化。傳統(tǒng)CA模型的轉(zhuǎn)換規(guī)則中的參數(shù)是固定的，它們不隨模擬過程中的時間和空間改變而變化。在模擬時間尺度跨度較大的復(fù)雜地理過程時，模擬的結(jié)果難以反映實際情況［8］。例如，在城市發(fā)展過程中，由于資源環(huán)境約束條件變化和政策等的調(diào)整，使得城市空間格局的演變規(guī)律在不同經(jīng)濟發(fā)展階段有所變化。利用CA進行模擬預(yù)測時，靜態(tài)的參數(shù)組合往往難以反映這種非線性地理過程。

動態(tài)調(diào)整模型參數(shù)的方法被認(rèn)為能夠較好地適應(yīng)非線性地理模擬，目前已有部分學(xué)者開展相關(guān)研究。文獻［3］提出利用肉眼判斷的方法來調(diào)整SLEUTH－CA模型參數(shù)值。此方法固定其他參數(shù)，反復(fù)調(diào)整某一參數(shù)尋求較優(yōu)模擬結(jié)果［3］。當(dāng)參數(shù)較多時，無法使用該方法，而且該方法有一定的主觀性。文獻［3］中又提出了計算搜索算法。利用計算機計算各個參數(shù)組合的模擬結(jié)果與實際結(jié)果的吻合度，確定最優(yōu)的參數(shù)組合［12］。但這種試驗需要高性能工作站運算幾百小時。文獻［6］提出案例推理（CBR）的CA模型，模型能夠根據(jù)新的案例動態(tài)更新轉(zhuǎn)換規(guī)則。但此方法忽略了模擬結(jié)果的誤差積累。文獻［12］提出基于數(shù)據(jù)同化的CA模型。模型考慮到模型誤差的積累，但是沒有動態(tài)調(diào)整轉(zhuǎn)換規(guī)則中的參數(shù)?？傮w而言，已有的研究難以在模擬過程中動態(tài)調(diào)整參數(shù)和控制誤差積累，這是影響CA模型在復(fù)雜非線性地理模擬中應(yīng)用的主要原因。

近年來興起的數(shù)據(jù)同化提供了一種能夠利用觀測數(shù)據(jù)自動調(diào)整并優(yōu)化模型狀態(tài)和參數(shù)的方法。它能夠融合不同來源、不同精度的直接或間接觀測數(shù)據(jù)，并且集成模型和各種觀測因子，自動調(diào)整模型模擬軌跡［13－14］。在數(shù)據(jù)同化的方法中，非線性濾波方法占有極其重要的地位，其中，集合卡爾曼濾波（EnKF）方法使用最為廣泛，在參數(shù)優(yōu)化估計中也有較多的應(yīng)用［15－16］。文獻［17］對一個二維非線性海風(fēng)模型進行參數(shù)優(yōu)化估計；文獻［18］提出了雙集合卡爾曼濾波方法優(yōu)化水文模型參數(shù)；文獻［19］在三維的云模型中對模型參數(shù)進行優(yōu)化；文獻［20］把聯(lián)合狀態(tài)矩陣的方法應(yīng)用到核污染擴散模型中研究模型參數(shù)變化。運用集合卡爾曼濾波方法能夠很好地在模擬過程中動態(tài)修正模型參數(shù)，并且能夠在初始參數(shù)不準(zhǔn)確的情況下快速捕捉參數(shù)特征使得被優(yōu)化的參數(shù)收斂于真實值，這為CA模型動態(tài)調(diào)整模型參數(shù)與狀態(tài)提供了可能。

本文提出了一種基于集合卡爾曼濾波（EnKF）動態(tài)優(yōu)化CA模型參數(shù)的方法。該方法先構(gòu)建由參數(shù)和狀態(tài)組成的聯(lián)合狀態(tài)矩陣，再利用觀測數(shù)據(jù)更新聯(lián)合矩陣，最后把更新后的狀態(tài)和參數(shù)重新輸入到模型中進行下一次模擬迭代。文中以廣東省東莞市為試驗區(qū)，模擬其城市擴張，分析了單參數(shù)優(yōu)化情況下參數(shù)的有效性，并討論了在多參數(shù)情況下參數(shù)收斂情況和城市模擬情景。

2 方 法

2.1 MCE－CA模型

在模擬復(fù)雜的城市演變過程中，元胞自動機以其“自下而上”、自組織性等優(yōu)點，被認(rèn)為是一種非常有效的模擬復(fù)雜城市系統(tǒng)的模型［3－4，11－21］。與數(shù)學(xué)模型相比，CA模型具有簡潔的形式并能很好地結(jié)合柵格GIS數(shù)據(jù)。通過定義合適的轉(zhuǎn)換規(guī)則，便可模擬出城市的發(fā)展情景［3－7］。研究表明，影響城市擴展的變量主要是空間距離變量，如離商業(yè)中心的最短距離、離道路的最短距離等，而這些變量對應(yīng)的參數(shù)值（權(quán)重）反映了對模型的“貢獻”度［4－7，21］。文獻［21］運用這些空間變量并提出了多準(zhǔn)則的判斷（MCE）來表達其他土地類型轉(zhuǎn)變?yōu)槌鞘杏玫氐母怕?/p>

式中，Pg是元胞（i，j）全局性開發(fā)概率；Si，j是元胞當(dāng)前狀態(tài)；z是描述元胞（i，j）開發(fā)特征的向量

式中，a、b為回歸模型的系數(shù)；x為距離變量。

全局性開發(fā)概率Pg中考慮到各種空間距離變量對元胞轉(zhuǎn)化的影響，在微觀層次上，鄰域影響是一個非常重要的因素。它通過一個3×3的核來計算元胞在空間上的相互影響，其計算為

式中，Ωtij是鄰域函數(shù)；condition是一個條件函數(shù)；如果Si，j是城市用地，則返回真（用1表示），否則返回假（用0表示）。

在實際中，還應(yīng)考慮單元約束條件。有些土地類型，如陡峭的山地等不適宜城市發(fā)展，而平緩的地區(qū)則很適宜城市發(fā)展，這里用連續(xù)的條件約束來表達

式中，con是土地適宜性評估值（0～1），用連續(xù)值可提高模型的科學(xué)性與模擬的真實性。

綜合考慮了全局發(fā)展概率，鄰域范圍和約束條件后，任意元胞在t時刻發(fā)展為城市用地的概率可由下式來表達

求出元胞發(fā)展概率后，還需要判斷此元胞是否發(fā)展為城市用地。文中采用閾值法，表達如下

式中，St＋1（ij）是元胞（i，j）在t＋1時刻的狀態(tài)；γ為隨機變量；β＝1／K，K為迭代數(shù)。若閾值pthreshold過大，城市形態(tài)過于集中；閾值過小，城市形態(tài)過于分散。根據(jù)文獻［4，21］提出的方案，并通過試驗，確定閾值為0.65。

MCE－CA模型是模擬城市發(fā)展的經(jīng)典模型。在模型模擬預(yù)測時，其假設(shè)條件是城市發(fā)展模式不隨時間變化。實際上，由于政治和經(jīng)濟等因素，城市發(fā)展的模式會改變。例如，某城市前期是鎮(zhèn)中心擴張發(fā)展為主；后期則以沿公路發(fā)展為主。在這種情況下，若模型用前期的轉(zhuǎn)換規(guī)則來模擬整個時期的城市擴張，模型誤差不僅不會消除，還會隨著模擬迭代的進行而傳遞下去。因此，考慮引入集合卡爾曼濾波算法到CA模型中，通過同化觀測數(shù)據(jù)動態(tài)更新轉(zhuǎn)換規(guī)則中的參數(shù)，調(diào)整模型運行軌跡，以減少誤差積累。

2.2 基于EnKF的數(shù)據(jù)同化算法

集合卡爾曼濾波（EnKF）是文獻［16］根據(jù)隨機動力預(yù)報理論而提出的一種順序同化算法。該算法首先用Monte Carlo方法生成初始集合，然后用集合預(yù)報的思想對初始集合進行預(yù)測得出預(yù)測集合，再把預(yù)測集合和觀測集合輸入到集合卡爾曼濾波更新方程中得到更新后預(yù)測集合即分析集合，最后用分析集合作初始集合預(yù)測到下一個觀測時刻，如此循環(huán)。其實現(xiàn)步驟如下：

（1）初始階段，模型狀態(tài)集合初始化。

式中，n為模型狀態(tài)變量的個數(shù)；N為集合的大小；下標(biāo)0為時間刻度。

（2）預(yù)測階段，將分析集合（或者初始集合）輸入到模型中，并預(yù)測下一個觀測時刻k＋1得到的預(yù)測狀態(tài)集合。

式中，M為模型算子，本文指元胞自動機模型；Xkf＋1是k＋1時刻的預(yù)測值；Xak是k時刻的分析值，若k＝0，Xak為初始集合X0；ωk是期望為0；方差為Wk的高斯白噪聲；Wk為模型的誤差方差矩陣。

（3）分析階段，融合預(yù)測狀態(tài)集合和觀測集合，并根據(jù)集合卡爾曼濾波更新公式，求得分析集合Xak＋1

式中，Yk＋1是k＋1時刻的狀態(tài)變量的觀測值；H為觀測算子；vk＋1是期望為0，方差為R的高斯白噪聲；N為集合個數(shù)。Kk＋1為k＋1時刻的卡爾曼增益矩陣，其計算如下

式中，R是觀測誤差協(xié)方差矩陣；Pfk＋1為狀態(tài)誤差協(xié)方差矩陣。

在應(yīng)用中，常常直接計算Pfk＋1HT與HPfk＋1HT來減少運算量。

（4）判斷是否到結(jié)束時刻，若未到，返回步驟（2），否則結(jié)束。

集合卡爾曼濾波運用集合預(yù)報的思想，把誤差統(tǒng)計量的預(yù)報隱含在擾動的集合變量中，避免了傳統(tǒng)卡爾曼濾波等需要專門對誤差協(xié)方差進行計算的問題，也解決了卡爾曼濾波應(yīng)用在非線性系統(tǒng)中的切線性近似的問題。集合卡爾曼濾波不僅簡化了矩陣計算，還能利用隨機樣本的預(yù)報值估計協(xié)方差得到較高的準(zhǔn)確度。

2.3 CA模型的參數(shù)優(yōu)化

數(shù)據(jù)同化方法能夠充分利用觀測值等信息提高模型模擬精度。引入數(shù)據(jù)同化方法到CA模型中，筆者進行了初步探索并提出了基于數(shù)據(jù)同化的CA模型［12］。在該模型的基礎(chǔ)上，進一步考慮參數(shù)組合對模擬結(jié)果的影響，故引入集合卡爾曼濾波的聯(lián)合狀態(tài)矩陣方法控制由參數(shù)引起的模型誤差。該方法主要包括如下幾個步驟，其流程如圖1。

首先，根據(jù)2.1節(jié)的描述構(gòu)建起CA模型并確定需要進行估計的參數(shù)。確定的參數(shù)包括各個空間距離變量的系數(shù)，如公式（2）所示。本文中包括4個，分別是brailways、bexpressways、btowns和broads。運用蒙特卡羅方法生成N（N為EnKF中的集合個數(shù)）組不同的參數(shù)組合。生成隨機參數(shù)組合包含兩個必要部分，參數(shù)初始估計值與參數(shù)隨機擾動噪聲。隨機參數(shù)集合表示如下

式中，bi為擾動后的參數(shù)值為參數(shù)初始估計值；ei為擾動噪聲，擾動噪聲的分布符合均值為0，標(biāo)準(zhǔn)差為真實值與初始估計值的差的絕對值。

圖1 動態(tài)優(yōu)化CA模型參數(shù)流程圖Fig.1 Flowchart of parameter optimization in CA model

其次，把生成的N組不同的參數(shù)值輸入到CA模型中進行城市模擬，可得出N組不同的模擬結(jié)果。此時CA模型可以用如下形式描述

式中，f為定義元胞從時刻t到t＋1時刻狀態(tài)的轉(zhuǎn)換函數(shù)；St是t時刻的元胞狀態(tài)；Ωt是鄰域函數(shù)；brailwaysi、bexpresswaysi、btownsi和broadsi指第i組參數(shù)組合中各個空間距離變量的系數(shù)值。

再次，整合N組參數(shù)的模擬結(jié)果到狀態(tài)矩陣中，以計算狀態(tài)誤差協(xié)方差。在計算誤差協(xié)方差前，考慮到CA模型模擬結(jié)果只有兩種形式（發(fā)展與不發(fā)展），不能直接參與計算。所以，本文考慮把研究區(qū)分成若干個規(guī)則的方格，每個方格包含m×m個像元（元胞），用城市發(fā)展密度參與數(shù)據(jù)同化運算。每個方格的城市發(fā)展密度為

獲得每個方格的城市發(fā)展密度后，便可用發(fā)展密度來計算誤差協(xié)方差。同理，觀測值也是用發(fā)展密度來計算。

最后，采用文獻［17］中的方法把參數(shù)等同于狀態(tài)，整合N組CA模型參數(shù)組合及其模擬結(jié)果的發(fā)展密度到同一個聯(lián)合狀態(tài)矩陣Xj里，即把參數(shù)brailwaysi、bexpresswaysi、btownsi和broadsi作為狀態(tài)增加到原有的狀態(tài)矩陣X中，構(gòu)成一個包含參數(shù)和狀態(tài)的聯(lián)合狀態(tài)矩陣Xj。再將此聯(lián)合狀態(tài)矩陣Xj代替X輸入集合卡爾曼濾波方法的更新方程（9）中。通過同化觀測數(shù)據(jù)可得到優(yōu)化后的參數(shù)值，并將其輸入CA模型，以進行下一時刻的模擬迭代，如此循環(huán)。

通過重復(fù)上述幾個步驟便可實現(xiàn)在模擬過程中動態(tài)調(diào)整模型參數(shù)，進而改變CA模型的轉(zhuǎn)換規(guī)則。此外，在集合卡爾曼濾波的分析階段，也可根據(jù)更新方程得出的最優(yōu)分析場并按照文獻［12］的方法對模擬結(jié)果進行修正。

3 試驗和分析

3.1 數(shù)據(jù)和模型參數(shù)

東莞市是珠江三角洲城市化發(fā)展中速度最快的城市之一。期間，該區(qū)域土地利用變化劇烈、城市擴張速度快，而且還經(jīng)歷多次產(chǎn)業(yè)布局調(diào)整，城市發(fā)展表現(xiàn)出復(fù)雜的非線性過程。選擇該區(qū)域具有典型的研究意義。

選用模擬起始數(shù)據(jù)為1993年東莞市遙感TM分類數(shù)據(jù)。城市發(fā)展的距離空間變量有到鐵路的距離、到高速公路的距離、到道路的距離和到鎮(zhèn)中心的距離（圖2）。運用SPSS軟件中的邏輯回歸方法對1993—1995年訓(xùn)練樣本進行運算得出的各個距離變量的權(quán)重（公式（2）中的bk）作為初始參數(shù)值。為了觀察參數(shù)調(diào)整后模型誤差情況以及方便對所提出的方法進行對比驗證，選擇觀測資料基于如下考慮：若引入多景遙感觀測資料來觀察參數(shù)收斂情況，往往會產(chǎn)生這種情況，縱然獲取的參數(shù)準(zhǔn)確，仍難模擬出“飛地”、“開發(fā)區(qū)”等人為干預(yù)的城市發(fā)展用地。而且此時模擬誤差主要是由CA根本規(guī)則——“鄰域發(fā)展規(guī)則”引起的，難以剝離參數(shù)誤差與鄰域規(guī)則誤差進行比較驗證。若使用較短模擬時間，期間的觀測資料又非常難獲取?？紤]到CA模型在模擬短時段城市發(fā)展有較高的精度（本文為87.7%），故文中用模擬的數(shù)據(jù)代替無法獲取的歷史觀測數(shù)據(jù)，同時，模擬的觀測數(shù)據(jù)也給參數(shù)收斂結(jié)果的驗證提供了理想的真值。獲取觀測數(shù)據(jù)的具體過程如下，首先用SPSS軟件中的邏輯回歸方法對1993—1997年訓(xùn)練樣本進行運算，得出的各個距離變量的權(quán)重并作為參數(shù)的真值（見表1）。然后把這些權(quán)重輸入CA模型中，進行模擬并得出模擬結(jié)果，最后把模擬結(jié)果轉(zhuǎn)換成發(fā)展密度的形式，并選擇如圖3所示的30個觀測點數(shù)據(jù)作為觀測數(shù)據(jù)。

圖2 空間距離變量Fig.2 Spatial variables of CA model

表1 參數(shù)初始值和初始方差Tab.1 Initial parameter values and variances

集合（ensemble size）過大，運算量會相當(dāng)巨大且對結(jié)果改善不明顯；集合過小，體現(xiàn)不出統(tǒng)計特征并且誤差較大［13，22］，本文選擇集合大小為30。由于傳統(tǒng)CA模型的模擬精度約為80%，故設(shè)定模型誤差為0.2。觀測誤差取決于觀測數(shù)據(jù)的精度（87.7%），故設(shè)定觀測誤差為0.13。

除了設(shè)置起始模擬數(shù)據(jù)、初始參數(shù)等相關(guān)參數(shù)外，還需注意生成擾動參數(shù)的方差和參數(shù)演化控制。擾動參數(shù)的方差大小對模型模擬結(jié)果的誤差統(tǒng)計有極大影響。文中用參數(shù)的標(biāo)準(zhǔn)差等于初始誤差，初始誤差等于初始估計值與真值之差的絕對值。在早期的濾波同化應(yīng)用中，同化時刻的參數(shù)演化常常是將隨機擾動加入原參數(shù)的集合里生成新的參數(shù)集合。此方法往往會由于參數(shù)本身的特性導(dǎo)致濾波發(fā)散和信息丟失［18］。文獻［23］提出了核平滑進行參數(shù)更新以避免濾波發(fā)散與信息丟失。其參數(shù)核平滑表示如下

為了檢驗聯(lián)合狀態(tài)矩陣方法在CA模擬中的適用性，進行了單參數(shù)優(yōu)化試驗。在真實的模型參數(shù)優(yōu)化中，往往是多參數(shù)同時優(yōu)化，所以進一步實現(xiàn)了多參數(shù)優(yōu)化，并分析了參數(shù)收斂效果與參數(shù)優(yōu)化前后的城市模擬情景。

圖3 觀察點空間分布Fig.3 Position of observation

3.2 單參數(shù)優(yōu)化

在單參數(shù)優(yōu)化過程中，僅僅被優(yōu)化的參數(shù)按公式（15）進行高斯擾動，其他參數(shù)保持真值不變。單參數(shù)試驗中，模型各個參數(shù)優(yōu)化結(jié)果如圖4和圖5所示。除了記錄參數(shù)優(yōu)化過程，圖中還給出了模擬結(jié)果的均方根誤差（RMS）。在圖例中，“未同化”與“同化后”分別表示未進行參數(shù)優(yōu)化與經(jīng)過參數(shù)優(yōu)化后的結(jié)果。雖然各個參數(shù)逼近真值的速度有所差別，但是都能夠在較少的同化次數(shù)后逼近真值，并在真值的較小范圍內(nèi)波動。各個參數(shù)在收斂于真實參數(shù)值的時間上顯示了一定的相似性，都是在10次同化操作前后收斂于各自的真實參數(shù)值。由于各個參數(shù)對觀測數(shù)據(jù)的敏感程度不同，所以收斂速度有一所差別。從RMS的變化情況可以注意到，同化后的誤差均比沒有進行同化的要低。參數(shù)“普通公路”與“鎮(zhèn)中心”的誤差積累約為未同化的1／3，而參數(shù)“高速公路”與“鐵路”誤差只有未同化的1／6。通過動態(tài)優(yōu)化參數(shù)和狀態(tài)，模型能夠較好地控制模型誤差積累。

圖4 優(yōu)化的參數(shù)值與RMS變化（1）Fig.4 Parameter value and RMS in parameter optimization（1）

圖5 估計的參數(shù)值與RMS變化（2）Fig.5 Parameter value and RMS in parameter optimization（2）

在單參數(shù)優(yōu)化試驗可以發(fā)現(xiàn)，本方法能夠在初始參數(shù)不準(zhǔn)確的情況下快速獲得參數(shù)特征，并使得被優(yōu)化的參數(shù)收斂于真實值。

3.3 多參數(shù)聯(lián)合優(yōu)化

為了完整地檢驗集合卡爾曼濾波在優(yōu)化模型參數(shù)的效果，分別進行了兩組雙參數(shù)和一組四參數(shù)同時優(yōu)化的試驗。在雙參數(shù)優(yōu)化中，被優(yōu)化的參數(shù)進行高斯擾動，其他參數(shù)保持真值。兩組雙參數(shù)優(yōu)化后的模擬精度見表2，參數(shù)收斂情況如圖6所示。結(jié)果表明，模擬精度有一定程度的提高，參數(shù)收斂情況與單參數(shù)優(yōu)化相似。經(jīng)過較少次的更新操作后，各個優(yōu)化的參數(shù)均能較快地逼近真值，并在相對較小的誤差范圍內(nèi)波動。同時，模型誤差RMS均有較大程度的降低。

圖6 雙參數(shù)聯(lián)合優(yōu)化的結(jié)果Fig.6 Parameter values and RMS for two spatial variables

當(dāng)把優(yōu)化的參數(shù)擴展到所有的參數(shù)集合時，各個參數(shù)都進行高斯擾動，而且其初始值均偏離真值。通過同化觀測資料，筆者得到的結(jié)果見圖7與表2。總體上，各個參數(shù)的收斂速度迅速，模型誤差較低，模擬精度有所改善。在優(yōu)化結(jié)果中，“高速公路”、“鐵路”與“鎮(zhèn)中心”3個參數(shù)均能較快地收斂，然而參數(shù)“普通公路”卻一直波動。這主要是由于兩個原因：①與其他參數(shù)相比，參數(shù)“普通公路”比其他值要大，這就意味著“貢獻”大，參數(shù)的敏感性強，觀測值對參數(shù)的變化影響較大，故在一定時間段內(nèi)出現(xiàn)振蕩，需要較多的更新次數(shù)來捕捉參數(shù)特征；②在進行數(shù)據(jù)同化前，將CA模擬結(jié)果轉(zhuǎn)換成發(fā)展密度，是用均值代替各個元胞狀態(tài)值，這在一定程度上降低了參數(shù)的敏感性。

對比模擬結(jié)果圖（圖8）可以發(fā)現(xiàn)，未進行參數(shù)優(yōu)化的模擬結(jié)果在市區(qū)中南部的大嶺山、市區(qū)東南部寮步以及研究區(qū)的東南部塘廈均有模擬過度的情況。這是由于轉(zhuǎn)換規(guī)則的誤差引起的，東莞市城市發(fā)展前期是圍繞城鎮(zhèn)中心進行擴張發(fā)展，后期沿公路進行發(fā)展。未進行優(yōu)化的CA模型因不能在模擬過程中改變轉(zhuǎn)換規(guī)則，會一直以圍繞城市中心發(fā)展的模式模擬城市演化過程，因此，在某些城鎮(zhèn)中心會產(chǎn)生模擬過度的情況。引入集合卡爾曼濾波后，模型能夠不斷修正轉(zhuǎn)換規(guī)則中的參數(shù)。修正后的參數(shù)能夠更好地反映真實情況與各個空間距離變量的關(guān)系，因此，參數(shù)優(yōu)化后的模擬結(jié)果基本上不存在過度發(fā)展的情況。模擬結(jié)果與真實參數(shù)的模擬結(jié)果更接近，也與實際情況很吻合。對比RMS與城市模擬結(jié)果，優(yōu)化后的模擬結(jié)果均有一定程度的改善。而且各個參數(shù)收斂效果良好，在經(jīng)過較少次的更新優(yōu)化操作后，均能夠在真值相對較小的誤差范圍內(nèi)波動。

表2 多參數(shù)優(yōu)化后的精度對比Tab.2 Comparison of precision and kappa coefficient for multiple parameters optimization

圖7 四參數(shù)聯(lián)合優(yōu)化的結(jié)果Fig.7 Parameter values and RMS for four spatial variables

隨著參與優(yōu)化的參數(shù)越多，與單參數(shù)優(yōu)化相比，部分參數(shù)需要的觀測數(shù)據(jù)會越多。在圖6和圖7中，某些時刻的參數(shù)值會突然偏離真值，這主要是由于觀測誤差引起的，一些觀測數(shù)據(jù)在某些時刻會由于觀測誤差而偏離真實值，從而導(dǎo)致同化結(jié)果也會有一定的誤差。雖然部分參數(shù)需要較多的同化次數(shù)來收斂于真值，但是從總體上，觀察數(shù)據(jù)的選擇和觀測誤差的估計是合理的，從而獲取到準(zhǔn)確的信息，較好地反映了參數(shù)自動優(yōu)化的整個過程。

圖8 四參數(shù)優(yōu)化的模擬結(jié)果Fig.8 Simulation results of four model parameters optimization

4 結(jié) 論

本文采用集合卡爾曼濾波方法更新CA模型參數(shù)與模型狀態(tài)，并以東莞市為試驗區(qū)進行模擬試驗，得到如下結(jié)論：

（1）本文方法通過構(gòu)建由模型參數(shù)和狀態(tài)組成的聯(lián)合狀態(tài)矩陣，并利用觀測數(shù)據(jù)更新此聯(lián)合狀態(tài)矩陣，得出更新后的參數(shù)與狀態(tài)，并將其重新輸入到模型中進行下一時刻的模擬。該方法有效地實現(xiàn)了模型參數(shù)和狀態(tài)的動態(tài)優(yōu)化，提高了模型在復(fù)雜非線性地理模擬中適應(yīng)能力。

（2）單參數(shù)和多參數(shù)聯(lián)合優(yōu)化試驗表明，該方法能夠很好地優(yōu)化CA模型參數(shù)值并快速收斂于真實參數(shù)值，而且還能夠有效地減少模型誤差。在多參數(shù)優(yōu)化試驗中，模擬結(jié)果的模擬精度與Kappa系數(shù)均有一定程度的提高。

（3）由于城市發(fā)展模式的改變，各個空間距離變量在模擬過程中對模擬結(jié)果的“貢獻”會變化，即，空間變量對CA模型轉(zhuǎn)換概率的影響會改變。同時空間變量對應(yīng)權(quán)重的不同會引起參數(shù)敏感度不同，其參數(shù)收斂速度會有所差別。

（4）空間距離變量對模擬結(jié)果有重要的影響，在如何選用合適方法以找尋最優(yōu)空間變量組合和發(fā)展密度網(wǎng)格大小對模型參數(shù)的收斂等方面需要作進一步研究。

［1］ WHITE R，ENGELEN G.Cellular Automata and Fractal Urban Form：A Cellular Modelling Approach to the Evolution of Urban Land Use Patterns ［J］.Environment and Planning A，1993，25（8）：1175－1199.

［2］ TOBLER W R.A Computer Movie Simulating Urban Growth in the Detroit Region［J］.Economic Geography，1970，46：234－240.

［3］ CLARKE K C，HOPPEN S，GAYDOS L.A Self－modifying Cellular Automata Model of Historical Urbanization in the San Francisco Bay Area［J］.Environment and Planning B，1997，24（2）：247－261.

［4］ LI X，YEH A G O.Modelling Sustainable Urban Development by the Integration of Constrained Cellular Automata and GIS［J］.International Journal of Geographical Information Science，2000，14（2）：131－152.

［5］ LI X，YEH A G O.Data Mining of Cellular Automata＇s Transition Rule［J］.International Journal of Geographical Information Science，2004，18（8）：723－744.

［6］ LI X，LIU X P.An Extended Cellular Automaton Using Case－based Reasoning for Simulating Urban Development in a Large Complex Region ［J］.International Journal of Geographical Information Science，2006，20（10）：1109－1136.

［7］ LIU Xiaoping，LI Xia.Fisher Discriminant and Automatically Getting Transition Rule of CA ［J］.Acta Geodaetica et Cartographica Sinica，2007，36（1）：112－118.（劉小平，黎夏.Fisher判別及自動獲取元胞自動機的轉(zhuǎn)換規(guī)則［J］.測繪學(xué)報，2007，36（1）：112－118.）

［8］ CLARKE K C，BRASS J A，RIGGAN P J.A Cellular Automata Model of Wildfire Propagation and Extinction［J］.Photogrammetric Engineering ＆ Remote Sensing，1994，60：1355－1367.

［9］ COUCLELIS H.Of Mice and Men：What Rodent Populations Can Teach Us About Complex Spatial Dynamics［J］.Environment and Planning A，1988，20：99－109.

［10］ CHEN Jianping，DING Huoping，WANG Gongwen，et al.Desertification Evolution Modeling through the Integration of GIS and Cellular Automata［J］.Journal of Remote Sensing，2004，8（3）：254－260.（陳建平，丁火平，王功文，等.基于GIS和元胞自動機的荒漠化預(yù)測模型［J］.遙感學(xué)報，2004，8（3）：254－260.）

［11］ CLARKE K C，GAYDOS L.Loose－coupling a Cellular Automata Model and GIS：Long－term Urban Growth Prediction for San Francisco and Washington／Baltimore［J］.International Journal of Geographical Information Science，1998，12（7）：699－714.

［12］ ZHANG Yihan，LI Xia，LIU Xiaoping，et al.The CA Model Based on Data Assimilation［J］.Journal of Remote Sensing，2011，15（3）：475－491.（張亦漢，黎夏，劉小平，等.基于數(shù)據(jù)同化的元胞自動機［J］.遙感學(xué)報，2011，15（3）：475－491.）

［13］ REICHLE R H，MCLAUGHLIN D B，ENTEKHABI D.Hydrologic Data Assimilation with the Ensemble Kalman Filter［J］.Monthly Weather Review，2002，130（1）：103－114.

［14］ LI Xin，HUANG Chunlin，CHE Tao，et al.Development of a Chinese Land Data Assimilation System：Its Progress and Prospects［J］.Progress in Natural Science，2007，17（2）：163－173.（李新，黃春林，車濤，等.中國陸面數(shù)據(jù)同化系統(tǒng)研究的進展與前瞻［J］.自然科學(xué)進展，2007，17（2）：163－173.）

［15］ LI Xin，BAI Yulong.A Bayesian Filter Framework for Sequential Data Assimilation ［J］.Advance in Earth Science，2010，25：515－522.（李新，擺玉龍.順序數(shù)據(jù)同化的Bayes濾波框架［J］.地球科學(xué)進展，2010，25：515－522.）

［16］ EVENSEN G.Sequential Data Assimilation with a Nonlinear Quasigeostrophic Model Using Monte－Carlo Methods to Forecast Error Statistics ［J］.Journal of Geophysical Research，1994，99（C5）：10143－10162.

［17］ AKSOY A，Mesoscale Ensemble－based Data Assimilation and Parameter Estimation ［D］.Texas：Texas A＆M University，2005.

［18］ MORADKHANI H，SOROOSHIAN S，GUPTA H V，et al. Dual State－parameter Estimation of Hydrological Models Using Ensemble Kalman Filter［J］.Advances in Water Resources，2005，28（2）：135－147.

［19］ TONG M J，XUE M.Simultaneous Estimation of Microphysical Parameters and Atmospheric State with Simulated Radar Data and Ensemble Square－root Kalman Filter.Part II：Parameter Estimation Experiments ［J］.Monthly Weather Review，2008（136）：1649－1668.

［20］ ZHENG D Q，LEUNG J K C，LEE B Y.Online Update of Model State and Parameters of a Monte Carlo Atmospheric Dispersion Model by Using Ensemble Kalman Filter［J］.Atmospheric Environment，2009，43（12）：2005－2011.

［21］ WU F，WEBSTER C J.Simulation of Land Development through the Integration of Cellular Automata and Multicriteria Evaluation ［J］.Environment and Planning B：Planning and Design，1998，25（1）：103－126.

［22］ CROW W T，WOOD E F.The Assimilation of Remotely Sensed Soil Brightness Temperature Imagery into a Land Surface Model Using Ensemble Kalman Filtering：A Case Study Based on ES－TAR Measurements during SGP97［J］.Advance Water Resources，2003，26：137－149.

［23］ LIU Fang.Bayesian Time Series：Analysis Methods Using Simulation Based Computations［D］.Durham：Duke University，2000.