一种基于OMI数据的中国区域对流层O3的反演方法

徐军 卢毅敏

0 引言

臭氧(O3)是大气中的一种痕量气体,约占大气总量的十万分之一[1]．平流层臭氧是地球的“保护伞”,能够使人类免受紫外线的干扰,在大气环境变化中发挥着重要的作用[2-3]．与平流层臭氧相比,高浓度的对流层臭氧严重影响着动植物的生长和人们的身体健康,会引发人们的呼吸性系统疾病,尤其对婴幼儿健康的损害更为严重[4-6]．对流层臭氧还是一种重要的温室气体,其增温效果仅次于CO2和SO2[7]．近年来,随着国家在PM2.5治理上的投入,PM2.5质量浓度逐渐降低,空气质量逐渐改善,但臭氧已逐渐代替PM2.5成为我国大部分地区的首要污染物[8-9]．因此,为了能够精确地评估对流层臭氧变化对人类身体健康、动植物生存的危害，以及对气候变化等问题产生的影响,需要一套精度较高、时间序列较长的臭氧数据来建模与分析,而现有的对流层臭氧产品无法满足较高时空分辨率的监测要求．

相对于地面站点数据,卫星能够获取到高空间分辨率、高精度的全球的臭氧数据．从20世纪70年代开始,美国就通过臭氧总量绘图光谱仪(Total Ozone Mapping Spectrometer,TOMS)获取全球尺度上的臭氧时空变化信息[10]．随着空间遥感技术的不断进步,一系列新型环境遥感卫星的应用也越来越广泛．由于臭氧监测仪(Ozone Monitoring Instrument,OMI)具有更高的空间分辨率,1 d能够获取覆盖全球的数据,且具有较长的时间序列,因此得到了广泛的使用．但OMI数据自2006年末开始就出现了数据缺失的现象,且随着时间的推移,缺失内容越来越多,严重影响了数据质量,因此需要对数据进行修复[11]．目前,针对臭氧数据修复的研究方法较少,常见的数据修复方法就是空间插补法,如:样条函数法[12]、反距离加权法[13]、克里金插值法[14]、贝叶斯插值法[15-16]等．此外也有学者利用多元线性回归[17]、卷积神经网络[18]、生成对抗网络[19-20]等进行多元数据的融合,从而实现对遥感数据的修复．基于OMI数据的特性,在短时间、小范围内的臭氧变化波动较小,因此本文根据时空拟合的方法实现对缺失数据的修复[21-22]．

对流层臭氧的反演方法有对流层云切片法[23]、卫星观测角差异法[24]、直接反演法(最优估计算法)[25]．卫星观测角差异法不能到达中国区域的对流层顶部,无法准确地获得平流层部分的臭氧,直接反演法结果会偏高,而云反演法一般适用于热带,由于云的物理和微物理特性随季节变化,不能将卫星云产品用于中国区域对流层臭氧浓度的反演．为此,本文选择对流层臭氧残差法实现对流层臭氧总量的提取,从而获取具有更高精度的对流层臭氧数据,为对流层臭氧污染的评估与防治提供重要的数据．

1 研究数据与方法

1.1 研究数据

1.1.1 OMI臭氧数据

臭氧监测仪(OMI)传感器是NASA于2004年发射的Aura卫星上搭载的传感器,主要用于监测痕量气体的变化,其空间分辨率为13 km×24 km,时间分辨率为1 d,可以提供每日的全球范围的臭氧观测结果．OMI可以观测臭氧柱总量和廓线数据、气溶胶、云量等,以及NO2、SO2等痕量气体．OMI产品共有Level1B、Level2、Level2G、Level3四种等级的产品,本文选择基本OMI-TOMS V8算法的Level3产品OMTO3e,它是空间分辨率为0.25°×0.25°,时间分辨率为1 d的臭氧总量产品．为了避免云量数据的干扰,本文剔除了云量大于0.3的数据格点．

图1 MLS臭氧廓线数据与OMI臭氧廓线数据的比较(2019-01-01)Fig.1 MLS ozone profile (a) and OMI ozone profile (b) on Jan.1st,2019

1.1.2 MLS臭氧廓线数据

微波临边探测仪(Microwave Limb Sounder,MLS)是搭载在Aura卫星上的临边观测仪,用于提供温度和O3、NO2、CO等痕量气体的廓线数据．该数据的垂直分辨率约为3 km,水平分辨率约为165 km×3 km,从地面(1 000 hPa)开始将数据分为55层,每层代表不同的气压值．虽然OMI也提供了臭氧廓线数据,其从地面1 013～0.3 hPa将大气层分为18层,因此MLS数据相较于OMI臭氧廓线数据具有更高的垂直分辨率(图1),其中MLS臭氧廓线(图1a)为体积混合比的形式,它定义为某一气体组分占整体空气体积的分数;OMI臭氧廓线(图1b)用多布森(Dobson)单位表示,它表示在标准大气状态下10-3cm臭氧层的厚度,单位为DU．

Froidevaux等[26]最早对MLS臭氧廓线数据的精度进行验证发现,MLS数据在平流层各个压强范围内的误差在1%左右,而对流层的误差则大于20%,因此本文所使用的数据主要是平流层范围内的数据．

1.1.3 对流层顶高度数据

对流层顶的高度数据主要选用美国国家环境预报中心和美国国家大气研究中心的NCEP/NCAR对流层顶高度再分析资料的每日数据,该数据的空间分辨率为2.5°×2.5°,时间分辨率为1 d,下载自http:∥www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.tropopause.html．

1.1.4 OMI/MLS对流层臭氧数据

OMI/MLS对流层臭氧数据是根据TOR(Tropospheric Ozone Residue)方法反演得到的对流层臭氧月平均数据集(https:∥acd-ext.gsfc.nasa.gov/Data_services/cloud_slice/new_data.html),该数据的空间分辨率为1.25°×1°,覆盖了60°S～60°N之间的全球对流层臭氧数据,且该数据与地基数据相比具有较好的一致性．

1.2 研究方法

1.2.1 时空加权法

由于OMI臭氧数据可以获取每日的臭氧柱总量数据,且在短时间、小范围内,臭氧总量的变化只有轻微的波动,因此可用时空拟合的方法来对臭氧总量进行修复．时空拟合的核心思想就是利用该时刻相邻两个时刻的数据作为参考数据(图2),图2中t±1时刻表示与t时刻相邻两个时相的数据,白底表示数据缺失,灰底表示数据有效,t时刻的红色三角形表示待修复像元,蓝色的三角形为参考像元,且t中的缺失区域能完全被t-1和t+1两个时刻的有效值覆盖．通过与目标数据建立相应的回归关系对缺失值进行预测,然后依据两个参考像元对目标数据的贡献程度进行加权求和,最终得到目标数据的拟合值．将目标数据的位置作为中心点,通过多组实验对比之后,确定一个以目标数据为中心点21×21的矩阵,并将这个矩阵里面的数据确定参考像元的个数N,若N>11,则将这些数据建立回归关系对目标数据进行拟合,否则将该目标数据判断为离群点(具体流程如图3所示)．若矩阵过大,则与目标像元的空间差异会增大，若矩阵过小,则参与修复的像元数达不到要求,这都将影响到数据修复的精度．

图2 参考像元的选取Fig.2 Reference pixel selection

I′t±1(x,y)=at±1×It±1(x,y)+bt±1,

(1)

其中:t表示待修复数据获取时间；t±1表示目标数据的参考数据；It±1(x,y)表示参考数据位于目标位置的臭氧柱总量;I′t±1(x,y)表示根据目标数据建立回归关系得到的拟合值;(x,y)表示目标像元所在的位置;a,b均是通过最小二乘法拟合所得到的值.

由于两个参考数据与目标数据所用的参考像元并不相同，即两个参考数据对目标像元有着不同的贡献度，因此目标像元的拟合值可以表示为

(2)

其中，ω1,ω2由地理加权回归获得:

(3)

图3 时空加权法流程Fig.3 Flow chart of spatial-temporal weighted fitting method

通过上述方法大致能修复90%左右数据,其余的缺失数据则可用t-1或t+1时刻该位置具有有效值的数据建立一组回归关系的拟合值代替目标像元的拟合值．

1.2.2 对流层臭氧残差法

对流层臭氧残差法的实质就是利用臭氧柱总量数据减去平流层臭氧数据,从而得到对流层臭氧总量数据．本文先根据NCEP/NCAR对流层顶高度数据和MLS臭氧廓线数据进行空间匹配,将臭氧廓线数据所在的经纬度与NCEP/NCAR对流层顶高度数据所在的格网进行匹配,从而得到该位置的对流层顶高度数据,再计算平流层臭氧柱总量(Stratosphere Column Ozone,SCO,其量值记为ηSCO),且由于MLS数据是臭氧体积比的形式,因此需要先将MLS数据与对流层臭氧数据统一为Dobson(DU)单位,因此其计算公式[27]可以表示为

(4)

其中，Vx为臭氧体积比的形式，P1为对流层顶的气压，0.215为平流层顶的气压值，单位为hPa，n为对流层顶到平流层顶气压带的条数，ΔP表示为两个相邻高度气压之差.之后统一SCO和臭氧柱总量(Column Ozone,CO,其量值记为ηCO)数据的空间尺度,将平流层的臭氧柱总量数据利用反距离加权法插值成与OMI臭氧柱总量对应的0.25°×0.25°的格网数据,再利用对流层臭氧残差法得到对流层臭氧柱总量(Troposphere Column Ozone,TCO,其量值记为ηTCO),具体流程如图4所示．

ηTCO=ηCO-ηSCO.

(5)

图4 对流层臭氧残差法Fig.4 Flow chart of tropospheric ozone residual method

2 中国区域对流层臭氧的反演

2.1 臭氧柱总量的修复试验

图5 2010年12月23日原始数据(缺失掩膜)Fig.5 Raw O3 data on Dec.23rd,2010 (missing mask)

本文所选择的研究区是根据中国区域经纬度范围所裁剪出的矩形区域(73.0°～ 135.5°E,3°～ 53.75°N),为了保证试验结果的可靠性,分别选取了2005年1月、4月、7月、10月的数据进行了4组试验进行精度比较．按照2010年12月23日(图5)的缺失掩膜进行修复,即利用该区域的缺失值来代替2005年4 d的缺失数据．其中,缺失数据的栅格点的个数为15 761个,缺失的经度范围最大达到了10°．

定量评价可以很好地评判数据的修复结果,表1所示的是对4组模拟试验利用反距离加权法、克里金插值法、时空加权法3种不同数据修复方法的均方根误差(RMSE)、平均绝对误差(MAE)的评价结果．反距离加权法和克里金插值法的评价结果具有一定的相似性,2种修复方法的RMSE和MAE的差距较小．时空加权法相较于前2种方法的RMSE和MAE均较小,说明其具有更好的修复效果．因此,说明时空拟合法具有更好的稳定性,修复数据的误差更小．

表1 2005年4组模拟数据修复方法的评价指标

图6—7、8—9、10—11、12—13所示的分别是对OMI臭氧柱总量数据修复结果和与之对应的3种不同修复方法与原始数据对比的散点图．通过2005年这4 d(图6—13)的数据可以发现中国区域的臭氧柱总量在200～460 DU之间,且高值区域集中于我国的东北部地区．同样,该数据也说明了臭氧柱总量在冬季、春季较高,而夏季、秋季较低的情况．随着纬度的上升,臭氧柱总量也在逐渐上升,而西部的青藏高原地区由于海拔较高,相较于同纬度的中东部地区,其臭氧总量相对较少．1月、4月的臭氧柱总量相对于7月、10月的臭氧柱总量在我国中南部地区差异较小,而在我国东三省1月和4月的臭氧柱总量明显高于7月和10月．

图6 2005年1月不同修复方法对比 a.原始数据;b.反距离加权法;c.克里金插值法;d.时空拟合法Fig.6 Performance comparison between different repair methods for data of January 2005 a.original data;b.inverse distance weighted interpolation;c.Kriging interpolation; d.the proposed spatial-temporal weighted fitting approach

图7 2005年1月不同修复方法散点图Fig.7 Scatter diagram of original data and repaired data for January 2005

图6a—d展示了2005年1月2日原始数据,以及通过掩膜裁剪后的利用反距离加权法、克里金插值法、时空加权法3种不同方法的修复结果．图6b和c分别是反距离加权和克里金插值的修复结果,可以看出,利用自身的空间相关性对臭氧总量数据进行修复在空间上出现了明显的不连续性,且对缺失区域的数据做了一定程度平滑处理．图6d是通过时空加权的方法实现对数据修复的,其修复结果相对于上述2种修复方法具有明显的优势,且修复结果与原始数据(图6a)基本一致,也没有出现不连续的现象．图7是3种不同修复方法对缺失区域的修复结果与原始数据的比较,反距离加权法相较于时空加权法和克里金插值法更为分散,说明其与真实值的偏离程度较大．因此,从定性和定量的角度考量,时空加权法相较于反距离加权法和克里金插值法具有更高的精度,修复效果最好．

图8 2005年4月不同修复方法对比 a.分别代表原始数据;b.反距离加权法;c.克里金插值法;d.时空拟合法Fig.8 Performance comparison between different repair methods for data of April 2005 a.original data;b.inverse distance weighted interpolation;c.Kriging interpolation; d.the proposed spatial-temporal weighted fitting approach

图9 2005年4月不同修复方法散点图Fig.9 Scatter diagram of original data and repaired data for April 2005

图10 2005年7月不同修复方法对比 a.分别代表原始数据;b.反距离加权法;c.克里金插值法;d.时空拟合法Fig.10 Performance comparison between different repair methods for data of July 2005 a.original data;b.inverse distance weighted interpolation;c.Kriging interpolation; d.the proposed spatial-temporal weighted fitting approach

图11 2005年7月不同修复方法散点图Fig.11 Scatter diagram of original data and repaired data for July 2005

图9、11、13分别是以2005年4月2日、7月2日、10月2日的数据作为原始数据,利用2010年12月23日的缺失区域作为掩膜的修复结果．通过图8—13可以发现,4月2日的修复结果除我国东北区域外未出现明显的空间不连续现象,而7月2日和10月2日的修复结果表明无论是反距离加权法还是克里金插值法在缺失的数据带上均存在明显的空间不连续现象,而时空加权法均能够保证数据在空间上的连续性,而且与原始数据具有相同的时空变化特征．通过散点图的修复结果与原始数据的对比可以发现,除4月时空拟合法相较于反距离加权法的数据较为分散以外,其余相较于反距离加权法和克里金插值法的离散程度更低,即与原始数据的偏差较小．虽然4月和7月的数据离散程度相对较大,但其与克里金插值相比离群点的个数较少,且偏差相对较小,说明了时空加权法相对于其余两种修复方法具有更好的效果．

图12 2005年10月不同修复方法对比 a.原始数据;b.反距离加权法;c.克里金插值法;d.时空拟合法Fig.12 Performance comparison between different repair methods for data of October 2005 a.original data;b.inverse distance weighted interpolation;c.Kriging interpolation; d.the proposed spatial-temporal weighted fitting approach

图13 2005年10月不同修复方法散点图Fig.13 Scatter diagram of original data and repaired data for October 2005

通过以上4组对比试验可以发现:从定性的角度考虑,时空加权法保持了数据在空间上的连续性,且与原始数据的偏离程度较小;从定量的角度考虑,无论是RMSE还是MAE的评价结果显示,时空加权法从时空变化的角度去考虑数据的修复,相较于反距离加权法和克里金插值法考虑到空间的因素的误差均较小,说明时空加权法对于OMI臭氧柱总量的修复效果较好．

图14为根据时空加权法对2019年臭氧柱总量修复的结果．总体上呈现出高纬度地区的臭氧柱总量高,低纬度地区少,西部地区相对于东部地区而言的臭氧柱总量较低．全年的臭氧柱总量位于220～400 DU之间,夏季最大值较低,而在中国区域内整体处于一个高值状态;冬季,我国南方地区的臭氧柱总量整体处于一个低值状态,北方的臭氧柱总量普遍偏高．

图14 2019年中国区域各月臭氧柱总量Fig.14 Monthly ozone column over China in 2019

2.2 对流层臭氧浓度的反演

图15为根据对流层臭氧残差法反演的2019年中国区域对流层臭氧的月际变化空间分布,可以看出除冬季(1、2、12月)以外,其余月份青藏高原地区的对流层臭氧相对于我国其他地方而言处于一个明显的低值区域,这主要是因为青藏高原地区的海拔较高,对流层的高度较低,因此其对流层臭氧柱量处于一个低值范围内．

图15 2019年对流层臭氧的月平均分布Fig.15 Monthly average distribution of tropospheric ozone in 2019

由图15可见,我国各地区的对流层臭氧柱总量差距较大,最高值与最低值相差100 DU左右,且总体上我国东部地区的对流层臭氧柱量低于中部地区和北部地区,这主要是受到地形因素的影响所致．由于气象条件和动力作用的影响,使我国北方地区的对流层臭氧柱总量常年处于较高水平;另一方面,南方经济发达地区主要为沿海城市,所产生的污染物能快速地被来自海上的大量的清洁空气所稀释．因此,同为经济发达地区的长三角和珠三角地区的对流层臭氧柱量明显低于京津冀地区．由图15可以看出,青藏高原范围内的臭氧柱总量处于一个明显的低值区域,这主要是由于青藏高原地区海拔较高,大部分地区的海拔高度都在4 000 m以上,臭氧柱总量的范围较窄;另一个原因是因为大气对流层与平流层之间的物质输送,稀释了对流层臭氧的含量,从而导致对流层臭氧柱总量的下降．

从臭氧的时间变化角度来说,夏季的臭氧柱总量较高,春秋季次之,冬季的对流层臭氧处于一个低值范围．这主要是因为夏季的温度高,加速了臭氧前体物的光化学反应,从而产生了较为严重的光化学污染．而东部沿海地区的臭氧柱总量低于中部地区,这主要是夏季风带来海上洁净的空气稀释了臭氧所致．

2.3 与OMI/MLS官方对流层臭氧产品的对比

将反演出的0.25°×0.25°对流层臭氧产品按照1°×1.25°的格网进行划分,从而求取每个格网的均值,再将其与OMI/MLS的对流层臭氧产品格网进行一一匹配,将拟合的结果与OMI/MLS官方臭氧产品进行对比(表2),其相关系数R基本在0.6以上,最高值达到了0.82,说明了对流层臭氧残差法拟合的结果与OMI/MLS的官方臭氧产品有着较高图16的相似性,反演的效果较好．

表2 2019年各月对流层臭氧残差法与 OMI/MLS官方产品的相关系数

为中国区域2019年夏季的OMI/MLS官方臭氧产品对流层臭氧的时空变化,可以看出,通过对流层臭氧残差法拟合的结果与OMI/MLS臭氧廓线数据的时空分布结果大概一致,但OMI/MLS臭氧廓线数据的官方产品空间分辨率较低,只针对月平均数据进行了反演,部分区域的反演结果存在着一定的数据缺失等问题．并且,OMI/MLS官方产品选用的是OMIO3PR数据进行反演,该廓线数据虽然与OMI数据属于同一种卫星产品,但其垂直分辨率较低,误差较大,从而导致反演的结果也存在着较大的误差．除此以外,OMI/MLS官方臭氧产品覆盖范围为60S°～60°N之间,而MLS臭氧廓线数据可以覆盖除南北极以外的所有区域．因此,首先利用时空拟合法对臭氧总量数据进行修复从而提高臭氧总量数据的精度,再将修复完整的OMIO3e臭氧数据和MLS臭氧廓线数据相结合,根据对流层臭氧残差法得到对流层臭氧数据集在一定程度上弥补了OMI/MLS官方臭氧产品在空间分辨率以及数据缺失等方面的缺陷,且可以及时监测对流层臭氧的日变化．

图16 2019年夏季(6月、7月、8月)OMI/MLS臭氧廓线数据Fig.16 OMI/MLS ozone profile data in June (left),July (middle),and August (right) of 2019

3 结论与展望

运用时空拟合的方法对OMI缺失数据进行修复,并与反距离加权法、克里金等插值方法从定性和定量的角度进行对比分析,再根据对流层臭氧残差法反演中国区域的对流层臭氧柱量,得到的结论主要有:

1)对于不同日期的对比试验,时空拟合法相对于克里金和反距离加权的修复方法从定量的角度其RMSE和MAE均较小,从定性的角度也更为平滑,没有了条带性的差异．

2)运用对流层臭氧残差法构建的对流层臭氧估算模型相对于其他的对流层臭氧的反演方法具有一定的优势,且比较适合中国区域的反演,其反演的结果与OMI/MLS臭氧廓线数据有着较高的相关性．由于夏季的温度较高,光化学反应加剧了臭氧前体物的转化,从而导致夏季的对流层臭氧浓度较高．

3)对流层臭氧残差法反演的结果与OMI/MLS官方臭氧产品相比较为相近,其相关系数R均在0.6以上,最高可达0.82,说明反演的结果具有一定的代表性．

本文所选择的时空拟合法以及对流层臭氧残差法虽然在一定程度上能够得到较为精确的信息,但随着数据缺失条带的增大,时空拟合法对臭氧总量数据的修复效果也越来越差．另外,MLS数据虽然具有更高的精度,但与OMI传感器测量臭氧总量的原理不尽相同,因此还需对臭氧廓线数据进行系统的误差校正,以此来进一步提高对流层臭氧浓度反演的精度．