一種基于集合最優(yōu)插值的排放源快速反演方法

吳煌堅(jiān) 林偉 孔磊 唐曉 王威 王自發(fā) ,7 陳松蹊 ,

1 北京大學(xué)光華管理學(xué)院，北京 100871

2 北京大學(xué)數(shù)學(xué)科學(xué)學(xué)院，北京 100871

3 北京大學(xué)統(tǒng)計(jì)科學(xué)中心，北京 100871

4 中國科學(xué)院大氣物理研究所大氣邊界層物理和大氣化學(xué)國家重點(diǎn)實(shí)驗(yàn)室，北京 100029

5 中國科學(xué)院大學(xué)，北京 100049

6 中國環(huán)境監(jiān)測總站，北京 100012

7 中國科學(xué)院區(qū)域大氣環(huán)境研究卓越創(chuàng)新中心，福建廈門 361021

1 引言

大氣污染源排放清單是大氣污染預(yù)報(bào)最為重要的基礎(chǔ)數(shù)據(jù)之一，同時(shí)也是大氣化學(xué)傳輸模式的重要不確定性來源。Moore and Londergan（2001）定量研究了空氣質(zhì)量模擬的不確定性來源，結(jié)果表明排放源不確定性最大，占預(yù)報(bào)不確定性的25%～50%，高于氣象、邊界條件、初始濃度的影響。Beekmann and Derognat（ 2003） 和 Hanna et al.（1998）的研究同樣表明，排放源是模擬不確定性最重要的來源之一。傳統(tǒng)自下而上的源清單構(gòu)建方法通過排放因子和活動水平推算排放量（Hao et al.,2002）。該方法需要統(tǒng)計(jì)大量排放個(gè)體信息，包含每個(gè)點(diǎn)源的準(zhǔn)確位置。但受有限的資源和實(shí)地信息制約，排放清單更新緩慢且包含較大不確定性 （ Zhang et al.,2009; 曹 國 良 等 ,2011; 魏 巍 等 ,2011）。近年來發(fā)表的排放清單中（Zhang et al.,2009; 魏巍等,2011; Kurokawa et al.,2013; Li et al.,2017），細(xì)顆粒物（PM2.5）、可吸入顆粒物 （PM10）的排放不確定性在100%左右，一氧化碳（CO）和揮發(fā)性有機(jī)物（VOCs）的不確定性多在50%以上；部分清單的SO2不確定性雖可低至12%，但近年來脫硫設(shè)備和散煤管控等措施的實(shí)施使得SO2排放大幅下降，更新緩慢的排放清單難以滿足高精度預(yù)報(bào)預(yù)警和減排評估的需求。

隨著觀測技術(shù)的發(fā)展，越來越多研究者開始利用觀測數(shù)據(jù)結(jié)合源反演方法估算排放清單（Streets et al.,2013）。作為先進(jìn)的同化方法之一，集合卡爾曼濾波（EnKF）（Evensen,2003）已經(jīng)成為校驗(yàn)和優(yōu)化清單的重要方法。朱江和汪萍（2006）使用理想試驗(yàn)探討了利用EnKF開展排放源反演的可行性。Tang et al.（2013）使用該方法反演了北京及周邊的CO排放，并對比反演前后的模擬結(jié)果，發(fā)現(xiàn)源反演可將驗(yàn)證站點(diǎn)偏差降低48%。Miyazaki et al.（2017）通過同化多個(gè)衛(wèi)星資料反演多年的地面NOx排放，發(fā)現(xiàn)全球NOx排放總量在2005年[47.9 Tg(N) a?1]至 2014 年 [47.5 Tg(N) a?1]之間僅有少量變化。Peng et al.（2017）使用 EnKF 方法同時(shí)同化初始場和排放源，發(fā)現(xiàn)排放源反演可大幅提高長三角和珠三角的夜間PM2.5預(yù)報(bào)效果。

雖然基于EnKF的排放源反演方法可有效提高源清單精度，但該方法需多次運(yùn)行大氣化學(xué)傳輸模式，龐大的計(jì)算需求限制了其應(yīng)用范圍，使其無法為實(shí)時(shí)空氣質(zhì)量預(yù)報(bào)快速更新排放清單。為減少EnKF計(jì)算量，Evensen（2003）提出其次優(yōu)方案，即集合最優(yōu)插值（EnOI）。該方案利用不同時(shí)刻的模擬結(jié)果構(gòu)建集合狀態(tài)，避免多次運(yùn)行預(yù)報(bào)模式，從而大幅降低計(jì)算量。該方法在海洋資料同化中得到廣泛應(yīng)用（Counillon and Bertino,2009; Oke et al.,2010; Kaurkin et al.,2016）。張金譜等（2014）將EnOI方法應(yīng)用于珠三角的PM10、SO2和NO2濃度同化，使驗(yàn)證站點(diǎn)的均方根誤差降低32%～42%。然而常規(guī)的EnOI方法未擾動大氣化學(xué)傳輸模式，難以建立排放強(qiáng)度與污染物濃度之間的誤差協(xié)方差矩陣，故而無法直接應(yīng)用于排放源反演。為建立該聯(lián)系，Wang et al.（2016）對排放源進(jìn)行擾動，使用擾動后模擬的集合數(shù)據(jù)計(jì)算排放與濃度之間的誤差協(xié)方差矩陣，進(jìn)而反演黑碳排放。該方案每次反演均需人工合成反演時(shí)段的平均氣象場，并用其驅(qū)動24 h的集合模擬。對于長生命周期的污染物，該方法還需使用更多的計(jì)算資源延長集合模擬時(shí)段，從而使排放源擾動信息在模擬濃度中得到充分體現(xiàn)。

為降低源反演的計(jì)算量，本文發(fā)展了一種基于EnOI的排放源快速反演方法。該方法使用歷史集合數(shù)據(jù)構(gòu)建誤差協(xié)方差矩陣，再根據(jù)常規(guī)的空氣質(zhì)量觀測和模擬濃度反演排放源。相比基于EnKF的源反演方法和 Wang et al.（2016）的方法，本方法無需多次使用大氣化學(xué)傳輸模式模擬反演時(shí)段，取得可重復(fù)使用的歷史集合數(shù)據(jù)集后，僅需一次常規(guī)空氣質(zhì)量模擬即可快速更新排放清單。本文將該方法應(yīng)用于主要大氣污染物之一的CO排放源反演。雖然環(huán)境CO濃度很少超過世界衛(wèi)生組織的警戒值 （8 h平均濃度不超過10 mg m?3）（https://apps.who.int/iris/handle/10665/260127[2020-04-10]），但它是對流層臭氧的重要前體物，同時(shí)是大氣氧化性的主要控制因素（Dekker et al.,2017）。CO 的主要清除方式是被OH自由基氧化，其生命周期在幾周至幾個(gè)月之間（Holloway et al.,2000），長短適中的生命周期使其成為污染傳輸?shù)某Ｓ檬聚櫸铮∟aeher et al.,2001）。

2 方法和數(shù)據(jù)

2.1 基于集合最優(yōu)插值（EnOI）的源反演方法

在EnOI的多數(shù)應(yīng)用中，狀態(tài)變量與觀測變量相同，然而在源反演中，待更新的狀態(tài)變量為排放強(qiáng)度，與觀測的污染物濃度不同，故而無法將EnOI直接應(yīng)用于排放源反演。針對于此，本文借鑒EnKF的參數(shù)估計(jì)方法，將濃度向量與排放源強(qiáng)度的組合作為狀態(tài)變量：

其中，Xi為第i個(gè)狀態(tài)變量，Ci和Ei分別為第i個(gè)集合成員的濃度和排放強(qiáng)度，m為集合成員個(gè)數(shù)。

EnOI的狀態(tài)更新方法與EnKF基本一致，其核心區(qū)別在于EnKF多次使用大氣化學(xué)傳輸模式模擬同化時(shí)段，而EnOI使用歷史的模擬樣本構(gòu)建集合。在以往研究中，此集合常為不同時(shí)刻的模擬結(jié)果。然而此方案未擾動排放源，無法建立排放與濃度之間的誤差協(xié)方差矩陣。鑒于此，本文使用擾動排放源模擬得到的歷史集合數(shù)據(jù)構(gòu)建EnOI所需集合：

其中，為構(gòu)建的第i個(gè)背景狀態(tài)變量，Cb和Eb分別為反演時(shí)段的模擬污染物濃度和排放強(qiáng)度，和分別為歷史集合中第i個(gè)成員的模擬濃度和排放強(qiáng)度（歷史集合數(shù)據(jù)集在2.3節(jié)詳細(xì)介紹），和分別為歷史集合的濃度和排放的平均值。

構(gòu)建狀態(tài)集合后，EnOI可以使用與EnKF相同的狀態(tài)更新方法。本文采用 Sakov and Oke（ 2016）提出的確定性EnKF（DEnKF）方法進(jìn)行排放源反演。相比EnKF方法，它無需擾動觀測，避免在擾動過程中引入誤差。同時(shí)還支持基于舒爾（ Schur）積的局地化方法，且計(jì)算上更為高效。Sun et al.（2009）對比了幾種常見的集合濾波算法，結(jié)果表明DEnKF最穩(wěn)健，且在小集合的同化中精度最高。該方法的狀態(tài)更新公式如下：

其中，R為觀測誤差協(xié)方差矩陣，本文假設(shè)不同站點(diǎn)間的觀測誤差不相關(guān)，站點(diǎn)誤差估計(jì)采用李飛等 （2019）基于模擬分辨率的方案。P為背景場的誤差協(xié)方差矩陣，通過下式計(jì)算：

其中，? 代表舒爾積， ρ為局地化矩陣，用于抑制長距離虛假相關(guān)，其元素個(gè)數(shù)與AbAbT相同，各元素取值采用基于格點(diǎn)間距離的局地化函數(shù)（Gaspari and Cohn,1999）。Ab為背景狀態(tài)增量矩陣，表示集合成員與集合均值的差異：

其中矩陣1是單位矩陣I的同型矩陣，且各元素皆為1。

使用公式（4）可以更新集合狀態(tài)的均值，而狀態(tài)增量矩陣的更新公式如下：

其中Aa為更新后的狀態(tài)增量矩陣。結(jié)合公式（4）和公式（7）可得各集合成員的分析狀態(tài)：

其中為第i個(gè)集合成員的分析狀態(tài)，包含更新后的排放源。

該方法的核心優(yōu)勢在于可重復(fù)利用歷史集合數(shù)據(jù)，反演時(shí)段無需復(fù)雜的集合模擬和設(shè)置，僅需一次常規(guī)的空氣質(zhì)量模擬，即可反演模擬時(shí)段的排放源。故而該方法的計(jì)算量顯著低于基于EnKF的排放源反演方法，有助于拓展排放源反演的應(yīng)用范圍，提高排放源的更新速度。

2.2 數(shù)值模式介紹

本文的空氣質(zhì)量模擬采用王自發(fā)等（2006）發(fā)展的嵌套網(wǎng)格空氣質(zhì)量預(yù)報(bào)模擬系統(tǒng)（NAQPMS）。該系統(tǒng)涉及的物理化學(xué)過程主要包括平流、擴(kuò)散過程，干、濕沉降過程，氣相、液相和非均相化學(xué)過程。其空間結(jié)構(gòu)采用三維歐拉輸送模型，以地形追隨坐標(biāo)為垂直坐標(biāo)，可同時(shí)模擬區(qū)域和城市尺度沙塵、PM2.5、PM10、SO2、NOx、CO、O3、NH3等多種污染物（王自發(fā)等,2008）。其中氣相化學(xué)模塊 采 用 CBM-Z 機(jī) 制 （ Zaveri and Peters,1999），考慮了71種化學(xué)反應(yīng)物種，133個(gè)化學(xué)反應(yīng)；液相化學(xué)模塊基于RADM的反應(yīng)機(jī)制，包含22種氣體和氣溶膠（Chang et al.,1987）。NAQPMS 的氣溶膠包括硫酸鹽、硝酸鹽、銨鹽和二次有機(jī)氣溶膠。該系統(tǒng)使用 ISORROPIA（Nenes et al.,1998）計(jì)算氣溶膠的無機(jī)化學(xué)組分在顆粒態(tài)與氣態(tài)之間的分配及氣溶膠的含水量。

本文中，NAQPMS模擬所需氣象場由中尺度天氣預(yù)報(bào)模式（WRF）提供。氣象模擬采用圖1所示的兩層嵌套設(shè)置。其中第一層嵌套（D1）的水平分辨率為45 km×45 km，覆蓋了亞洲大部分區(qū)域；第二層嵌套（D2）的分辨率為15 km×15 km，覆蓋了除南海南部以外的中國區(qū)域。本研究以WRF連續(xù)模擬36 h、取后24小時(shí)的方式為NAQPMS模擬提供每日的氣象場驅(qū)動，再將每日氣象場拼接得到模擬時(shí)段的完整氣象場。NAQPMS模擬采用單層嵌套設(shè)置，模擬D2層15 km分辨率的各污染物濃度。

圖1 區(qū)域設(shè)置Fig.1 Domain configuration

本文使用清華大學(xué)提供的2010年MIX清單作為先驗(yàn)排放源（Li et al.,2017）。該清單集成了MEIC（中國大陸）、REAS2（日本、臺灣）、PKU-NH3（中國氨排放清單）、CAPSS（韓國）等亞洲各地區(qū)的排放清單，涵蓋10種大氣污染物和溫室氣體，網(wǎng)格分辨率為 0.25°（緯度）×0.25°（經(jīng)度）。該清單被用于東亞模式比較計(jì)劃第三期（MICSAsia III）的模擬研究，并且被聯(lián)合國半球大氣污染傳輸計(jì)劃（HTAP）采用（http://www.meicmodel.org[2020-04-10]）。

2.3 歷史集合數(shù)據(jù)集

本文通過NAQPMS模式模擬的歷史集合數(shù)據(jù)集建立濃度與排放源之間的聯(lián)系，進(jìn)而使用觀測濃度反演排放源。為使該數(shù)據(jù)集反映排放誤差對濃度的影響，需對排放源進(jìn)行擾動。若直接對每個(gè)網(wǎng)格的排放施加擾動，遠(yuǎn)小于格點(diǎn)數(shù)的擾動樣本將包含過大的采樣誤差。鑒于此，本文采用Evensen （2003）的偽隨機(jī)擾動方法生成500組均值為0、標(biāo)準(zhǔn)差為1的偽隨機(jī)擾動樣本。該樣本在每個(gè)空間格點(diǎn)上服從標(biāo)準(zhǔn)正態(tài)分布，每個(gè)樣本在空間上均較為平滑，且格點(diǎn)間樣本的相關(guān)性隨距離遞減。當(dāng)格點(diǎn)間距離達(dá)到解相關(guān)距離時(shí)，樣本的理想相關(guān)系數(shù)為e?1。若解相關(guān)距離過大，反演算法難以差異化調(diào)整臨近排放源；若解相關(guān)系數(shù)過小，有限的樣本將產(chǎn)生顯著的遠(yuǎn)距離虛假相關(guān)。本文使用試錯(cuò)法在不產(chǎn)生顯著虛假相關(guān)的前提下減小解相關(guān)距離，最終將解相關(guān)距離設(shè)為150 km。為減少NAQPMS模擬次數(shù)，本文采用Evensen（2004）提出的算法將500組擾動樣本壓縮為50組。以上先過量采樣再壓縮的方法，可以在不增加樣本量的前提下減小抽樣誤差，規(guī)避小概率樣本，使最終生成的樣本在每個(gè)格點(diǎn)上更接近標(biāo)準(zhǔn)正態(tài)分布。根據(jù)先驗(yàn)排放源的不確定性（Li et al.,2017），本文將壓縮后的擾動樣本均值調(diào)整為1，將標(biāo)準(zhǔn)差按表1調(diào)整。使用調(diào)整后的樣本擾動先驗(yàn)排放源，得到集合模擬所需的擾動排放清單，即：

其中，Eo為先驗(yàn)排放源，和分別為第i個(gè)集合成員的擾動樣本和擾動排放源。將擾動后的50組排放源分別使用NAQPMS模擬，得到對應(yīng)的50 組模擬濃度，取其時(shí)間平均，即得（i=1,…,50）。

表1 排放源擾動系數(shù)標(biāo)準(zhǔn)差Table 1 Standard deviations of the emission perturbation coefficients

在歷史集合數(shù)據(jù)集的模擬中，本研究采用與2.2節(jié)相同的模式和區(qū)域設(shè)置。模擬時(shí)段為2013年12月17日至2014年1月31日和2014年12月17日至2015年1月31日。從兩時(shí)段中分別取2014年1月和2015年1月的模擬月均濃度和擾動排放源構(gòu)建EnOI所需的背景狀態(tài)集合。值得指出的是，該數(shù)據(jù)可重復(fù)應(yīng)用于EnOI源反演，反演時(shí)段所采用的模式和區(qū)域設(shè)置均無需與歷史數(shù)據(jù)集相同。

2.4 試驗(yàn)設(shè)置

基于EnOI的源反演方法根據(jù)模擬濃度誤差調(diào)整排放源。為反演2015年1月全國的CO排放，本文首先使用先驗(yàn)排放源模擬2014年12月25日至2015年1月31日的全國空氣質(zhì)量，再使用兩組集合數(shù)據(jù)集分別進(jìn)行兩組反演試驗(yàn)。反演試驗(yàn)1使用2014年1月歷史集合數(shù)據(jù)集構(gòu)建濃度與排放之間的誤差協(xié)方差矩陣，而反演試驗(yàn)2使用2015年1月的集合數(shù)據(jù)集。兩組試驗(yàn)均根據(jù)2015年1月同化站點(diǎn)的CO觀測模擬濃度差異，反演2015年1月全國15 km分辨率的CO排放。反演中，觀測數(shù)據(jù)來自中國環(huán)境監(jiān)測總站，站點(diǎn)分布與人口分布基本一致，主要集中在胡煥庸線以東地區(qū)。為保障數(shù)據(jù)質(zhì)量，本研究先采用一種基于概率分布的自適應(yīng)離群值檢測方法剔除異常數(shù)據(jù)（Wu et al.,2018），再排除有效數(shù)據(jù)較少的站點(diǎn)，得到全國共1456個(gè)站點(diǎn)的觀測數(shù)據(jù)集。為客觀評估反演效果，本研究在包含兩個(gè)以上觀測站點(diǎn)的城市中隨機(jī)選取一個(gè)站點(diǎn)作為獨(dú)立驗(yàn)證站點(diǎn)，最終從全國1456個(gè)站點(diǎn)中選取349個(gè)站點(diǎn)作為驗(yàn)證站點(diǎn)，其余1107個(gè)站點(diǎn)作為同化站點(diǎn)。兩類站點(diǎn)的分布如圖2所示。

為驗(yàn)證源反演效果，本文使用反演的排放源再次模擬2014年12月25日至2015年1月31日的全國空氣質(zhì)量。通過對比兩套排放源在獨(dú)立驗(yàn)證站點(diǎn)的CO模擬結(jié)果，定量評估基于EnOI的源反演對CO模擬的改進(jìn)效果。

圖2 空氣質(zhì)量監(jiān)測站點(diǎn)分布Fig.2 Distribution of air quality monitoring stations

3 結(jié)果

本文首先使用先驗(yàn)排放源模擬2015年1月全國空氣質(zhì)量。圖3a對比了地面月均CO模擬和觀測濃度的空間分布。由圖可見，采用先驗(yàn)排放源模擬的CO濃度在全國大部分地區(qū)顯著低于觀測值，特別是在山西、吉林、遼寧等地和珠三角周邊的中國南部地區(qū)。全國349個(gè)驗(yàn)證站點(diǎn)的月均濃度偏差達(dá)?0.74 mg m?3。圖4 使用箱型圖對比了全國驗(yàn)證站點(diǎn)每日的觀測和模擬濃度。由圖可見，模擬的CO濃度變化趨勢與觀測較為一致，峰值均出現(xiàn)在4日、10日和25日，低值均出現(xiàn)在6日和27日。但模擬濃度的中位數(shù)和變化范圍普遍低于觀測濃度。圖5給出了北京、太原、武漢、廣州4個(gè)城市驗(yàn)證站點(diǎn)的進(jìn)一步分析（圖6a展示了2010年1月的先驗(yàn)CO排放強(qiáng)度）。4個(gè)城市中，北京市的CO模擬相關(guān)系數(shù)最高，達(dá)0.81，日均濃度偏差為?0.35 mg m?3；太原市的模擬濃度存在顯著低估，偏差達(dá)?2.16 mg m?3；廣州市也有明顯低估，偏差為?0.69 mg m?3；武漢市的模擬結(jié)果則有明顯高估，偏差為 1.20 mg m?3。

模擬濃度低于觀測反映出先驗(yàn)排放強(qiáng)度可能低于真實(shí)排放，反之亦然。反演試驗(yàn)1根據(jù)2015年1月全國1115個(gè)站點(diǎn)的模擬誤差，結(jié)合2014年1月的集合數(shù)據(jù)集，使用EnOI方法反演得到2015年1月的排放源調(diào)整系數(shù)（如圖6b）。由于觀測站點(diǎn)未覆蓋所有區(qū)域，故而臺灣、新疆、西藏、青海和內(nèi)蒙古的部分區(qū)域未進(jìn)行反演調(diào)整。但除臺灣外，大部分高排放地區(qū)均得到有效反演，未反演區(qū)域的先驗(yàn)排放強(qiáng)度普遍較低，對反演誤差影響較小。反演后大部分地區(qū)的排放強(qiáng)度均有所增加，全國2015年1月的CO排放量為46.6 Tg，比先驗(yàn)排放量（21.2 Tg）高120%。在CO濃度被顯著低估的山西和新疆等地，反演的排放強(qiáng)度可比先驗(yàn)排放源高3倍以上。

為評估源反演效果，本文使用反演的排放源再次模擬全國CO濃度。由圖3b可見，源反演可顯著提高CO模擬濃度，從而將偏差幅度由0.74 mg m?3降至 0.01 mg m?3。其中，山西、甘肅、兩廣地區(qū)和東北三省的CO模擬偏差均有顯著降低。與此同時(shí)，源反演可將全國驗(yàn)證站點(diǎn)的平均均方根誤差降低18%，相關(guān)系數(shù)提升0.03。在圖4的箱形圖對比中，反演后模擬濃度的中位數(shù)和變化范圍均顯著增大，與觀測高度一致。在4城市驗(yàn)證站點(diǎn)的濃度對比（圖5）中，源反演將CO模擬的相關(guān)系數(shù)提高0.02～0.20。在太原市驗(yàn)證站點(diǎn)中，源反演將CO 濃度提高 1.97 mg m?3，偏差幅度降低 91%。與此同時(shí)，廣州市驗(yàn)證站點(diǎn)的偏差幅度也大幅降低，降幅達(dá)88%。在CO濃度被高估的武漢，源反演可降低其排放強(qiáng)度和模擬濃度，同時(shí)大幅抑制9～11日的虛假峰值。

圖3 （a）反演前和（b）反演試驗(yàn)1得到的2015年1月CO模擬濃度分布。反演試驗(yàn)1使用2014年1月的集合數(shù)據(jù)集，結(jié)合反演前2015年1月的模擬濃度誤差，估計(jì)2015年1月的CO排放。圓點(diǎn)顏色表示驗(yàn)證站點(diǎn)的觀測月均濃度。驗(yàn)證站點(diǎn)的平均偏差（BIAS）、均方根誤差（RMSE）和標(biāo)準(zhǔn)差（r）列于圖中Fig.3 Simulated CO concentrations (a) before and (b) after inversion test 1 in Jan 2015.Inversion test 1 estimates the CO emission in Jan 2015 with the ensemble dataset of Jan 2014 and the simulated error of Jan 2015 before inversion.Colors of circles indicate the observed monthly mean concentrations at validation sites.Biases (BIAS),root-mean-square errors (RMSE),and correlations (r) over all validation sites are also shown

為減小計(jì)算量，本文的源反演方法忽略了歷史時(shí)段與反演時(shí)段的誤差協(xié)方差差異，使用歷史集合數(shù)據(jù)集構(gòu)建表征排放與濃度切線性關(guān)系的誤差協(xié)方差矩陣。歷史時(shí)段與反演時(shí)段的氣象條件差異可影響誤差協(xié)方差矩陣，進(jìn)而影響反演結(jié)果。然而，反演中使用的污染物平均濃度仍取自反演時(shí)段的觀測和模擬結(jié)果，因此反演結(jié)果可以體現(xiàn)反演時(shí)段的氣象和污染特征，從而限制兩個(gè)時(shí)段氣象條件差異對結(jié)果的影響。為定量評估該影響，本文在反演試驗(yàn)2中使用2015年1月的集合數(shù)據(jù)集再次構(gòu)建誤差協(xié)方差矩陣并反演排放源。結(jié)果發(fā)現(xiàn)，使用2015年1月集合數(shù)據(jù)集反演的全國CO排放量為46.0 Tg，僅比使用2014年1月集合數(shù)據(jù)集反演的結(jié)果小1%。在太原市和廣州市的驗(yàn)證站點(diǎn)中（圖5），兩套反演源的模擬均方根誤差差異均小于0.01 mg m?3。在全國349個(gè)驗(yàn)證站點(diǎn)中，使用反演時(shí)段集合數(shù)據(jù)集的平均均方根誤差為 0.78 mg m?3，僅比使用歷史集合數(shù)據(jù)集低 0.02 mg m?3，占反演前均方根誤差的2%。以上結(jié)果表明，2014年1月與2015年1月的氣象差異對反演結(jié)果影響有限，使用歷史集合數(shù)據(jù)集代替反演時(shí)段集合數(shù)據(jù)集可有效反演排放源。

圖4 反演前后2015年1月全國驗(yàn)證站點(diǎn)的每日CO濃度箱型圖（圓點(diǎn)表示中位數(shù)，箱形的頂部和底部分別表示75和25百分位數(shù)）Fig.4 Boxplot of daily CO concentrations before and after emission inversion over all validation sites in Jan 2015 (central mark indicates the median,and the top and bottom edges of the box indicate the 25th and 75th percentiles,respectively)

圖5 反演前后站點(diǎn)的觀測和模擬時(shí)間序列。反演試驗(yàn)1和反演試驗(yàn)2分別表示使用2014年1月和2015年1月集合數(shù)據(jù)集進(jìn)行反演，得到的2015年1月模擬結(jié)果。各站點(diǎn)的偏差（BIAS）、均方根誤差（RMSE）和標(biāo)準(zhǔn)差（r）列于圖中Fig.5 Time series of simulated and observed CO concentrations before and after emission inversion with corresponding biases (BIAS),root-meansquare errors (RMSE),and correlations (r).Inversion test 1 and inversion test 2 represent the concentrations in Jan 2015 by the inversions with ensemble dataset of Jan 2014 and Jan 2015,respectively

綜上所述，使用先驗(yàn)排放源模擬的2015年1月CO濃度存在顯著偏差。本文使用基于EnOI的源反演方法同化站點(diǎn)觀測數(shù)據(jù)，結(jié)合2014年1月的歷史集合數(shù)據(jù)集，快速反演全國15 km分辨率的CO排放。該方案反演的全國CO排放僅比使用2015年1月集合數(shù)據(jù)集反演的結(jié)果高1%，表明歷史時(shí)段與反演時(shí)段的氣象條件差異只能有限地影響月均CO排放的反演結(jié)果。使用反演的排放源再次模擬，可大幅降低CO模擬偏差，同時(shí)顯著降低均方根誤差，提高相關(guān)系數(shù)，表明基于EnOI的快速源反演方法可有效降低排放源的不確定性，提高空氣質(zhì)量的模擬和預(yù)報(bào)精度。

圖6 （a）CO 先驗(yàn)排放強(qiáng)度與（b）反演調(diào)整系數(shù)Fig.6 (a) A priori CO emissions and (b) adjustment factors from inversion

4 結(jié)論

本文發(fā)展了一種基于EnOI的大氣排放源反演方法。該方法根據(jù)模擬的濃度誤差調(diào)整排放強(qiáng)度，可降低排放清單的不確定性，為空氣質(zhì)量預(yù)報(bào)快速更新排放源。相比基于EnKF的源反演方法，本文方法無需多次使用大氣化學(xué)傳輸模式模擬反演時(shí)段，僅需一次未擾動排放源的常規(guī)模擬即可快速更新排放清單，將排放源反演計(jì)算量降低至少一個(gè)數(shù)量級。

本文將該方法應(yīng)用于2015年1月全國15 km分辨率的CO排放源反演。相比基于2010年1月的先驗(yàn)CO排放，反演的全國2015年1月CO排放比先驗(yàn)排放源高120%，且排放強(qiáng)度在大部分地區(qū)均有所增強(qiáng)。在CO模擬濃度被顯著低估的山西、新疆等地，反演的排放強(qiáng)度可達(dá)先驗(yàn)源的3倍以上。使用2014年1月集合數(shù)據(jù)集反演的全國CO排放僅比使用2015年1月集合數(shù)據(jù)集的反演結(jié)果高1%，表明歷史時(shí)段與反演時(shí)段的氣象條件差異只能有限地影響月均CO排放的反演結(jié)果。為評估排放源反演效果，本文對比了反演前后全國349個(gè)獨(dú)立驗(yàn)證站點(diǎn)的CO濃度，發(fā)現(xiàn)使用歷史數(shù)據(jù)集反演的排放源可顯著降低模擬偏差，將全國349個(gè)獨(dú)立驗(yàn)證站點(diǎn)的平均低估從 0.74 mg m?3降至0.01 mg m?3，同時(shí)將均方根誤差降低18%。

本文將基于EnOI的污染源反演方法成功應(yīng)用于CO排放源反演，后續(xù)研究可嘗試將其應(yīng)用于黑碳、NO2和一次顆粒物等已使用EnKF成功反演的物 種 （ Peng et al.,2017; Miyazaki et al.,2017;Houweling et al.,2017）。然而作為 EnKF 的次優(yōu)方案，EnOI未使用觀測數(shù)據(jù)在線更新排放源。排放源的誤差可影響線性化的誤差協(xié)方差估計(jì)，故而使用本文方法反演強(qiáng)非線性的污染源可能存在較大的不確定性（如使用臭氧濃度反演揮發(fā)性有機(jī)物排放）。此外，本文方法使用歷史集合數(shù)據(jù)集估計(jì)反演時(shí)段的誤差協(xié)方差矩陣，兩時(shí)段的氣象條件差異可影響反演結(jié)果，使其難以適用于極端天氣或高時(shí)間分辨率（如一天）的污染源反演。盡管難以適用于上述極端情形，本文方法可將污染源反演的計(jì)算量降低一個(gè)數(shù)量級并且適用于多種污染物的月均排放源反演。此外，為進(jìn)一步提高反演精度，可使用更新的排放源再次模擬，并根據(jù)新的模擬濃度誤差再次使用基于EnOI的源反演方法，通過反復(fù)迭代以應(yīng)對源反演中普遍存在的非線性和偏差問題。后續(xù)研究可進(jìn)一步對比相同計(jì)算量下，基于EnKF的源反演方法與基于EnOI的迭代源反演的反演精度，以做出更優(yōu)選擇。