馮帆王自發(fā)唐曉

1中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心/超級(jí)計(jì)算中心，北京1001902中國(guó)科學(xué)院大氣物理研究所大氣邊界層物理和大氣化學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京1000293中國(guó)科學(xué)院大學(xué)，北京100049

?

一個(gè)基于打靶法的大氣污染源反演自適應(yīng)算法

馮帆1, 2, 3王自發(fā)2唐曉2

1中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心/超級(jí)計(jì)算中心，北京100190
2中國(guó)科學(xué)院大氣物理研究所大氣邊界層物理和大氣化學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京100029
3中國(guó)科學(xué)院大學(xué)，北京100049

摘 要污染源反演對(duì)大氣污染預(yù)報(bào)及控制有重要意義。目前普遍采用的源反演統(tǒng)計(jì)方法存在對(duì)觀測(cè)誤差、源清單先驗(yàn)估計(jì)誤差敏感等弱點(diǎn)?；诖虬蟹ㄋ枷氲母鞣N算法以其精度高、程序簡(jiǎn)單、實(shí)用性強(qiáng)的特點(diǎn)被廣泛應(yīng)用于系統(tǒng)控制領(lǐng)域。本文提出的基于打靶法思想的大氣污染源反演自適應(yīng)算法在精度高、算法簡(jiǎn)明的基礎(chǔ)上彌補(bǔ)了統(tǒng)計(jì)方法的不足：能處理源清單中的大誤差、初值大誤差、觀測(cè)值在個(gè)別時(shí)間點(diǎn)的大誤差；無(wú)需先驗(yàn)分布假設(shè)及誤差估計(jì)。本文還以簡(jiǎn)單模型的理想試驗(yàn)為例，展示了該自適應(yīng)算法的計(jì)算效果。

關(guān)鍵詞打靶法 污染源 反演 自適應(yīng) 系統(tǒng)控制

馮帆，王自發(fā)，唐曉. 2016. 一個(gè)基于打靶法的大氣污染源反演自適應(yīng)算法 [J]. 大氣科學(xué), 40 (4): 719?729. Feng Fan, Wang Zifa, Tang Xiao. 2016. Development of an adaptive algorithm based on the shooting method and its application in the problem of estimating air pollutant emissions [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 40 (4): 719?729, doi:10.3878/j.issn.1006-9895.1505.15113.

1 引言

近年來(lái)，大氣污染模式的研究得到了很大發(fā)展。污染排放源信息對(duì)模式至關(guān)重要。然而，建立污染源清單是一個(gè)具有較大不確定性的復(fù)雜過(guò)程，涉及對(duì)經(jīng)濟(jì)發(fā)展?fàn)顩r、人口密度、能源結(jié)構(gòu)、工廠分布、交通流量、車型比例等多方面的調(diào)查研究。因此，根據(jù)模式的模擬值和監(jiān)測(cè)站的污染物濃度觀測(cè)值來(lái)反演污染源信息對(duì)大氣污染預(yù)報(bào)及控制均有重要意義。

大氣污染預(yù)報(bào)是已知原因求結(jié)果，屬于正問(wèn)題，即給定了微分方程和初、邊值條件，求滿足給定條件的解；而大氣污染控制是由結(jié)果求原因，屬于反問(wèn)題，即排放源等參數(shù)作為未知量出現(xiàn)在微分方程中，我們需利用觀測(cè)值等附加信息反求未知參數(shù)（劉峰等，2003；楊一帆和張凱山，2013）。多種方法已被應(yīng)用于反演污染源，比如逆向軌跡法（蔡旭暉等，2002）、遺傳算法（陳軍明等，2002）、伴隨方程法（劉峰等，2003；郭少冬等，2010）、四維變分法（Mendoza-Dominguez and Russell, 2001）、集合卡爾曼濾波法和集合卡爾曼平滑法（朱江和汪萍，2006；Tang et al., 2013）以及蒙特卡洛空間反演法（楊一帆和張凱山，2013）等。

目前，源反演普遍采用統(tǒng)計(jì)類方法，但其存在對(duì)觀測(cè)誤差、源清單先驗(yàn)估計(jì)誤差敏感等弱點(diǎn)。而基于打靶法思想的算法恰好能彌補(bǔ)統(tǒng)計(jì)法的不足。打靶法思想的本質(zhì)是將微分方程兩時(shí)刻點(diǎn)間的演化問(wèn)題轉(zhuǎn)化為控制問(wèn)題來(lái)求解，通過(guò)調(diào)整控制參量使觀測(cè)時(shí)刻點(diǎn)的模擬值接近觀測(cè)值。在源反演問(wèn)題中，這個(gè)控制參量就是源排放速率?；诖虬蟹ㄋ枷氲母鞣N算法以其精度高、程序簡(jiǎn)單、實(shí)用性強(qiáng)的特點(diǎn)被廣泛應(yīng)用于系統(tǒng)控制領(lǐng)域（向明華和劉云青，1990；屈香菊，1992；范慕輝和石鐵君，1994；凌復(fù)華和李繼躍，1988；Ning et al., 2004；吳文海等，2005；張洪波等，2008；齊笳羽等，2011）。

本文針對(duì)大氣污染源反演問(wèn)題提出一個(gè)基于打靶法思想的自適應(yīng)算法。該算法在精度高、算法簡(jiǎn)明的基礎(chǔ)上還彌補(bǔ)了許多統(tǒng)計(jì)方法的不足：能處理源清單中的大誤差、初值大誤差、觀測(cè)值在個(gè)別時(shí)間點(diǎn)的大誤差；無(wú)需先驗(yàn)分布假設(shè)及誤差估計(jì)。本文還以簡(jiǎn)單模型的理想試驗(yàn)為例，給出了該自適應(yīng)算法的計(jì)算效果。

2 源反演自適應(yīng)算法設(shè)計(jì)

從污染物濃度觀測(cè)值估計(jì)污染物排放速率可作為系統(tǒng)控制問(wèn)題進(jìn)行考慮。其中，污染物排放速率為控制變量。關(guān)于打靶法思想介紹，見附錄A。基于打靶法的思路，我們可以將觀測(cè)值作為靶子，設(shè)法尋找使模擬值與觀測(cè)值的距離小于容忍度的排放速率。

該算法設(shè)計(jì)的困難主要來(lái)自以下兩方面：第一，直觀上，當(dāng)模擬值大于觀測(cè)值時(shí)，我們需減小排放速率；當(dāng)模擬值小于觀測(cè)值時(shí)，我們需增大排放速率。那么如何增減能使射擊次數(shù)盡量小且程序簡(jiǎn)明；第二，除了源清單的誤差，實(shí)際中還不可避免地會(huì)遇到初值誤差、觀測(cè)誤差。當(dāng)這些誤差較大時(shí)，算法能否及時(shí)糾錯(cuò)，使模擬值盡快主動(dòng)回到正常位置（而不是被動(dòng)地等待長(zhǎng)時(shí)間模式預(yù)熱，即Spin Up）。增強(qiáng)算法的健壯性，使算法具有糾錯(cuò)功能的必要性說(shuō)明見附錄B。

因此，采用何種打靶策略能夠立刻對(duì)上述可能出現(xiàn)的大誤差進(jìn)行糾錯(cuò)并且較為快速、準(zhǔn)確地進(jìn)行源反演是本算法設(shè)計(jì)要解決的核心問(wèn)題。這里，進(jìn)行源反演可利用的唯一信息是誤差（Err）＝模擬值－觀測(cè)值。如何利用誤差 Err對(duì)排放速率進(jìn)行高效校正是十分重要的。判定何時(shí)需要對(duì)大誤差及時(shí)糾錯(cuò)以及如何糾錯(cuò)也是要解決的關(guān)鍵問(wèn)題。

2.1 源反演的打靶問(wèn)題描述

若不考慮模式誤差，大氣污染物濃度在某一空間點(diǎn)上的抽象模型可以表示為

其中，C是污染物濃度的向量（其分量表示不同污染物種的濃度）；是污染源排放速率的向量；算子F表示污染物的輸送、化學(xué)過(guò)程。

假定在t1, ..., tk時(shí)刻有污染物的濃度觀測(cè)向量Y1, ..., Yk，則在時(shí)間段［ti, ti+1］i=0, 1, …, k?1的污染物濃度演化過(guò)程可視為打靶問(wèn)題：在該時(shí)間段上控制Q，使得

其中，Tol為設(shè)定的允許誤差界限。于是，源反演問(wèn)題轉(zhuǎn)化成了分時(shí)間段的打靶問(wèn)題。Q 在每個(gè)時(shí)間段上為常值向量。

2.2 源反演的打靶算法設(shè)計(jì)

基于上面的敘述和分析，我們?cè)诒竟?jié)給出各時(shí)間段［ti, ti+1］上的自適應(yīng)算法流程。

在每個(gè)時(shí)間段［ti, ti+1］上，算法分為三大步：初始化、若干次射擊和保存本時(shí)間段的結(jié)果。若干次射擊又分成4小步，其中步驟（c）（用的當(dāng)前值在［ti, ti+1］上解初值問(wèn)題）即射擊。

在每次射擊后更新s維誤差向量Err＝C (ti+1;（其中C只提取與Y相同的物種分量），為下一次射擊前校正向量作準(zhǔn)備[見2.2.1節(jié)中步驟（d.1）]。

注意到只有通過(guò)有效的觀測(cè)信息，我們才有可能得到有效的源反演值。當(dāng)初值或觀測(cè)值在個(gè)別時(shí)間點(diǎn)存在大誤差時(shí)，反演信息必然無(wú)效。由于觀測(cè)值在若干個(gè)時(shí)間點(diǎn)可能存在大誤差，這使得系統(tǒng)在自然狀態(tài)下會(huì)經(jīng)歷若干個(gè)Spin Up，導(dǎo)致反演信息在Spin Up期間均無(wú)效。因此，主動(dòng)采取措施自動(dòng)糾錯(cuò)、增強(qiáng)算法的健壯性、使污染物濃度和源反演結(jié)果立刻回到正常位置是十分必要的（糾錯(cuò)的必要性圖例說(shuō)明見附錄B）。

當(dāng)初值或觀測(cè)值在個(gè)別時(shí)間點(diǎn)存在大誤差時(shí)，我們的算法恰恰能自動(dòng)糾錯(cuò)，使污染物濃度和源反演結(jié)果迅速回到正常位置。這就大大增強(qiáng)了算法的健壯性，從觀測(cè)值中盡可能多地反演出有效源排放速率[見2.2.1節(jié)中步驟（d.2）]。實(shí)際上，為了得到最終的有效信息，我們將對(duì)結(jié)果進(jìn)行后處理：去掉糾錯(cuò)階段的源排放速率（見本文數(shù)值試驗(yàn)中圖1e–f’）。

2.2.1 算法

（2）當(dāng)s維向量|Err|存在分量大于允許誤差界限Tol且ShootingCount小于MaxShooting（最大射擊次數(shù)）時(shí)，步驟分為：（a）ShootingCount＝ ShootingCount＋1；（b）用設(shè)定的經(jīng)驗(yàn)值初始化的當(dāng)前值在［ti, ti+1］上解初值問(wèn)題，得到各物種濃度計(jì)算方法如下：

For j = 1到s

（d.2）用Err(j) 更新Q (j)：[公式（3）為核心，其他均為了增強(qiáng)算法健壯性、為了糾錯(cuò)而設(shè)計(jì)]為正：

否則

（3）保存本時(shí)間段［ti, ti+1］的計(jì)算結(jié)果，為下一時(shí)間段作準(zhǔn)備；保存經(jīng)過(guò)若干次射擊最終接受的和用給下一時(shí)間段賦初值。

2.2.2 算法須注意問(wèn)題

（1）源清單中各物種的排放速率為非負(fù)值。在運(yùn)行此算法之前，我們需要對(duì)源清單進(jìn)行預(yù)處理：如果某物種有觀測(cè)值且該物種在源清單中某時(shí)刻的排放速率為零，則用一個(gè)小的正數(shù)代替零。這樣才能使步驟（d.2）的更新公式起作用。

圖1 試驗(yàn)1的結(jié)果：（a）濃度C1；（b）濃度C2；（b’）為（b）C2的局部放大；（c）濃度C3；（d）濃度C4；（e）Q1源排放速率；（e’）為（e）Q1去掉糾錯(cuò)階段；（f）Q2源排放速率；（f’）為（f）Q2去掉糾錯(cuò)階段Fig. 1 Results of Experiment 1: (a) Concentration C1, (b) concentration C2, (b’) partial enlargement of (b), (c) concentration of C3, (d) concentration of C4, (e) emissions rate of Q1, (e’) panel (e) without error correction phase, (f) emissions rate of Q2, (f’) panel (f) without error correction phase

圖1 （續(xù)）Fig. 1 (Continued)

（2）在試驗(yàn)中，我們對(duì)各物種使用了相同的允許誤差界限 Tol。實(shí)際上，為了進(jìn)一步減小射擊次數(shù)，不同的物種可以使用不同的 Tol：值域范圍大的物種取值相對(duì)大，值域范圍小的物種取值相對(duì)小。Tol 不宜過(guò)小。

（3）2.2.1節(jié)中步驟（d.2）的 CheckCritical 為設(shè)定的較小正值。當(dāng)源排放速率已經(jīng)非常小時(shí)，若仍有Err (j)＞Tol，那么我們決定使用負(fù)值糾錯(cuò)。實(shí)際上，這負(fù)值在后處理中作為糾錯(cuò)階段的源排放速率將被去掉。（見圖1e–f’）

3 簡(jiǎn)單模型的數(shù)值試驗(yàn)

在本節(jié)中，我們采用下面描述空氣質(zhì)量的簡(jiǎn)單無(wú)量綱模型進(jìn)行試驗(yàn)：

其中，Ci是第i個(gè)物種的濃度，是隨時(shí)間變化的污染源排放速率。試驗(yàn)在t?[0,12]上進(jìn)行，共3個(gè)試驗(yàn)，具體設(shè)置見表1。在試驗(yàn)中，我們假設(shè)真實(shí)初始濃度在試驗(yàn)1、試驗(yàn)2中，假設(shè)真實(shí)排放速率在試驗(yàn)3中，假設(shè)真實(shí)排放速率在我們的試驗(yàn)中不考慮模式誤差。試驗(yàn)各時(shí)間段的射擊次數(shù)見表2。

在本文的數(shù)值試驗(yàn)中，經(jīng)驗(yàn)參量 Qcoeff第一個(gè)分量對(duì)應(yīng)Q1，第二個(gè)分量對(duì)應(yīng)Q2。源清單的排放速率均為

試驗(yàn)1的目的是檢驗(yàn)本算法對(duì)污染物濃度的反演效果和污染源反演能力。從試驗(yàn)1的結(jié)果我們看到：本算法可使用較少射擊次數(shù)反演源（表2和表3），并且能處理源清單中的大誤差、初值大誤差、觀測(cè)值在個(gè)別時(shí)間點(diǎn)的大誤差（圖1；圖1b’在時(shí)刻6的觀測(cè)值存在大誤差）。反演后，沒(méi)有用觀測(cè)值反演的物種濃度也有可能得到改善（圖1d）。在源排放速率變化很小的情況下，我們從反演結(jié)果仍能夠看到源變化趨勢(shì)（圖1f’）。另外，試驗(yàn)1與附錄B展示了負(fù)值糾錯(cuò)功能的必要性。

試驗(yàn)2除了觀測(cè)時(shí)間間隔不同，其他均與試驗(yàn)1相同（表1）。對(duì)比試驗(yàn)1與試驗(yàn)2，其他設(shè)置不變，只在時(shí)間上加密觀測(cè)數(shù)據(jù)不一定得到更準(zhǔn)確的反演源（圖2），而且射擊次數(shù)可能會(huì)明顯增加（表2）。

表1 試驗(yàn)列表Table 1 The detailed setup of the experiments

表2 射擊次數(shù)列表Table 2 Shooting times over each time span

表3 試驗(yàn)1的源反演過(guò)程Table 3 The shooting process of emission in Experiment 1

圖2 試驗(yàn)2的結(jié)果：（a）濃度C1；（b）濃度C2；（b’）為（b）C1的局部放大；（c）Q1源排放速率；（c’）為（c）Q1去掉糾錯(cuò)階段；（d）Q2源排放速率；（d’）為（d）Q2去掉糾錯(cuò)階段Fig. 2 Results of Experiment 2: (a) Concentration C1, (b) concentration of C2, (b’) partial enlargement of (b), (c) emissions rate of Q1, (c’) panel (c) without error correction phase, (d) emissions rate of Q2, (d’) panel (d) without error correction phase

試驗(yàn)3給出了在初值、觀測(cè)值誤差不太大情況下的算法效果。反演后的源較好地與真值吻合（圖3）。

4 小結(jié)和討論

污染源反演問(wèn)題常作為統(tǒng)計(jì)估計(jì)問(wèn)題處理。比如，四維變分法、集合卡爾曼濾波是人們常用于反演的統(tǒng)計(jì)方法。四維變分法本質(zhì)上是使目標(biāo)函數(shù)極小化的泛函問(wèn)題。集合卡爾曼濾波是一種用蒙特卡羅短期集合預(yù)報(bào)來(lái)估計(jì)預(yù)報(bào)誤差協(xié)方差的四維變分方法。集合卡爾曼濾波無(wú)需模式反向積分，無(wú)需預(yù)報(bào)模式的切線性模式和伴隨模式，因此與四維變分法相比，更容易實(shí)現(xiàn)、可移植性強(qiáng)。然而，四維變分、集合卡爾曼濾波等統(tǒng)計(jì)方法存在下面幾個(gè)弱點(diǎn)：先驗(yàn)分布的高斯假設(shè)；線性的分析方程在模式為強(qiáng)非線性時(shí)，對(duì)結(jié)果造成不利影響；對(duì)觀測(cè)誤差較敏感；污染源排放的先驗(yàn)誤差估計(jì)對(duì)反演結(jié)果有較大影響（朱江等，2006）。

本文提出的基于打靶法思想的大氣污染源反演自適應(yīng)算法本質(zhì)上是系統(tǒng)控制問(wèn)題，與統(tǒng)計(jì)量無(wú)關(guān)，且計(jì)算過(guò)程中用到的主要信息為觀測(cè)值。只要觀測(cè)值本身合理、有效，就可以反演出可信的源排放速率。如果我們具備可靠的觀測(cè)誤差信息，則可利用該信息先對(duì)觀測(cè)數(shù)據(jù)作前處理校正，然后再應(yīng)用本文設(shè)計(jì)的打靶算法進(jìn)行源反演。

該基于打靶法思想的大氣污染源反演算法彌補(bǔ)了四維變分、集合卡爾曼濾波等統(tǒng)計(jì)方法的不足，具有如下特點(diǎn)：（1）能處理源清單中的大誤差、初值大誤差、觀測(cè)值在個(gè)別時(shí)間點(diǎn)的大誤差；（2）無(wú)需先驗(yàn)分布假設(shè)及誤差估計(jì)。在實(shí)際中，對(duì)污染源進(jìn)行反演是一個(gè)很復(fù)雜的問(wèn)題。污染源與觀測(cè)站點(diǎn)的空間分布不均勻性、可觀測(cè)的物種少于需要反演的污染源物種等局限性會(huì)使反演結(jié)果不唯一。在今后的工作中，我們會(huì)考慮通過(guò)有效地增加約束條件或信息來(lái)處理這一問(wèn)題。

圖2 （續(xù)）Fig. 2 (Continued)

氣象場(chǎng)誤差、模式分辨率誤差、化學(xué)機(jī)制誤差等因素均可造成模擬與觀測(cè)不一致。如果忽略這些模式誤差、僅通過(guò)調(diào)整排放源來(lái)減小模擬和觀測(cè)之間的不一致，那么源反演結(jié)果將會(huì)出現(xiàn)偏差。集合卡爾曼濾波/平滑考慮了各種誤差對(duì)源反演的影響（Tang et al., 2013）。如何將集合卡爾曼濾波/平滑與本文提出的基于打靶法思想的自適應(yīng)算法結(jié)合起來(lái)解決源反演問(wèn)題也是我們今后將要探討的。

圖3 試驗(yàn)3的結(jié)果：（a）濃度C1；(b) 濃度C2；（c）Q1源排放速率；（d）Q2源排放速率Fig. 3 Results of Experiment 2: (a) Concentration C1, (b) concentration C2, (c) emissions rate of Q1, (d) emissions rate of Q2

致謝 感謝中國(guó)科學(xué)院大氣物理研究所朱江研究員及審稿專家們對(duì)本文提出的寶貴意見。

參考文獻(xiàn)（References）

蔡旭暉, 邵敏, 蘇芳. 2002. 甲烷排放源逆向軌跡反演模式研究 [J]. 環(huán)境科學(xué), 23 (5): 19–24. Cai Xuhui, Shao Min, Su Fang. 2002. A backward trajectory inversion model for methane emission over Beijing area [J]. Environ. Sci. (in Chinese), 23 (5): 19–24, doi:10.3321/j.issn: 0250-3301.2002.05.004.

陳軍明, 徐大海, 朱蓉. 2002. 遺傳算法在點(diǎn)源擴(kuò)散濃度反演排放源強(qiáng)中的應(yīng)用 [J]. 氣象, 28 (9): 12–16. Chen Junming, Xu Dahai, Zhu Rong. 2002. Application of genetic algorithms to point-source inversion [J]. Meteor. Mon. (in Chinese), 28 (9): 12–16, doi:10.3969/j.issn.1000-0526.2002.09.003.

范慕輝, 石鐵君. 1994. 打靶法在力法解超靜定梁中的應(yīng)用 [J]. 河北工學(xué)院學(xué)報(bào), 23 (4): 32–37. Fan Muhui, Shi Tiejun. 1994. The shooting method’s application in solving statically indeterminate beams [J]. J. Hebei Inst. Technol. (in Chinese), 23 (4): 32–37.

郭少冬, 楊銳, 蘇國(guó)鋒, 等. 2010. 基于伴隨方程和 MCMC方法的室內(nèi)污染源反演模型研究 [J]. 應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào), 18 (4): 695–704. Guo Shaodong, Yang Rui, Su Guofeng, et al. 2010. Investigation of the inversion modeling for indoor contaminant source based on the adjoint equation and MCMC method [J]. J. Basic Sci. Engin. (in Chinese), 18 (4):695–704, doi:10.3969/j.issn.1005-0930.2010.04.017.

凌復(fù)華, 李繼躍. 1988. 打靶法在常微分方程邊界層型奇異攝動(dòng)問(wèn)題中的應(yīng)用 [J]. 應(yīng)用數(shù)學(xué)和力學(xué), 9 (7): 609–615. Ling Fuhua, Li Jiyue. 1988. Shooting method in singular perturbation problem of ordinary differential equations with boundary layers [J]. Applied Mathematics and Mechanics (in Chinese), 9 (7): 609–615.

劉峰, 胡非, 王鍔一. 2003. 用伴隨方程研究空氣污染的優(yōu)化控制 [J].環(huán)境科學(xué)學(xué)報(bào), 23 (4): 472–475. Liu Feng, Hu Fei, Wang Eyi. 2003. Study on optimal control of air pollution using adjoint equation [J]. Acta Scientiae Circumstantiae (in Chinese), 23 (4): 472–475, doi:10.3321/j. issn:0253–2468.2003.04.012.

Mendoza-Dominguez A, Russell A G. 2001. Estimation of emission adjustments from the application of four-dimensional data assimilation to photochemical air quality modeling [J]. Atmos. Environ., 35 (16): 2879–2894, doi:10.1016/S1352-2310(01)00084-X.

Ning Jiping, Han Qun, Chen Zhiqiang, et al. 2004. A powerful simple shooting method for designing multi-pumped fibre Raman amplifiers [J]. Chin. Phys. Lett., 21 (11): 2184–2187, doi:10.1088/0256-307X/21/11/ 030.

齊笳羽, 祝金川, 李成仁, 等. 2011. 打靶法實(shí)現(xiàn)Lorenz-Haken激光混沌的周期和穩(wěn)態(tài)控制 [J]. 半導(dǎo)體光電, 32 (1): 132–134. Qi Jiayu, Zhu Jinchuan, Li Chengren, et al. 2011. Realization of periods and steadystate control for Lorenz-Haken laser chaos with shooting method [J]. Semiconductor Optoelectronics (in Chinese), 32 (1): 132–134.

屈香菊. 1992. 直接多重打靶法在軌跡優(yōu)化方面的應(yīng)用 [J]. 飛行力學(xué), 10 (1): 13–21, 39. Qu Xiangju. 1992. The application of the multiple shooting algorithm to trajectory optimization [J]. Flight Dynamics (in Chinese), 10 (1): 13–21, 39.

Tang Xiao, Zhu Jiang, Wang Zifa, et al. 2013. Inversion of CO emissions over Beijing and its surrounding areas with ensemble Kalman filter [J]. Atmos. Environ., 81: 676–686, doi:10.1016/j.atmosenv.2013.08.051.

吳文海, 周志剛, 韓強(qiáng), 等. 2005. 基于打靶法的戰(zhàn)機(jī)攻擊導(dǎo)引方法研究[J]. 飛機(jī)設(shè)計(jì), (1): 47–51. Wu Wenhai, Zhou Zhigang, Han Qiang, et al. 2005. Study on attack trajectory of pursuiting fighter based on shooting method [J]. Aircraft Design (in Chinese), (1): 47–51, 10.3969/ j.issn.1673-4599.2005.01.011.

向明華, 劉云青. 1990. 邊值打靶法的改進(jìn)與應(yīng)用 [J]. 哈爾濱船舶工程學(xué)院學(xué)報(bào), 11 (4): 449–456. Xiang Minghua, Liu Yunqing. 1990. Modification and application of the multiple shooting method [J]. Journal of Harbin Shipbuilding Engineering Institute (in Chinese), 11 (4): 449–456.

楊一帆, 張凱山. 2013. 突發(fā)型大氣污染源位置識(shí)別反演問(wèn)題的數(shù)值模擬 [J]. 環(huán)境科學(xué)學(xué)報(bào), 33 (9): 2388–2394. Yang Yifan, Zhang Kaishan. 2013. Numerical simulation on source identification of accidentally occurring air pollution [J]. Acta Scientiae Circumstantiae (in Chinese), 33 (9): 2388–2394.

張洪波, 鄭偉, 湯國(guó)建. 2008. 采用打靶法設(shè)計(jì)考慮地球扁率的機(jī)動(dòng)軌道 [J]. 宇航學(xué)報(bào), 29 (4): 1177–1181. Zhang Hongbo, Zheng Wei, Tang Guojian. 2008. Maneuver trajectory design with J2correction based on shooting procedure [J]. Journal of Astronautics (in Chinese), 29 (4): 1177–1181, doi:10.3873/j.issn.1000-1328.2008.04.017.

朱江, 汪萍. 2006. 集合卡爾曼平滑和集合卡爾曼濾波在污染源反演中的應(yīng)用 [J]. 大氣科學(xué), 30 (5): 871–882. Zhu Jiang, Wang Ping. 2006. Ensemble Kalman smoother and ensemble Kalman filter approaches to the joint air quality state and emission estimation problem [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 30 (5): 871–882, doi:10.3878/j.issn.1006-9895.2006.05.16.

附錄A 打靶法基本思想簡(jiǎn)介

在科學(xué)與工程計(jì)算中，會(huì)經(jīng)常遇到邊值問(wèn)題。打靶法源于將邊值問(wèn)題化為初值問(wèn)題來(lái)解。對(duì)于二階常微分方程邊值問(wèn)題：

可以化為初值問(wèn)題來(lái)求解：

我們?cè)O(shè)法確定 A的值，使得初值問(wèn)題［公式（A-2）］的解(;)y t A在tb=的值(;)y b A滿足

其中，Tol為我們?cè)试S的誤差界限。

于是，我們得到了邊值問(wèn)題［公式（A-1）］的近似解(;)y t A。

顯然，圖A-1打靶法思想是系統(tǒng)控制問(wèn)題：A為控制變量，b為t = b時(shí)的系統(tǒng)觀測(cè)值，通過(guò)調(diào)節(jié)A使模擬值(;)y b A接近觀測(cè)值b。

圖A-1 打靶法示意圖：將邊值問(wèn)題化為初值問(wèn)題來(lái)解Fig. A-1 Sketch of the shooting method: Transforming a boundary value problem into an initial value problem

附錄B 算法中糾錯(cuò)功能的必要性——增強(qiáng)健壯性

去掉負(fù)值糾錯(cuò)功能，如圖所示，C1的濃度值從時(shí)刻10開始才逐漸Spin Up到正常位置（圖B-1a）；C2的濃度值需要更長(zhǎng)的時(shí)間才能 Spin Up到正常位置（圖B-1b）；Q1源排放速率從時(shí)刻11開始才逐漸提供有效信息（圖B-1c）。Q2源排放速率需要更長(zhǎng)的時(shí)間才能提供有效信息（時(shí)刻12以前均為無(wú)效信息）。

對(duì)比圖1與圖B-1，在初值或觀測(cè)值誤差較大的情況下，去掉負(fù)值糾錯(cuò)功能的算法提供的有效源反演信息很少。因此，考慮負(fù)值糾錯(cuò)功能使得算法健壯性大大增強(qiáng)。

圖B-1 試驗(yàn)1去掉負(fù)值糾錯(cuò)功能的結(jié)果：（a）濃度C1；（b）濃度C2；（c）Q1源排放速率Fig. B-1 The results of Experiment 1 without error correction function: (a) Concentration C1; (b) concentration C2; (c) emissions rate of Q1

資助項(xiàng)目 中國(guó)科學(xué)院戰(zhàn)略性先導(dǎo)灰霾專項(xiàng)XDB05030200

Funded by Strategic Priority Research Program (B) of the Chinese Academy of Sciences (Grant XDB05030200)

文章編號(hào)1006-9895(2016)04-0719-11 中圖分類號(hào) X11 P409

文獻(xiàn)標(biāo)識(shí)碼A

doi:10.3878/j.issn.1006-9895.1505.15113

收稿日期2015-01-19；網(wǎng)絡(luò)預(yù)出版日期 2015-09-21

作者簡(jiǎn)介馮帆，女，1981年出生，博士研究生，主要從事應(yīng)用數(shù)學(xué)建模及算法研究。E-mail: feng_fan_100@hotmail.com

Development of an Adaptive Algorithm Based on the Shooting Method and Its Application in the Problem of Estimating Air Pollutant Emissions

FENG Fan1, 2, 3, WANG Zifa2, and TANG Xiao2
1 Supercomputing Center, Computer Network Information Center, Chinese Academy of Sciences, Beijing 100190
2 State Key Laboratory of Atmospheric Boundary Layer Physics and Atmospheric Chemistry, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029
3 University of Chinese Academy of Sciences, Beijing 100049

AbstractThe inversion of pollutant emissions is important in air pollution prediction and control. The statistical methods generally adopted possess weaknesses, such as sensitivity for observation error and prior estimation of emissions. Given their simplicity, high precision, and practicality, algorithms based on the shooting method are widely used in the field of systems control. In this paper, we derive a shooting method–based adaptive algorithm to estimate air pollutant emissions. Besides its high precision and simple procedure, this adaptive algorithm compensates for the weaknesses of statistical methods. It is able to deal with the large level of error in the emissions inventory, initial conditions, and observations; a prior distribution assumption and error estimation are not required. Simulations of a simple system are presented to illustrate the effectiveness of this adaptive algorithm.

KeywordsShooting method, Pollutant emissions, Inversion, Adaptive algorithm, Systems control