趙慶志，蘇 靜，楊鵬飛，姚宜斌

1. 西安科技大學測繪科學與技術(shù)學院，陜西 西安 710054； 2. 武漢大學測繪學院,湖北 武漢 430079

氣溶膠是指懸浮在氣體中的固體和液體微粒與氣體載體共同組成的多相體系，其動力學直徑大約是0.001～100 μm[1]。盡管氣溶膠在地球大氣中的組成非常小，但它對大氣物理過程、大氣環(huán)境和氣候變化的影響較大[2]。氣溶膠作為各類氣象(雨、雪、冰雹等)形成的凝結(jié)核，在降水預報中發(fā)揮著重要的作用[3-4]。同時，氣溶膠是大氣污染的重要組成部分，在沙塵暴期間，污染物與煙霧、顆粒物密切相關(guān)[5-6]。此外，氣溶膠會吸收太陽短波輻射和地球發(fā)出的長波輻射，影響地球大氣系統(tǒng)的能量平衡和全球氣候變化。在各種影響氣候變化的因素中，氣溶膠的氣候效應目前最不清楚[7-8]。氣溶膠光學厚度(aerosol optical depth，AOD)定義為介質(zhì)的消光系數(shù)在垂直方向上的積分，用于描述氣溶膠對光的削減作用。它是氣溶膠最重要的參數(shù)之一，是表征大氣渾濁程度的關(guān)鍵物理量，也是確定氣溶膠氣候效應的重要因素。因此，對AOD進行準確監(jiān)測和估計對于研究全球輻射平衡和區(qū)域及全球氣候變化具有重要意義[9]。

地面實測AOD是廣泛應用于大氣空氣質(zhì)量研究的重要數(shù)據(jù)源之一。如氣溶膠機器人網(wǎng)絡(aerosol robotic network,AERONET)和中國氣溶膠遙感網(wǎng)絡(China aerosol remote sensing network,CARSNET)，均可提供站點上高精度的AOD實測數(shù)據(jù)。但AERONET和CARSNET監(jiān)測站點數(shù)量有限，空間分布過于稀疏，且主要集中于主城區(qū)，因此，無法在大尺度上獨立描述AOD的空間變化[10]。與地面實測數(shù)據(jù)相比，衛(wèi)星觀測數(shù)據(jù)可有效彌補監(jiān)測站點空間分辨率低的缺陷，且能夠反映AOD的連續(xù)空間變化信息[11]。如搭載在Terra和Aqua衛(wèi)星上的中尺度分辨率成像光譜儀(moderate-resolution imaging spectroadiometer,MODIS)，可以反演得到空間分辨率為3～10 km的全球氣溶膠光學厚度數(shù)據(jù)。但由于極軌衛(wèi)星系統(tǒng)自身的局限性，AOD產(chǎn)品的時間分辨率較低，加上云污染，限制了遙感衛(wèi)星在高時間分辨率AOD變化研究中的應用[12]。此外，衛(wèi)星遙感技術(shù)獲取的AOD也可同化到數(shù)值再分析模型中，這大大提高了同化模型的精度[13-14]。美國國家航空航天局同化了1980年至今的現(xiàn)代回顧性分析(modern-era retrospective analysis for research and applications,version 2，MERRA-2)數(shù)據(jù)，對大氣氣溶膠進行了再分析研究和應用。MERRA-2比衛(wèi)星觀測技術(shù)提供了更好的AOD估計，覆蓋范圍廣[15-16]。然而，由于同化系統(tǒng)的系統(tǒng)偏差，從MERRA-2獲得的AOD在嚴重霧霾條件下誤差較大[17]。除上述方法，還可以通過建立AOD與一些氣象參數(shù)間的關(guān)系模型，反演得到高精度的AOD信息。文獻[18]利用非線性主成分分析等方法得到影響AOD的幾個綜合指標，并聯(lián)合氣象參數(shù)建立AOD的預測模型。文獻[19]建立了AOD與空氣質(zhì)量指數(shù)之間的線性回歸模型，發(fā)現(xiàn)兩者在春夏季存在良好的相關(guān)性。然而，上述模型所需建模參數(shù)較多，建模成本高且時間分辨率低，難以應用于實際。

大氣水汽是對流層中一個重要的氣象參數(shù)，其含量和變化可利用大氣可降水量(precipitable water vapor,PWV)來反映。PWV是指地面到對流層頂部單位橫截面積空氣柱中所有水汽轉(zhuǎn)化為液態(tài)水的含量[20]。自文獻[21]提出全球?qū)Ш叫l(wèi)星系統(tǒng)(Global Navigation Satellite System,GNSS)氣象學以來，基于GNSS技術(shù)的PWV反演逐漸成為獲取大氣水汽的主要手段之一，已具有全天候、高精度、高時間分辨率、低成本等優(yōu)點[22-23]。相關(guān)研究表明，大氣水汽也會影響氣溶膠的變化[24-25]。文獻[26]發(fā)現(xiàn)GNSS PWV與MODIS AOD具有良好的相關(guān)性，其相關(guān)系數(shù)高達0.636。文獻[27]提出了一種基于GNSS的PM2.5預測模型，該模型考慮了GNSS反演的天頂濕延遲(zenith wet delay,ZWD)。文獻[28]進一步分析了GNSS反演的天頂干延遲(zenith total delay，ZTD)與PM2.5的關(guān)系，發(fā)現(xiàn)兩者之間存在良好的相關(guān)性。上述研究中的ZWD或ZTD可以通過計算進一步得到PWV，這表明，GNSS反演的大氣參數(shù)(ZTD、ZWD或PWV)與PM2.5具有很好的關(guān)系，可用于AOD監(jiān)測的相關(guān)研究。然而，目前基于GNSS PWV的AOD建模研究較少。

因此，本文建立了兩種基于GNSS PWV的AOD自適應預測模型。提出的模型顧及了相鄰歷元AOD間的時間自相關(guān)性，且模型系數(shù)能夠自適應更新。第1種模型直接利用GNSS PWV對550 nm的AOD進行建模，簡稱TAF(total AOD forecast)模型。第2種模型在TAF模型的基礎(chǔ)上，進一步考慮了不同類型AOD對PWV的敏感性，簡稱FTAF(five type-based AOD forecast)模型。與傳統(tǒng)模型對比，本文提出的方法僅利用GNSS觀測數(shù)據(jù)構(gòu)建AOD模型，數(shù)據(jù)易于獲取，節(jié)省成本，建模方法簡單且精度較高，可為基于GNSS的AOD監(jiān)測和預警研究提供理論支持。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)描述

1.1.1 CMONOC

中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(crustal movement observation network of China,CMONOC)是中國地殼運動觀測網(wǎng)絡的延續(xù)，始建于1997年，從2010年開始分批次安裝調(diào)試，于2011年投入使用，共建有260個GNSS基準站，為GNSS氣象學等研究領(lǐng)域提供了高精度的觀測資料[29-30]。本文選取2012—2018年京津冀地區(qū)共16個GNSS測站進行研究，首先獲取16個GNSS測站的PWV數(shù)據(jù)。該數(shù)據(jù)集是基于GNSS ZTD和氣象參數(shù)計算得到，其中COMNOC的ZTD數(shù)據(jù)利用PANDA軟件計算得到，溫度和氣壓參數(shù)來自EAR5再分析產(chǎn)品。圖1給出了選取GNSS站點和無線電探空站點(radiosonde,RS)在京津冀地區(qū)的地理分布。

1.1.2 ECMWF ERA5

ERA5是歐洲中期天氣預報中心ECMWF(European Centre for Medium-Range Weather Forecasts)的第5代數(shù)據(jù)集，包含了1979年至今的數(shù)據(jù)集。與第4代產(chǎn)品ERA-Interim相比，ERA5增加了許多新變量，使用了更先進的同化技術(shù)，結(jié)合改進的濕度分析、衛(wèi)星數(shù)據(jù)校正，使得ERA5可提供時間分辨率為1 h的大氣、陸地和海洋氣候變量的估計值[31]。本文選取ERA5數(shù)據(jù)集中的氣溫和氣壓數(shù)據(jù)，其時空分辨率分別為1 h和0.25°×0.25°，該數(shù)據(jù)集可通過https:∥www.ecmwf.int/下載。

圖1 GNSS和RS站點在京津冀地區(qū)的地理分布Fig.1 Geographical distribution of GNSS and RS stations in the Beijing-Tianjin-Hebei Region

1.1.3 CAMS AOD

本文使用的AOD產(chǎn)品是由ECMWF最新的綜合預報系統(tǒng)所衍生的哥白尼氣溶膠監(jiān)測服務(Copernicus aerosol monitoring service,CAMS)系統(tǒng)生成。該數(shù)據(jù)集包括2004年至今的數(shù)據(jù)，時間分辨率為3 h[32]。數(shù)據(jù)集主要包括AOD 550 nm，BC AOD 550(AOD black carbon 550 nm)，DU AOD 550(AOD dust 550 nm)，OM AOD 550(AOD organic mater 550 nm)，SS AOD550(AOD sea salt 550 nm), SU AOD 550(AOD sulfate 550 nm)。其中，AOD 550 nm值等于其他5種AOD類型值的總和。本文選取上述6種類型的AOD進行分析，相應的數(shù)據(jù)可以從www.soda-pro.com下載。

1.1.4 無線電探空

無線電探空技術(shù)是目前氣象領(lǐng)域最常用的PWV探測方法之一。探空觀測是通過釋放搭載由天然乳膠制成的探空氣球，借助安裝在探空氣球上的傳感器來探測地面距高空30 km處不同高度上的氣溫、氣壓、濕度和風速等氣象數(shù)據(jù)的觀測技術(shù)，利用探空觀測數(shù)據(jù)，通過積分的方法可得到站點上的PWV，相應數(shù)據(jù)可從ftp:∥ftp.ncdc.noaa.gov/pub/data/igra下載。

1.2 GNSS PWV獲取

由于ERA5提供的氣象參數(shù)和GNSS站點所在高度不同，首先，需在垂直方向上對ERA5的氣壓、氣溫數(shù)據(jù)進行高程改正，如式(1)和式(2)所示。然后，在水平方向上根據(jù)雙線性插值得到GNSS站點上的溫度和氣壓

P=P0(1-0.000 022 6(h-h0))5.225

(1)

dT/dh=0.006 5C/m

(2)

式中，P0和P分別代表格網(wǎng)點上的氣壓和站點上的氣壓數(shù)據(jù)；h0和h分別代表格網(wǎng)高和站點高；dT和dh分別為站點與格網(wǎng)間氣溫和高程的差值。

利用GNSS站點上的P計算天頂靜力學延遲(zenith hydrostatic delay，ZHD)[33]

(3)

式中，P為地面氣壓，單位為hPa；φ為緯度，單位為rad；H為大地高，單位為km。然后，在ZTD中減去ZHD，即可得到ZWD。獲得ZWD后，根據(jù)式(4)即可獲得PWV

(4)

2 AOD自適應預測模型構(gòu)建

2.1 PWV精度評估

為驗證計算的GNSS PWV的精度，選取與GNSS測站并址的兩個RS站進行驗證。首先利用高程擬合經(jīng)驗公式將RS上的PWV統(tǒng)一到GNSS測站高度上

PWVh1=PWVh2·exp(-(h1-h2)/2000)

(5)

式中，PWVh1和PWVh2分別表示在高度h1和h2上對應的PWV。

圖2給出了2012—2018年期間HELY站(37.4°N，114.7°E)、HEEY站(40.1°N，114.2°E)GNSS與RS的PWV的散點密度圖。由圖2可以看出，GNSS和RS的PWV具有良好的一致性，兩個站點上PWV的相關(guān)系數(shù)分別為0.82和0.87(Q<0.05)。一般認為，當兩個變量之間的相關(guān)系數(shù)達到0.8以上，就可稱這兩組變量為強相關(guān)，相關(guān)性越強，表示驗證數(shù)據(jù)的精度越高[35]，GNSS PWV與RS PWV相關(guān)性較強，而RS PWV屬于實測數(shù)據(jù)，精度較高，因此，GNSS PWV精度較好。統(tǒng)計結(jié)果發(fā)現(xiàn)，兩個站點上的平均RMS、bias和MAE分別為0.89、0.63、0.63和1.07、0.75、-0.76，兩個站點處的RMS均小于3 mm，進一步說明本文利用的GNSS PWV具有較高的精度[36]。

圖2 2012—2018年HELY站、HEEY站GNSS和RS的PWV的散點密度Fig. 2 Scatter density of PWV of GNSS and RS at HELY and HEEY stations from 2012 to 2018

2.2 AOD時間自相關(guān)性分析

AOD變化除了受PWV等氣象因素的影響，還受相鄰時段AOD的影響。為了便于描述，將相鄰時段AOD間的相關(guān)性定義為AOD的時間自相關(guān)性。本節(jié)分析了2012—2018年4個測站天平均AOD的時間自相關(guān)性。圖3給出了4個GNSS測站(HECD站、HELQ站、HEEY站和HETS站)相鄰時段(天)AOD的時間自相關(guān)性熱點圖。由熱力圖可以看出兩個屬性之間的相關(guān)系數(shù)，顏色越深說明相關(guān)性越強，相關(guān)系數(shù)達到0.5以上，可認為兩個屬性之間顯著相關(guān)[35]。由圖3可知，4個GNSS站點上的AOD與上一歷元的AOD值具有明顯的相關(guān)性，4個GNSS站點相鄰時段AOD的最大相關(guān)系數(shù)分別為0.75、0.82、0.84和0.74。此外，AOD的時間自相關(guān)性還表現(xiàn)出時間間隔越短，相關(guān)性越高的特點。因此，在建立AOD預測模型時應考慮上一歷元AOD的變化情況。

圖3 2012—2018年HECD站、HELQ站、HEEY站和HETS站相鄰時段AOD的相關(guān)性熱點Fig.3 Correlation hot map of AOD in adjacent periods of HECD, HELQ, HEEY and HETS stations from 2012 to 2018

2.3 AOD與PWV相關(guān)性分析

采用四分位間距法對16個GNSS站點上2012—2018年每3 h的AOD和PWV數(shù)據(jù)進行粗差剔除，并分析PWV與不同類型AOD間的相關(guān)性。圖4給出了2012—2018年期間，16個GNSS站點上不同類型AOD和PWV數(shù)據(jù)的相關(guān)性統(tǒng)計圖。由圖4可知，PWV與AOD、BC AOD和SU AOD均具有較好的相關(guān)性(0.587、0.501和0.614)，16個站點的平均相關(guān)系數(shù)均在0.5以上(P<0.05)。而PWV與DU AOD、OM AOD、SS AOD之間的相關(guān)性較差，其相關(guān)系數(shù)分別為0.137、0.209和0.224。上述結(jié)果表明，不同類型AOD對PWV的敏感性不同，造成這一現(xiàn)象的原因是不同類型氣溶膠的光學厚度占比不同。因此，在構(gòu)建AOD 550 nm模型時，應考慮PWV對不同類型AOD數(shù)據(jù)的影響。

2.4 TAF和FTAF模型構(gòu)建

本文首先建立直接利用GNSS PWV的AOD 550 nm預測模型，簡稱為TAF(total AOD forecast)模型，該模型考慮了AOD的時間自相關(guān)性，且能夠自適應更新模型系數(shù)。通過上述分析發(fā)現(xiàn)，AOD 550 nm可通過5種類型的AOD求和得到，且不同類型的AOD對PWV的敏感性不同。因此，本文在TAF的基礎(chǔ)上，進一步提出一種顧及5種類型AOD的AOD 550 nm預測(five-based total AOD forecast,FTAF)方法。該方法首先利用GNSS PWV分別針對5種類型的AOD建模；然后通過求和的方式得到最終的AOD 550 nm值。

圖4 2012—2018年P(guān)WV與不同類型AOD的相關(guān)性Fig.4 Correlation between PWV and different types of AOD from 2012 to 2018

2.4.1 TAF模型構(gòu)建

本文提出的TAF模型將PWV和上一歷元的AOD作為輸入?yún)?shù)，其表達式如下

AOD550(i)=a0+a1·PWV(i)+

a2·AODCAMS-AOD(i-1)

(6)

式中，a0、a1和a2為模型系數(shù)；PWV(i)表示第i時刻的PWV數(shù)據(jù)；AODCAMS-AOD(i-1)表示第i-1時刻的AOD數(shù)據(jù)。

通過式(7)得到初始AOD預測值與實際值之間的殘差dAOD，利用LS方法分析發(fā)現(xiàn)dAOD同時具有年周期和半年周期。因此，進一步考慮AOD的季節(jié)變化因素，建立如下殘差周期模型

(7)

式中,dAODmodel(i)為殘差的周期模型值；A0為dAOD的平均值；A1、A2分別為dAOD的年周期和半年周期系數(shù)；φ1、φ2均為相位；DOY為一年中的年積日。

由式(6)和式(7)可得最終的TAF模型計算公式

AOD550(i)=a0+a1·PWV(i)+a2·

AODCAMS-AOD(i-1)+dAODmodel(i)

(8)

式中，AOD550(i)為第i時刻TAF模型的預測值。

2.4.2 FTAF模型構(gòu)建

基于TAF建模方法，分別對5種類型的AOD建模，然后通過加權(quán)求和的方法獲取550 nm的AOD。其表達式如下

(9)

在AOD建模過程中，模型訓練數(shù)據(jù)的長度會影響模型的質(zhì)量。因此，在建模時，首先確定了模型訓練數(shù)據(jù)的長度，分別以1 a和2 a的訓練數(shù)據(jù)建立TAF模型和FTAF模型，發(fā)現(xiàn)訓練數(shù)據(jù)的時間窗口為1 a時，建模得到的AOD的RMS、bias和MAE值最小。因此，在下面試驗中，TAF模型和FTAF模型的訓練數(shù)據(jù)的窗口值均為1 a，窗口值每個月更新一次。

3 精度評估

3.1 AOD時序分析

為了驗證本文提出的兩種AOD預測模型精度，對TAF模型和FTAF模型得到的長時序AOD進行分析。圖5給出了2014—2018年BJFS站、HECC站、TJBD站和JIXN站4個站點上CAMS AOD和TAF模型、FTAF模型預測的AOD時序圖。由圖5可知，兩種模型預測的AOD數(shù)據(jù)與CAMS AOD數(shù)據(jù)均具有很好的符合性。相對于TAF模型，F(xiàn)TAF模型與CAMS AOD之間符合程度更好。圖6給出了2014—2018年選取站點上CAMS AOD和TAF AOD、FTAF AOD的殘差頻率直方統(tǒng)計圖。由圖6可知，與TAF模型相比較，F(xiàn)TAF模型的AOD殘差值較小、分布更集中在0附近。表1給出了2014—2018年期間，CAMS AOD和TAF AOD、FTAF AOD模型之間殘差的統(tǒng)計結(jié)果。由表1可知，基于TAF模型預測的2014—2018年站的AOD的RMS、bias和MAE分別為0.24、0.07和0.19，F(xiàn)TAF預測的2014—2018年的AOD的RMS、bias和MAE分別為0.22、0.05和0.17。因此，本文提出的TAF模型和FTAF模型都具有很好的精度，且FTAF模型稍優(yōu)于TAF模型。

圖5 2014—2018年BJFS、HECC、TJBD和JIXN站點的CAMS AOD和TAF模型、FTAF模型預測的AOD的序列圖Fig.5 Sequence diagram of AOD predicted by CAMS AOD and TAF model、TAF model of BJFS, HECC, TJBD and JIXN stations from 2014 to 2018

表1 2014—2018年期間，CAMS AOD與TAF模型、FTAF模型之間AOD差異的RMS、bias和MAE

3.2 模型內(nèi)外符合精度評定

為進一步驗證本文模型的內(nèi)外符合精度，圖7和圖8分別給出了2014—2018年期間16個GNSS站點TAF模型和FTAF模型的內(nèi)外符合精度。由圖7和圖8可知，本文提出的TAF模型和FTAF模型的內(nèi)外符合精度指標在16個測站上均較小，且FTAF模型的精度優(yōu)于TAF模型。TAF模型和FTAF模型內(nèi)符合精度的平均STD、bias和MAE分別為0.23、0.08、0.19和0.21、0.05、0.17，外符合精度的平均RMS、bias和MAE分別為0.24、0.08、0.019和0.22、0.05、0.16。上述結(jié)果進一步驗證了本文提出的TAF模型和FTAF模型的可靠性。此外，隨著測站上AOD的降低(HEYY站至HECX站是按照AOD大小排序)，TAF模型和FTAF模型的精度越高，說明兩種模型均受AOD的影響，其值越低，預報精度越好。

圖6 2014—2018年期間BJFS站、HECD站、HETS站和TJBD站上的CAMS AOD與TAF模型、FTAF 模型的AOD殘差頻率Fig.6 Residual frequency of CAMS AOD and TAF AOD、FTAF AOD of BJFS, HECD, HETS and TJBD stations from 2014 to 2018

圖7 2014—2018年期間TAF模型和FTAF模型在16個GNSS站點上的內(nèi)符合精度的STD、MAE和biasFig.7 STD, MAE and bias of internal accuracy of TAF model and FTAF model at 16 GNSS stations during 2014—2018

圖8 2014—2018年期間TAF模型和FTAF模型在16個GNSS站點上的外符合精度的RMS、MAE和biasFig.8 STD, MAE and bias of external accuracy of TAF and FTAF model at 16 GNSS stations during 2014—2018

3.3 TAF模型和FTAF模型對比

由上述分析可以看出，F(xiàn)TAF模型的精度優(yōu)于TAF模型，這主要是由于FTAF模型考慮了不同類型AOD對PWV的敏感性。為了進一步分析FTAF模型相對于TAF模型的改善情況，圖9給出了16個GNSS站點上FTAF模型相對于TAF模型的內(nèi)外符合改善率分布。由圖9可以看出，相對于TAF模型，F(xiàn)TAF模型能夠有效改善GNSS測站上預測的AOD精度。統(tǒng)計結(jié)果顯示，2014—2018年期間，16個GNSS站點上STD和RMS的改善率為6%～12%，平均改善率分別為9.7%和9.8%。上述結(jié)果進一步表明FTAF模型具有更高的精度。

圖9 2014—2018年期間TAF模型和FTAF模型之間的RMS和STD改善率分布Fig.9 Distribution of RMS and STD improvement rate between TAF model and FTAF model from 2014 to 2018

4 結(jié) 論

本文提出了兩種基于GNSS PWV的AOD自適應預測模型。提出的模型以GNSS PWV作為唯一的外部輸入?yún)?shù)，考慮了AOD與相鄰時段AOD的時間自相關(guān)性且模型系數(shù)能夠自適應更新。第1種模型直接利用GNSS PWV對550 nm AOD建模，第2種模型進一步考慮了不同類型AOD對PWV的敏感性。選取京津冀地區(qū)16個GNSS站點2014—2018年的數(shù)據(jù)進行試驗，對提出兩種模型的精度進行評估。試驗結(jié)果表明，GNSS PWV與AOD具有相同的時序變化特征，二者間也展現(xiàn)出較好的相關(guān)性。兩種模型預測的AOD與CAMS AOD長時序均具有很好的一致性，其中FTAF模型精度優(yōu)于TAF模型。相較于TAF模型，F(xiàn)TAF模型的STD和RMS的改善率分別為9.7%和9.8%。本文提出的兩種AOD預測模型僅利用GNSS PWV作為輸入?yún)?shù)即可預測較高精度的AOD值，為大氣環(huán)境質(zhì)量監(jiān)測和研究提供了一種新的手段。

致謝：感謝ECMEF提供的ERA5數(shù)據(jù)集和CAMS AOD數(shù)據(jù)，感謝CMONOC提供的GNSS ZTD數(shù)據(jù)集，感謝美國國家氣候數(shù)據(jù)中心提供的IGRA2數(shù)據(jù)資料。