基于二次多项式的AOD数据融合方法

李广超，李如仁，赵阳阳，杨 震，卢月明

(1. 沈阳建筑大学，辽宁 沈阳 110168； 2. 中国测绘科学研究院，北京 100830； 3. 中国矿业大学，北京 100083)

气溶胶光学厚度(aerosol optical depth，AOD)作为大气气溶胶最重要的参数之一，是表征大气混浊度或大气中气溶胶总含量的重要物理量[1]；反映了整层大气或某特定大气层中颗粒物对太阳辐射的削弱程度，是评价大气颗粒污染和气溶胶辐射气候效应的一个关键参数[2]；它可以用来推算大气气溶胶含量，是确定大气气溶胶辐射气候效应的关键因子，气溶胶光学厚度也是大气污染物浓度的重要指示因子，其与PM10、NO2、SO2浓度之间具有较好的关系[3]。利用AOD推算大气气溶胶含量的研究已经比较成熟，全球大气气溶胶类型和含量变化与气候变化和大气环境污染密切相关，是气象学、环境学和医学研究关注的热点问题，严重影响着人们的生活[4-5]。

在气溶胶融合和卫星反演产品的验证等方面研究，线性回归分析应用广泛，许多学者都在利用线性回归拟合方法对大气气溶胶数据进行相关的处理研究，包括Aqua和Terra MODIS AOD数据拟合，这种方法在气溶胶方面的应用已经相当成熟，并取得了相对较好的成果。Liang等对美国华盛顿特区的4个地点，中国北京市、非洲尼日利亚和加拿大布拉茨湖利用新的AOD估计算法对MODIS AOD进行估计，并利用MODIS AOD估计值和AERONET AOD进行线性回归分析，取得了比较好的效果[6]；Lary等利用机器学习算法对Terra/Aqua MODIS AOD与AERONET AOD进行线性回归分析，有效提高了线性相关系数[7]；王伟齐等利用Terra卫星和Aqua卫星的MODIS 3 km气溶胶光学厚度产品，分析了地面PM2.5和两颗卫星AOD的时空分布特征，并在时空匹配的基础上，建立了AOD与PM2.5浓度之间的线性回归模型[8]；李龙在有效值位置处的MODIS/Terra AOD和MODIS/Aqua AOD，建立二者之间的线性关系模型[9]；，闫婧华和Qi等分别利用MODIS和MISR在550 nm通道观测的气溶胶光学厚度与AERONET观测数据匹配后进行线性回归分析[10-11]。

针对Aqua和Terra MODIS AOD数据研究方面，采用线性回归方法拟合数据结果精度不够高。朱绪超等利用Landsat 8 OLI遥感影像对塔里木河下游地区进行植被指数提取，并与野外实测的叶面积指数值进行回归分析，并利用效果最好的模型进行模拟和精度的验证，结果表明：在3种样本分析方式中，使用的6种植被指数与LAI相关性均很好，说明在塔里木河下游使用遥感影像提取植被指数来反演区域叶面积指数的方法是可行的，3种分析方式中相关性最好的均是二次多项式模型[12]。刘明亮等采用二次多项式拟合外方位元素模型，与国外的定向片法模型进行平差算法对比，结果表明基于二次多项式拟合外方位元素模型的平差方法有效地提高了影像的平面精度与高程精度，与定向片法模型相比，该算法简单，结果同样满足1∶5万摄影测量的要求[13]。因此，本文提出二次多项式回归算法对Aqua和Terra MODIS AOD数据进行拟合，利用二次多项式回归和线性回归算法拟合同一组MODIS AOD数据，采用均方根误差(mean absolute error，MAE)，平均绝对误差(root mean square error，RMSE)，线性相关系数(Pearson)R作为性能评价指标[9]。

1 数据准备和处理

本文所用的MODIS数据来自美国戈达德宇宙飞行中心数据分发系统(http:∥ladsweb.nascom.nasa.gov/data/search.html)。MODIS传感器搭载于美国国家航空航天局(National Aeronautics and Space Administration，NASA)地球观测系统(earth observing system，EOS)中Terra和Aqua卫星上[14]。Terra卫星和Aqua卫星分别于1999年12月和2002年5月发射升空。Terra卫星于当地时间上午10:30左右自北向南通过赤道，Aqua卫星则是于当地时间下午1:30左右自南向北通过赤道[15]。Terra和Aqua两颗卫星分别提供了MOD04_3K和MYD04_3K两种MODIS AOD产品数据。研究区域选用的地理范围为(35.5°N—43°N，113°E—120°E)的京津冀地区，MODIS AOD数据是2015年随机选择的一天和4—6月数据作为研究对象，频率为每天，空间分辨率为3 km。本文采用MODIS Collection 6 MOD04_3K和MYD04_3K数据集中参数名为“Optical_Depth_Land_And_Ocean”，波段为550 nm的2级AOD数据。

本研究中利用C#、ArcGIS Engine、Visual Studio 2013程序对MOD04_3K和MYD04_3K原始图像进行批量处理，通过重采样获取MODIS AOD值，并采用java语言程序对AOD数据进行相关处理，最终使用Python语言来实现数据的拟合。

2 方法原理

本研究中分别使用线性回归、二次多项式回归算法对数据进行拟合处理，二次多项式函数拟合的基本原理[16]如下。

设拟合二次多项式函数为

(1)

式中，a0、a1、a2为拟合系数；xi表示MOD04_3K的值；yi表示MYD04_3K的值。

各点到这条曲线的距离之和，即偏差平方和为

(2)

式中，R2为偏差平方和。

为了求得符合条件的ai值，对式(2)右边求ai偏导数，因而得到

(3)

(4)

(5)

将式(3)—式(5)左边进行化简，可以得到

(6)

(7)

(8)

把式(6)—式(8)表示成矩阵的形式，可以得到

(9)

将式(9)简化后可得到

(10)

也就是说XA=Y，那么A=(XTX)-1XTY，从而求出ai值。

对于二次多项式，则至少需要6个以上已知点来确定模型中的参数[17]。

3 MODIS AOD数据拟合结果与分析

为了测试二次多项式数据融合方法的有效性和避免拟合结果的偶然性，本文以线性回归作为对比方法，对一天和4—6月的Aqua和Terra MODIS AOD数据进行融合。试验方法使用Python语言实现，采用MAE、RMSE、R作为性能评价指标。

3.1 MODIS AOD数据的拟合

本文对随机选择一天和4—6月的Aqua和Terra MODIS AOD数据进行拟合分析。以MOD04_3K为X轴、MYD04_3K为Y轴，用线性回归和二次多项式回归算法拟合MODIS AOD数据，拟合结果及拟合方程如图1所示。

注：图中上方实线为y=x线，下方实线为拟合方程线。图1 不同卫星MODIS传感器上AOD数据相关性

3.2 拟合结果及讨论

试验结果得出的RMSE、MAE、R、二次多项式回归拟合方法比线性回归拟合方法RMSE提升百分率、MAE提升百分率、R提升百分率，见表1。

表1 回归分析得出的RMSE、MAE、R的值RMSE提升百分率、MAE提升百分率、R提升百分率

通过利用线性回归和二次多项式回归算法对Aqua和Terra MODIS AOD数据拟合结果可以看出，无论是对于一天或者是3个月的Aqua和Terra MODIS AOD数据，二次多项式回归算法拟合得出RMSE、MAE、R值都要优于线性回归拟合方法得出的值。对一天MODIS AOD数据的回归拟合中，二次多项式回归拟合得出的RMSE、MAE、R值分别为0.144 2、0.098 9、0.918 9；线性回归拟合得出的RMSE、MAE、R值分别为0.159 8、0.113 8、0.899 5；二次多项式回归拟合方法比线性回归拟合方法RMSE提升9.76%、MAE提升13.09%、R提升2.16%。对3个月MODIS AOD数据的回归拟合中，二次多项式回归拟合得出的RMSE、MAE、R值分别为0.191 2、0.126 9、0.913 1；线性回归拟合得出的RMSE、MAE、R值分别为0.202 2、0.135 8、0.896 7；二次多项式回归拟合方法比线性回归拟合方法RMSE提升5.44%、MAE提升6.55%、R提升1.83%。

4 结 论

本文运用MOD04_3K(MODIS/Terra)和MYD04_3K(MODIS/Aqua)两种MODIS AOD产品数据，使用二次多项式回归拟合方法和线性回归拟合方法对京津冀地区的AOD卫星产品开展了对比研究，并通过平均绝对误差、均方根误差、线性相关系数来对拟合结果来进行评价。对随机选择一天和4—6月的MODIS AOD数据进行拟合，结果显示：二次多项式回归拟合方法得出的RMSE、MAE、R值优于线性回归拟合方法得出的值。

研究结果表明：二次多项式回归拟合方法不仅在Aqua和Terra MODIS AOD数据拟合研究方面适用，而且对于同一组Aqua和Terra MODIS AOD数据，在一天数据的拟合试验中，二次多项式回归拟合方法比线性回归拟合方法得到的RMSE提升9.76%、MAE提升13.09%、R值提升2.16%；在3个月数据的拟合试验中，二次多项式回归拟合方法比线性回归拟合方法得到的RMSE提升5.44%、MAE提升6.55%、R值提升1.83%，因此利用二次多项式回归拟合Aqua和Terra MODIS AOD数据能提高MODIS AOD数据拟合的精度。

在后续的研究中，将进一步加入其他考虑因子，如湿度、温度、降雨之类的气象因子和人类活动、地表覆盖等因子，进而对AOD数据精度进行修正，得到高精度的AOD数据。

[1] 刘浩，高小明，谢志英,等.京津冀晋鲁区域气溶胶光学厚度的时空特征[J].环境科学学报，2015，35(5):1506-1511.

[2] 赵胡笳.中国东北城市地区大气气溶胶光学特性及其直接辐射效应研究[D].北京：中国气象科学研究院,2014.

[3] 蒋维东.AOD卫星观测资料三维变分同化试验研究[D].南京：南京信息工程大学,2015.

[4] 李晓静，高玲，张兴赢,等.卫星遥感监测全球大气气溶胶光学厚度变化[J].科技导报，2015，33(17):30-40.

[5] 罗云峰，吕达仁，周秀骥,等.30年来我国大气气溶胶光学厚度平均分布特征分析[J].大气科学，2002，26(6):721-730.

[6] LIANG S，ZHONG B，FANG H.Improved Estimation of Aerosol Optical Depth from MODIS Imagery over Land Surfaces[J].Remote Sensing of Environment，2006，104(4):416-425.

[7] LARY D J，REMER L A，MACNEILL D，et al.Machine Learning and Bias Correction of MODIS Aerosol Optical Depth[J]. IEEE Geoscience & Remote Sensing Letters，2009，6(4):694-698.

[8] 王伟齐,张成网,臧增亮,等.Terra和Aqua卫星MODIS 3 km AOD与北京PM_(2.5)对比分析[J].气象科学,2017,37(1):93-100.

[9] 李龙.基于泛克里金方法的气溶胶光学厚度融合[D].上海：华东师范大学，2015.

[10] QI Y L，GE J M，HUANG J P.Spatial and Temporal Distribution of MODIS and MISR Aerosol Optical Depth over Northern China and Comparison with AERONET[J].Chinese Science Bulletin，2013，58(20):2497-2506.

[11] 闫婧华.MISR和MODIS气溶胶光学厚度产品的数据融合及其在大气细颗粒物监测中的应用研究[D].青岛：中国海洋大学，2014.

[12] 朱绪超，袁国富，易小波,等.基于Landsat 8 OLI影像的塔里木河下游河岸林叶面积指数反演[J].干旱区地理(汉文版)，2014，37(6):1248-1256.

[13] 刘明亮，康倩.二次多项式拟合外方位元素卫星三线阵CCD影像平差[J].航天返回与遥感，2011，32(3):77-82.

[14] 马宗伟.基于卫星遥感的我国PM2.5时空分布研究[D].南京：南京大学，2015.

[15] 陈韵竹.气溶胶光学厚度的数据融合研究[D].上海：华东师范大学,2013.

[16] 黄征凯.确定区域似大地水准面方法的比较分析[D].桂林：桂林理工大学，2013.

[17] 程跃辉，王茂和.用二次多项式拟合曲线推算最终沉降量[J].华东公路，1999(6):46-49.