基于二次多項式的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.