利用GNSS PWV的AOD自适应预测方法

赵庆志，苏 静，杨鹏飞，姚宜斌

1. 西安科技大学测绘科学与技术学院，陕西 西安 710054； 2. 武汉大学测绘学院,湖北 武汉 430079

气溶胶是指悬浮在气体中的固体和液体微粒与气体载体共同组成的多相体系，其动力学直径大约是0.001～100 μm[1]。尽管气溶胶在地球大气中的组成非常小，但它对大气物理过程、大气环境和气候变化的影响较大[2]。气溶胶作为各类气象(雨、雪、冰雹等)形成的凝结核，在降水预报中发挥着重要的作用[3-4]。同时，气溶胶是大气污染的重要组成部分，在沙尘暴期间，污染物与烟雾、颗粒物密切相关[5-6]。此外，气溶胶会吸收太阳短波辐射和地球发出的长波辐射，影响地球大气系统的能量平衡和全球气候变化。在各种影响气候变化的因素中，气溶胶的气候效应目前最不清楚[7-8]。气溶胶光学厚度(aerosol optical depth，AOD)定义为介质的消光系数在垂直方向上的积分，用于描述气溶胶对光的削减作用。它是气溶胶最重要的参数之一，是表征大气浑浊程度的关键物理量，也是确定气溶胶气候效应的重要因素。因此，对AOD进行准确监测和估计对于研究全球辐射平衡和区域及全球气候变化具有重要意义[9]。

地面实测AOD是广泛应用于大气空气质量研究的重要数据源之一。如气溶胶机器人网络(aerosol robotic network,AERONET)和中国气溶胶遥感网络(China aerosol remote sensing network,CARSNET)，均可提供站点上高精度的AOD实测数据。但AERONET和CARSNET监测站点数量有限，空间分布过于稀疏，且主要集中于主城区，因此，无法在大尺度上独立描述AOD的空间变化[10]。与地面实测数据相比，卫星观测数据可有效弥补监测站点空间分辨率低的缺陷，且能够反映AOD的连续空间变化信息[11]。如搭载在Terra和Aqua卫星上的中尺度分辨率成像光谱仪(moderate-resolution imaging spectroadiometer,MODIS)，可以反演得到空间分辨率为3～10 km的全球气溶胶光学厚度数据。但由于极轨卫星系统自身的局限性，AOD产品的时间分辨率较低，加上云污染，限制了遥感卫星在高时间分辨率AOD变化研究中的应用[12]。此外，卫星遥感技术获取的AOD也可同化到数值再分析模型中，这大大提高了同化模型的精度[13-14]。美国国家航空航天局同化了1980年至今的现代回顾性分析(modern-era retrospective analysis for research and applications,version 2，MERRA-2)数据，对大气气溶胶进行了再分析研究和应用。MERRA-2比卫星观测技术提供了更好的AOD估计，覆盖范围广[15-16]。然而，由于同化系统的系统偏差，从MERRA-2获得的AOD在严重雾霾条件下误差较大[17]。除上述方法，还可以通过建立AOD与一些气象参数间的关系模型，反演得到高精度的AOD信息。文献[18]利用非线性主成分分析等方法得到影响AOD的几个综合指标，并联合气象参数建立AOD的预测模型。文献[19]建立了AOD与空气质量指数之间的线性回归模型，发现两者在春夏季存在良好的相关性。然而，上述模型所需建模参数较多，建模成本高且时间分辨率低，难以应用于实际。

大气水汽是对流层中一个重要的气象参数，其含量和变化可利用大气可降水量(precipitable water vapor,PWV)来反映。PWV是指地面到对流层顶部单位横截面积空气柱中所有水汽转化为液态水的含量[20]。自文献[21]提出全球导航卫星系统(Global Navigation Satellite System,GNSS)气象学以来，基于GNSS技术的PWV反演逐渐成为获取大气水汽的主要手段之一，已具有全天候、高精度、高时间分辨率、低成本等优点[22-23]。相关研究表明，大气水汽也会影响气溶胶的变化[24-25]。文献[26]发现GNSS PWV与MODIS AOD具有良好的相关性，其相关系数高达0.636。文献[27]提出了一种基于GNSS的PM2.5预测模型，该模型考虑了GNSS反演的天顶湿延迟(zenith wet delay,ZWD)。文献[28]进一步分析了GNSS反演的天顶干延迟(zenith total delay，ZTD)与PM2.5的关系，发现两者之间存在良好的相关性。上述研究中的ZWD或ZTD可以通过计算进一步得到PWV，这表明，GNSS反演的大气参数(ZTD、ZWD或PWV)与PM2.5具有很好的关系，可用于AOD监测的相关研究。然而，目前基于GNSS PWV的AOD建模研究较少。

因此，本文建立了两种基于GNSS PWV的AOD自适应预测模型。提出的模型顾及了相邻历元AOD间的时间自相关性，且模型系数能够自适应更新。第1种模型直接利用GNSS PWV对550 nm的AOD进行建模，简称TAF(total AOD forecast)模型。第2种模型在TAF模型的基础上，进一步考虑了不同类型AOD对PWV的敏感性，简称FTAF(five type-based AOD forecast)模型。与传统模型对比，本文提出的方法仅利用GNSS观测数据构建AOD模型，数据易于获取，节省成本，建模方法简单且精度较高，可为基于GNSS的AOD监测和预警研究提供理论支持。

1 数据与方法

1.1 数据描述

1.1.1 CMONOC

中国大陆构造环境监测网络(crustal movement observation network of China,CMONOC)是中国地壳运动观测网络的延续，始建于1997年，从2010年开始分批次安装调试，于2011年投入使用，共建有260个GNSS基准站，为GNSS气象学等研究领域提供了高精度的观测资料[29-30]。本文选取2012—2018年京津冀地区共16个GNSS测站进行研究，首先获取16个GNSS测站的PWV数据。该数据集是基于GNSS ZTD和气象参数计算得到，其中COMNOC的ZTD数据利用PANDA软件计算得到，温度和气压参数来自EAR5再分析产品。图1给出了选取GNSS站点和无线电探空站点(radiosonde,RS)在京津冀地区的地理分布。

1.1.2 ECMWF ERA5

ERA5是欧洲中期天气预报中心ECMWF(European Centre for Medium-Range Weather Forecasts)的第5代数据集，包含了1979年至今的数据集。与第4代产品ERA-Interim相比，ERA5增加了许多新变量，使用了更先进的同化技术，结合改进的湿度分析、卫星数据校正，使得ERA5可提供时间分辨率为1 h的大气、陆地和海洋气候变量的估计值[31]。本文选取ERA5数据集中的气温和气压数据，其时空分辨率分别为1 h和0.25°×0.25°，该数据集可通过https:∥www.ecmwf.int/下载。

图1 GNSS和RS站点在京津冀地区的地理分布Fig.1 Geographical distribution of GNSS and RS stations in the Beijing-Tianjin-Hebei Region

1.1.3 CAMS AOD

本文使用的AOD产品是由ECMWF最新的综合预报系统所衍生的哥白尼气溶胶监测服务(Copernicus aerosol monitoring service,CAMS)系统生成。该数据集包括2004年至今的数据，时间分辨率为3 h[32]。数据集主要包括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进行分析，相应的数据可以从www.soda-pro.com下载。

1.1.4 无线电探空

无线电探空技术是目前气象领域最常用的PWV探测方法之一。探空观测是通过释放搭载由天然乳胶制成的探空气球，借助安装在探空气球上的传感器来探测地面距高空30 km处不同高度上的气温、气压、湿度和风速等气象数据的观测技术，利用探空观测数据，通过积分的方法可得到站点上的PWV，相应数据可从ftp:∥ftp.ncdc.noaa.gov/pub/data/igra下载。

1.2 GNSS PWV获取

由于ERA5提供的气象参数和GNSS站点所在高度不同，首先，需在垂直方向上对ERA5的气压、气温数据进行高程改正，如式(1)和式(2)所示。然后，在水平方向上根据双线性插值得到GNSS站点上的温度和气压

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

(1)

dT/dh=0.006 5C/m

(2)

式中，P0和P分别代表格网点上的气压和站点上的气压数据；h0和h分别代表格网高和站点高；dT和dh分别为站点与格网间气温和高程的差值。

利用GNSS站点上的P计算天顶静力学延迟(zenith hydrostatic delay，ZHD)[33]

(3)

式中，P为地面气压，单位为hPa；φ为纬度，单位为rad；H为大地高，单位为km。然后，在ZTD中减去ZHD，即可得到ZWD。获得ZWD后，根据式(4)即可获得PWV

(4)

2 AOD自适应预测模型构建

2.1 PWV精度评估

为验证计算的GNSS PWV的精度，选取与GNSS测站并址的两个RS站进行验证。首先利用高程拟合经验公式将RS上的PWV统一到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的相关系数分别为0.82和0.87(Q<0.05)。一般认为，当两个变量之间的相关系数达到0.8以上，就可称这两组变量为强相关，相关性越强，表示验证数据的精度越高[35]，GNSS PWV与RS PWV相关性较强，而RS PWV属于实测数据，精度较高，因此，GNSS PWV精度较好。统计结果发现，两个站点上的平均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时间自相关性分析

AOD变化除了受PWV等气象因素的影响，还受相邻时段AOD的影响。为了便于描述，将相邻时段AOD间的相关性定义为AOD的时间自相关性。本节分析了2012—2018年4个测站天平均AOD的时间自相关性。图3给出了4个GNSS测站(HECD站、HELQ站、HEEY站和HETS站)相邻时段(天)AOD的时间自相关性热点图。由热力图可以看出两个属性之间的相关系数，颜色越深说明相关性越强，相关系数达到0.5以上，可认为两个属性之间显著相关[35]。由图3可知，4个GNSS站点上的AOD与上一历元的AOD值具有明显的相关性，4个GNSS站点相邻时段AOD的最大相关系数分别为0.75、0.82、0.84和0.74。此外，AOD的时间自相关性还表现出时间间隔越短，相关性越高的特点。因此，在建立AOD预测模型时应考虑上一历元AOD的变化情况。

图3 2012—2018年HECD站、HELQ站、HEEY站和HETS站相邻时段AOD的相关性热点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相关性分析

采用四分位间距法对16个GNSS站点上2012—2018年每3 h的AOD和PWV数据进行粗差剔除，并分析PWV与不同类型AOD间的相关性。图4给出了2012—2018年期间，16个GNSS站点上不同类型AOD和PWV数据的相关性统计图。由图4可知，PWV与AOD、BC AOD和SU AOD均具有较好的相关性(0.587、0.501和0.614)，16个站点的平均相关系数均在0.5以上(P<0.05)。而PWV与DU AOD、OM AOD、SS AOD之间的相关性较差，其相关系数分别为0.137、0.209和0.224。上述结果表明，不同类型AOD对PWV的敏感性不同，造成这一现象的原因是不同类型气溶胶的光学厚度占比不同。因此，在构建AOD 550 nm模型时，应考虑PWV对不同类型AOD数据的影响。

2.4 TAF和FTAF模型构建

本文首先建立直接利用GNSS PWV的AOD 550 nm预测模型，简称为TAF(total AOD forecast)模型，该模型考虑了AOD的时间自相关性，且能够自适应更新模型系数。通过上述分析发现，AOD 550 nm可通过5种类型的AOD求和得到，且不同类型的AOD对PWV的敏感性不同。因此，本文在TAF的基础上，进一步提出一种顾及5种类型AOD的AOD 550 nm预测(five-based total AOD forecast,FTAF)方法。该方法首先利用GNSS PWV分别针对5种类型的AOD建模；然后通过求和的方式得到最终的AOD 550 nm值。

图4 2012—2018年PWV与不同类型AOD的相关性Fig.4 Correlation between PWV and different types of AOD from 2012 to 2018

2.4.1 TAF模型构建

本文提出的TAF模型将PWV和上一历元的AOD作为输入参数，其表达式如下

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

a2·AODCAMS-AOD(i-1)

(6)

式中，a0、a1和a2为模型系数；PWV(i)表示第i时刻的PWV数据；AODCAMS-AOD(i-1)表示第i-1时刻的AOD数据。

通过式(7)得到初始AOD预测值与实际值之间的残差dAOD，利用LS方法分析发现dAOD同时具有年周期和半年周期。因此，进一步考虑AOD的季节变化因素，建立如下残差周期模型

(7)

式中,dAODmodel(i)为残差的周期模型值；A0为dAOD的平均值；A1、A2分别为dAOD的年周期和半年周期系数；φ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模型构建

基于TAF建模方法，分别对5种类型的AOD建模，然后通过加权求和的方法获取550 nm的AOD。其表达式如下

(9)

在AOD建模过程中，模型训练数据的长度会影响模型的质量。因此，在建模时，首先确定了模型训练数据的长度，分别以1 a和2 a的训练数据建立TAF模型和FTAF模型，发现训练数据的时间窗口为1 a时，建模得到的AOD的RMS、bias和MAE值最小。因此，在下面试验中，TAF模型和FTAF模型的训练数据的窗口值均为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数据与CAMS AOD数据均具有很好的符合性。相对于TAF模型，FTAF模型与CAMS AOD之间符合程度更好。图6给出了2014—2018年选取站点上CAMS AOD和TAF AOD、FTAF AOD的残差频率直方统计图。由图6可知，与TAF模型相比较，FTAF模型的AOD残差值较小、分布更集中在0附近。表1给出了2014—2018年期间，CAMS AOD和TAF AOD、FTAF AOD模型之间残差的统计结果。由表1可知，基于TAF模型预测的2014—2018年站的AOD的RMS、bias和MAE分别为0.24、0.07和0.19，FTAF预测的2014—2018年的AOD的RMS、bias和MAE分别为0.22、0.05和0.17。因此，本文提出的TAF模型和FTAF模型都具有很好的精度，且FTAF模型稍优于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 模型内外符合精度评定

为进一步验证本文模型的内外符合精度，图7和图8分别给出了2014—2018年期间16个GNSS站点TAF模型和FTAF模型的内外符合精度。由图7和图8可知，本文提出的TAF模型和FTAF模型的内外符合精度指标在16个测站上均较小，且FTAF模型的精度优于TAF模型。TAF模型和FTAF模型内符合精度的平均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。上述结果进一步验证了本文提出的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站点上的内符合精度的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模型对比

由上述分析可以看出，FTAF模型的精度优于TAF模型，这主要是由于FTAF模型考虑了不同类型AOD对PWV的敏感性。为了进一步分析FTAF模型相对于TAF模型的改善情况，图9给出了16个GNSS站点上FTAF模型相对于TAF模型的内外符合改善率分布。由图9可以看出，相对于TAF模型，FTAF模型能够有效改善GNSS测站上预测的AOD精度。统计结果显示，2014—2018年期间，16个GNSS站点上STD和RMS的改善率为6%～12%，平均改善率分别为9.7%和9.8%。上述结果进一步表明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 结 论

本文提出了两种基于GNSS PWV的AOD自适应预测模型。提出的模型以GNSS PWV作为唯一的外部输入参数，考虑了AOD与相邻时段AOD的时间自相关性且模型系数能够自适应更新。第1种模型直接利用GNSS PWV对550 nm AOD建模，第2种模型进一步考虑了不同类型AOD对PWV的敏感性。选取京津冀地区16个GNSS站点2014—2018年的数据进行试验，对提出两种模型的精度进行评估。试验结果表明，GNSS PWV与AOD具有相同的时序变化特征，二者间也展现出较好的相关性。两种模型预测的AOD与CAMS AOD长时序均具有很好的一致性，其中FTAF模型精度优于TAF模型。相较于TAF模型，FTAF模型的STD和RMS的改善率分别为9.7%和9.8%。本文提出的两种AOD预测模型仅利用GNSS PWV作为输入参数即可预测较高精度的AOD值，为大气环境质量监测和研究提供了一种新的手段。

致谢：感谢ECMEF提供的ERA5数据集和CAMS AOD数据，感谢CMONOC提供的GNSS ZTD数据集，感谢美国国家气候数据中心提供的IGRA2数据资料。