采用Landsat8產品算法流程的高分一號數據大氣校正

張曉月，李琳琳，王 瑩，張 琪，李國春

采用Landsat8產品算法流程的高分一號數據大氣校正

張曉月1,2,3，李琳琳2,3，王 瑩2,3，張 琪2,3，李國春4

（1.中國氣象局沈陽大氣環境研究所，沈陽 110166； 2. 遼寧省生態氣象和衛星遙感中心，沈陽 110166； 3.高分辨率對地觀測系統遼寧數據與應用中心，沈陽 110166； 4.沈陽農業大學農學院，沈陽 110866）

高分一號（GF-1）衛星搭載的傳感器，實現了高分辨率和寬幅成像能力的結合，使其在精準農業方面應用發揮重要作用。該文在6S大氣輻射模擬模型基礎上，參照LaSRC大氣校正流程，設計了GF-1衛星WFV/MSS數據從算法原理分析到編碼實現的大氣校正。算法應用大氣總傳輸率、水汽透過率和大氣后向半球反照率等參數和高程、大氣可降水以及臭氧含量等值，循環計算使GF-1像元紅藍通道值比等于來自MODIS數據的紅藍通道值比，求得GF-1像元氣溶膠，再將包含當日大氣可降水和臭氧含量的輔助數據文件等代入6S模型，得到大氣校正后的地表反射率。試驗表明，該大氣校正方法在有農田和林木等植被覆蓋的中低緯度大氣校正效果較好，對稀疏植被的荒漠裸地和建筑地表大氣校正效果相對稍差。比較GF-1 WFV/MSS數據與Landsat8(LC8) OLI數據的基于6S模型LaSRC流程算法的大氣校正結果，GF-1 WFV/MSS各傳感器與LC8 OLI大氣校正結果的相關系數為0.825～0.972，2種衛星數據大氣校正結果相關性高，其中WFV相較于MSS顯示出與LC8 OLI更相近的大氣校正結果。結果表明，應用6S模型原理參照LaSRC校正流程設計的自行估計數據逐像元水平氣溶膠參數的GF-1衛星數據大氣校正方法應用方便、可操作性強，適合生長季農林監測等陸面應用。

衛星；遙感；GF1 WFV/MSS；LC8 OLI；6S；LaSRC；大氣校正

0 引 言

高分一號衛星具備時間分辨率和空間分辨率雙高的特性，因其軌道穩定性和數據可持續性等諸多優勢在多領域取得了較好應用效果[1]。高分一號數據不僅局限在對地表地貌的觀測，而是越來越多用于定量應用遙感，農業方面主要在農作物種植區劃、農業用地精準提取、農作物分類識別、農田面積估算、農作物產量預報等領域發揮著重要作用[2]。但眾所周知衛星傳感器接收到的地面輻射受大氣分子、氣溶膠以及臭氧、水汽等氣體的散射、吸收影響會造成信號失真，因此在應用高分一號衛星遙感數據進行地表參數定量反演時，應盡可能將大氣噪聲從地面輻射信息中剝離，保留地面物體真實反射率信息。大氣校正是解決這一問題的重要手段，按照是否具備明確物理意義分為相對校正模型和絕對校正模型[3]。目前較為常用的大氣校正方法有圖像特征法、經驗線性法和大氣輻射傳輸理論法[4]。圖像特征法是一種相對校正方法，針對圖像特征消除大氣影響，常見的有直方圖均衡法[5]、暗目標法[6-7]、不變目標法[8]、內部平均法（IARR，internal average relative reflectance）[9]、平場域法（FF，flat field）[10]等。經驗線性法結合野外試驗，依賴于地面同步定標點實測數據[11]。基于輻射傳輸理論（RT，radiative transfer）的大氣校正方法具備明確的物理意義，該類方法以輻射傳輸理論為基礎，考慮大氣參數等影響，校正精度較高、應用范圍廣泛。

目前大氣校正處理應用較為廣泛的6S模型是美國馬里蘭大學地理系在5S[12]模型基礎上研發的[13]，該模型在氣體透過率、瑞利散射及氣溶膠光學厚度計算等方面進行了改進，在MODIS產品大氣校正計算方面表現尤為突出[14]。2005年6SV1版本公布，該版本較之前的模型主要增加了輻射極化的計算[15]，2015年更新的6SV2.1是最新版本。學者們應用6S模型進行了多種遙感數據的大氣校正研究，推進了模型的改進和應用。劉佳等[16]研究了基于6S模型的高分一號衛星16 m分辨率數據大氣校正，并將批量處理的效率提高了75.0%以上。孫林等[17]提出了地表反射率支持的GF-1 PMS氣溶膠光學厚度反演及大氣校正方法，通過對北京、太湖2個AERONET站點觀測的氣溶膠光學厚度驗證，得到反演地表反射率與實測地表反射率誤差低于0.015的結果，算法具備較高精度。孫元亨等[18]應用6S模型對GF-1 WFV數據進行了大氣校正，開展了GF-4 PMS與GF-1 WFV地表反射率及NDVI一致性的分析，校正前后相關系數分別為0.74和0.77，二者表現出較好的一致性。薛現光等[19]提出在6S模型基礎上，以氣溶膠厚度為自變量對大氣參數進行擬合的算法，對GF-1進行大氣校正，與FLAASH對比具有較好的一致性。胡勇等[20]利用MODIS大氣產品和6S模型對GF-1 WFV數據進行了大氣校正處理，不同地表覆蓋類型GF-1 WFV和MODIS反射率數據的交叉對比結果顯示出較好的一致性。

基于輻射傳輸理論基礎，近年來相關機構相繼研發了一系列大氣校正流程化平臺和系統，使大氣校正處理流程更便捷，科研人員利用這些算法軟件開展了一系列大氣校正相關研究。LOWTRAN模型1971年由美國空軍劍橋研究所研制，已發展至LOWTRAN7版本[21]。MODTRAN模型從LOWTRAN 模型發展而來，由美國空軍地球物理實驗室研制，將光譜分辨率由LOWTRAN的20提高到2 cm-1[22-23]。楊貴軍等[24]利用MODTRAN4和MODIS產品，對環境與減災小衛星HSI高光譜數據進行了大氣校正。鄭盛等[25]利用MODTRAN4模型對HJ-1衛星CCD數據進行逐像元大氣校正，并從光譜曲線、MODIS地表反射率產品、NDVI的3個方面比較了大氣校正效果。熊世為等[26]利用MODTRAN模型對HJ/CCD進行了大氣校正，校正結果與高精度MODIS地表反射率產品在植被、居民地、水體3類地物一致性較高。美國空軍實驗室1998年開始以MODTRAN為基礎發起了FLAASH的研發[27]，專業遙感數據處理軟件ENVI大氣校正模塊集成了FLAASH模型[28]，李衛國等[29]利用FLAASH進行了基于薄云霧去除的ETM+數據大氣校正。ATCOR[30]模型是以MODTRAN4為基礎開發的大氣校正模型，ERDAS軟件集成了ATCOR模型。常用的還有商用大氣校正軟件包ACORN[31]主要針對0.4～2.5m高光譜和多光譜數據進行大氣校正。

本文分析了基于6S模型和LaSRC（Landsat Surface Reflectance Code）流程進行GF-1數據大氣校正的適用性，進行了通道相關等改進處理和重新規劃，設計實現了針對GF-1數據的大氣校正算法和軟件平臺。LaSRC方法依賴于外部輔助數據，需以待校正數據同日的臭氧和水汽等大氣影響氣體產品作為初值，考慮到對時效性的要求和實時資料有時難以獲得，平臺提供了利用近6 a臭氧和水汽逐日值作為替代和使用實時值的2種方案。從而保證在無外源數據支持下仍然能夠提供大氣校正產品。鑒于中國境內暫無LaSRC產品，且LaSRC移植到RSD（Remote Sensing Desktop）平臺處理LC8 L1T數據能夠代表原LaSCR產品，本文以同地同時次LC8 OLI L1T的RSD大氣校正產品和GF1 WFV/MSS數據大氣校正處理結果對比分析二者大氣校正效果和精度差異，探索對國產高分衛星數據的高精度大氣校正處理平臺軟件工程化實現的實用方法。

1 原理方法

6SV模型是基于6S大氣輻射傳輸模型基礎上發展的極化矢量版本，美國USGS EROS項目針對Landsat8（LC8）數據開發了LaSRC大氣校正程序。該程序運行在Linux操作系統，用于美國境內LC8數據地表反射率產品反演，對于光譜波段較窄的Operational Land Imager(OLI)傳感器具有較好的適用性[32]。應用待校正LC8 OLI數據反演像元水平氣溶膠，以及水汽、臭氧含量從MODIS產品或NCEP數據獲取是LaSRC流程算法求算地表反射率的2個主要特點。以該流程對LC8 OLI數據進行大氣校正時，首先利用查找表和大氣總傳輸率、氣體透過率、大氣內反射率等6S模型參數從預定義的22個氣溶膠光學厚度查找表估算地表反射率初值；從MODIS產品提取紅藍通道地表反射率比值，根據該比值與NDVImir線性關系，由LC8 OLI地表反射率初值和線性斜率及截距估算LC8 OLI紅藍通道比，在3種埃氏系數等級和22個氣溶膠光學厚度下循環計算，找到使MODIS紅藍通道比和LC8 OLI紅藍通道比殘差最小時的反射率，求得氣溶膠；再利用MODIS產品通過空間插值得到LC8 OLI像元水平水汽、臭氧含量，利用DEM數據及氣壓與海拔高度關系計算并空間插值得到LC8 OLI像元水平大氣壓強，通過衛星參數文件得到像元幾何參數；將表觀反射率，氣溶膠以及相關參數輸入6S模型最終求得地表反射率（圖1）。

圖1 算法原理流程圖

1.1 6S模型算法

假設地表為均勻朗伯面，大氣中分子與氣溶膠散射及氣體吸收水平一致，衛星傳感器接收到的輻射由大氣路徑輻射、透過大氣的目標地物直接反射輻射、經由大氣散射的目標地物反射輻射、背景地物散射輻射以及背景地物多次地表散射輻射構成，衛星傳感器接收到的反射率表示為[33]

式中TOA為表觀反射率；為大氣內反射率，為地表反射率，為大氣后向半球反照率，()為下行輻射傳輸率，()為上行輻射傳輸率，為太陽天頂角，為傳感器天頂角。

大氣內反射率又稱大氣路徑輻射，大氣中同時存在大氣分子及氣溶膠散射和氣體吸收，簡單將大氣分子散射和氣體吸收的作用隔離開來分別考慮顯然與實際情況不符，則假設存在一種平均狀況，為便于計算認為有一半的水汽混合在氣溶膠中參與散射，綜合考慮大氣內各分子及氣體耦合作用，將大氣內反射率計算如下[34]

式中為瑞利散射反射率，Tg/2為半程水汽透過率。將式（1）帶入表達式并改進如下：

式中為大氣透過率，Tg為臭氧和其他氣體的透過率。上述式中，反射率、傳輸率、透射率無量綱，角度單位rad。

將式（4）帶入式（2），即可進行地表反射率的求算，地氣解耦是實現大氣校正或氣溶膠光學厚度反演的重要一環。應用該方法對LC8 OLI數據進行大氣校正時，氣溶膠的模擬由預先建立的查找表及MODIS產品估算。

1.2 光譜響應分析

從美國USGS網站及中國資源衛星應用中心網站下載LC8 OLI、GF1 WFV/MSS光譜響應數據，制作LC8 OLI、GF1 WFV/MSS可見光及近紅外波段光譜響應曲線（圖2），用于計算初始的分子光學厚度參數、臭氧傳輸參數、水汽傳輸的兩組參數和其他氣體的3組參數。

圖2 LC8 OLI、GF-1 WFV/MSS光譜響應曲線

比較二者光譜響應差異，GF-1的4個WFV傳感器及2個MSS傳感器光譜響應曲線相近，LC8 OLI紅光及近紅外波段波譜范圍相對較窄，藍綠波段光譜響應函數與GF-1 WFV/ MSS差異較小，雖然Landsat8 OLI與 GF-1 WFV/MSS對應通道并非完全相同，但波段設置較為相似的設計，使應用LC8 OLI適用的大氣校正算法進行GF1 WFV/MSS數據大氣校正成為較好的選擇。

應用積分算法對各傳感器不同波段光譜響應數據進行運算（式5），通過建立以LC8 OLI為基準的光譜匹配系數（式6），與USGS典型地物光譜庫數據進行對比，修訂GF-1 WFV/MSS光譜響應函數差異[35]。

式中1和2為波段范圍最小值和最大值，()為波長處光譜響應值，()為波長處反射率。計算Landsat8,i和GF1,i，并進行相應波段的比較求得光譜匹配系數。

1.3 像元幾何參數

本文進行大氣校正時，需要求算每個像元衛星方位角和天頂角、太陽方位角和天頂角。利用隨衛星元數據發放的參數文件中給定的仰俯、側滾、偏航等角度及地理定位等信息計算，對于一條掃描線任意列數據衛星方位角和天頂角與給定參數關系如下：

按照衛星過境時間及像元經緯度進行每一像元太陽方位角和天頂角的計算，需對影像進行精準的幾何（正射）校正，理論上認為對影像的幾何定位與大氣校正無既定處理順序[36]。本文首先對GF-1進行正射校正，根據校正后的幅寬與日地距離比得到整幅的太陽天頂角變化，然后根據幅寬像元數線性插值，得到逐像元太陽天頂角。太陽方位角對計算反射率影響較小且變化很小，這里使用了景中心數值。

1.4 表觀反射率

獲取像元表觀反射率首先需通過輻射定標系數將像元通道計數值轉換為像元表觀輻亮度，再由大氣層頂太陽輻照度、日地天文距離、太陽天頂角數據將大氣層頂表觀輻亮度轉換為大氣層頂表觀反射率，計算公式如下：

式中為表觀輻亮度，W/(m2·sr·m)；DN為像元通道計數值；scale和offset為輻射定標系數；為日地天文距離；為太陽輻照度，W/(m2·sr·m)；為太陽天頂角，rad。

1.5 大氣參數

1.5.1 氣溶膠光學厚度

地表反射率與氣溶膠的耦合關系使得在進行大氣校正時，氣溶膠的剝離存在一定難度，在一些大氣校正平臺中往往是根據觀測點地理位置及時間選擇一個氣溶膠模式，如限定緯度、季節和區域類型等。這種預定義的氣溶膠模式較粗糙，即使計算當日氣溶膠也多是整幅的氣溶膠而不是逐像元的氣溶膠值，代表性較差。

本文以LaSRC流程對GF-1數據進行大氣校正時，氣溶膠的估算采用引入MODIS產品為輔助數據的方法，首先計算MODIS產品紅藍波段地表反射率比值，根據紅藍波段比與NDVImir線性關系，以及NDVImir與NDVI的線性關系[36]，估算GF-1 WFV/MSS紅藍波段地表反射率比值，分別比較埃氏系數為1.0、1.75、2.5這3種情況和22個氣溶膠光學厚度時MODIS和GF-1 WFV/MSS數據紅藍波段地表反射率比值的殘差，由殘差最小時地表反射率求得氣溶膠。

以2016年8月26日遼寧境內條帶號121、行編號31的LC8數據為試驗，分別采集有植被和稀疏植被覆蓋地區樣本各300個，比較NDVI與NDVImir關系。

圖3 NDVI與NDVImir比較分析

圖3也可看出NDVI與NDVImir線性關系明顯，圖中有植被覆蓋地區NDVI與NDVImir數值均較稀疏植被覆蓋地區大，有植被覆蓋地區二者相關系數為0.958，線性回歸方程常數項和一次項分別為?0.262和1.097；稀疏植被覆蓋地區二者相關系數為0.912，線性回歸方程常數項和一次項分別為?0.284和1.128；2個方程均通過0.05水平顯著性檢驗。

由于計算輸出結果為波長550 nm處AOT，需根據AOT與波長之間的關系[37]求算其他各波段AOT數值，公式如下：

1.5.2 瑞利光學厚度

在獲取瑞利散射反射率時需以瑞利光學厚度（ROD）為輸入項，ROD是波長及海拔高度的函數[38-39]，波長一定時海平面處ROD計算如下

式中()為波長處ROD。為海拔高度，m；對于任意海拔高度處ROD計算為

1.5.3 大氣影響氣體含量

應用MODIS產品獲取逐日大氣水汽、臭氧含量，該數據空間分辨率5.5 km，將其與待校正數據進行像元尺度的空間匹配。當MODIS產品實時性受限，在所需處理時段數據無逐日MODIS產品數據時，對于臭氧、大氣可降水氣候值的計算采用2013－2018年近6 a MODIS產品逐日動態滾動平均值進行補充，雖然不是待校正當日數據，但相較于根據經緯度及時間的查找表參數確定方式仍具有較高精度。當平臺未檢測到與待校正數據同日MODIS產品輔助數據時，提示并自動調用動態平均值作為輸入量完成大氣校正處理。

1.5.4 大氣壓強

大氣壓強數據用于計算大氣后向半球反照率、氣體透過率、大氣總傳輸率、大氣內反射率等參數，影像各像元處大氣壓強通過海平面氣壓及海拔高度關系計算

式中為大氣壓強，hPa；0為海平面氣壓，hPa。

2 數據及處理

2.1 試驗數據

LC8 OLI數據通過美國USGS網站下載，數據級別為L1TP。GF-1 WFV/MSS數據來源于中國資源衛星應用中心網站，數據級別L1A。篩選2016—2017年GF-1數據，在比較GF-1與LC8大氣校正結果差異時，選取二者過境時傳感器采集數據有較大交叉覆蓋范圍的同日數據，LC8與GF-1均在地方時午前過境，選取同日數據大氣狀況差異小，便于對二者的大氣校正結果進行比較分析（見表1）。

表1 衛星數據集

2.2 數據處理

對于LC8 OLI數據大氣校正流程的實現，首先使用數據集提供的質量評價數據對遙感數據進行云區、云陰影、水體的檢驗，對識別出的像元缺漏和云等無法處理的數據進行標記隔離。通過屬性文件獲取定標信息進行輻射定標、獲取太陽角度進行太陽高度訂正，計算出表觀反射率。利用查找表估算初始地表反射率，通過NDVImir估算紅藍波段地表反射率比值，與MODIS產品紅藍波段地表反射率比值對比求算AOT，最后結合表觀反射率和大氣氣體含量以及其他參數完成大氣校正。對于GF-1 WFV/MSS數據大氣校正流程的實現，與LC8 OLI不同的是，GF-1 WFV/MSS為未經正射校正的L1A級數據，使用數據集提供的RPC參數結合全球30 m分辨率DEM數據進行正射校正，使二者進行精確幾何配準。其次GF-1 WFV/MSS的定標系數由中國資源衛星應用中心每年更新一次，數據期間可能有變化，應用對應的LC8 OLI數據對定標系數進行交叉驗證，提高定標的準確性。同時GF-1 WFV/MSS需要對側滾、仰俯和偏航角處理得到相關角度數據（式7、8）。經上述處理后，計算得到GF-1 WFV/MSS的表觀反射率，再通過NDVI求算紅藍波段地表反射率比值，完成后續大氣校正處理。最后將LC8 OLI和GF-1 WFV/MSS數據進行相同區域的裁剪和配準（圖4）。

圖4 數據處理主要流程

3 結果與分析

3.1 LC8 OLI大氣校正

LaSRC的開發是為進行LC8 L1G產品大氣校正的，目前USGS共享的數據為L1T級別，由于不同級別數據的像元重定位，使校正結果存在少量偏差。利用同一時次LC8 OLI數據以及其地表反射率產品對算法實現精度及適用性進行分析。從USGS網站下載同時次LC8 OLI數據（數據集名稱為：LC08_L1TP_030031_20170831_20180125_01_T1）和地表反射率產品（數據集名稱為：LC08_CU_015008_20170831_20180126_C01_V01_SR），對LC8 OIL數據進行大氣校正處理并與USGS LC8 OLI地表反射率產品逐波段比較（圖5）。

從LC8 OLI數據大氣校正結果與USGS LC8 OLI地表反射率產品中以100×100個像元的窗口采集樣本點各2 500個，比較二者相關性。從圖4中可以看出應用LC8 OIL數據的大氣校正結果與USGS LC8 OLI地表反射率產品可見光至短波紅外的第1～第7個波段相關關系均較好，決定系數為0.993～0.998，對LC8 OIL數據進行大氣校正算法與美國USGS發布的LC8 OLI地表反射率產品擬合度高。本地化LaSRC方法對于LC8 L1TP級別數據的大氣校正能夠代表USGS LC8 OLI地表反射率產品水平。

3.2 GF1 WFV/MSS大氣校正分析

對2016－2017年GF1 WFV/MSS試驗數據進行大氣校正處理。圖6為校正前后2017-09-07日遼寧境內GF-1 WFV/MSS數據1、2、3波段合成真彩色圖像，分析大氣校正處理對影像的目視清晰度的影響。

注：λ為波長，n為樣本容量，R2為決定系數。

圖6 GF-1大氣校正前后影像

從圖6中可以看出，由于大氣校正過程對大氣分子、氣溶膠、水汽、臭氧等噪聲的處理，大氣校正后影像清晰度提升，尤其是影像下部受大氣影響明顯區域，大氣校正后目視判讀性得到改善，這一效果是相對大氣校正無法實現的。

應用下墊面地類較為豐富的2017年6月5日GF-1 MSS數據為例，該數據下墊面類型包括：裸土、農田、林地、建筑、淡水、海水，且該時段農林地表植被生長利于進行大氣校正前后的地表分析，數據景中心大氣層臭氧質量濃度0.642 mg/L，水汽質量濃度0.96 mg/L，氣溶膠光學厚度0.35。分析大氣校正對通道數據的影響，從該數據提取大氣校正前后所有像元，計算逐通道每景影像所有像元混合的表觀反射率均值與地表反射率均值的差值，定義為校正量（以符號表示）。

從表2看出，對于GF-1 MSS傳感器的藍、綠、紅通道，表觀反射率與地表反射率差值逐漸減小，尤其藍通道降低最明顯、綠通道次之、紅通道降低最少；相較于可見光波段，近紅外波段大氣校正前后地表反射率表現恰好相反，各傳感器近紅外波段地表反射率較表觀反射率略有增加，且部分數值增加較多，算法適應性還需進一步研究。大氣校正處理對藍、綠通道影響較大，這與波長較短的藍、綠光易受大氣影響，長波輻射對大氣的穿透性較強的規律相一致。對于受大氣影響較大的藍、綠通道，大氣校正處理使各傳感器地表反射率相較于表觀反射率值降低了6.12%～8.2%和2.6%～4.29%，紅波段校正量為0.1%～2.18%，由此可見，數據定量應用時缺乏精確的大氣校正處理必將引起相應后續誤差。

表2 GF-1 WFV/MSS校正量

從2017年6月5日GF-1 MSS數據中采集裸土、農田、林地、建筑、淡水、海水區域大氣校正前后試驗影像樣本像元每種300個，分析大氣校正處理對不同地物類型影像光譜影響，比較校正前表觀反射率和校正后地表反射率差異（圖7）。

對于裸土、農田、林地、建筑4種地表下墊面類型，在可見光波段均表現出隨波長增加大氣校正處理使反射率差值差異減小，尤其在紅波段大氣校正處理與否反射率值變化不大的特征，這符合可見光短波波段易受大氣影響的規律。無論是否進行大氣校正處理，裸土及建筑可見光波段反射率均高于農田及林地、近紅外波段反射率均低于農田及林地。大氣校正前后，反射率隨波長變化曲線趨勢較為相似，無植被覆蓋的裸土及建筑大氣校正前各可見光波段反射率值相近，大氣校正后反射率隨波長逐漸增加；農田及林地大氣校正前可見光波段反射率隨波長增加減小，大氣校正使波長短的藍波段反射率明顯降低；4個陸表地類近紅外波段大氣校正后反射率值均略增加。水體反射率整體表現出與陸表相反的特征，兩類水體大氣校正前可見光波段反射率特征相似，均表現為隨波長降低，且海水反射率較淡水偏高，大氣校正處理使水體各波段反射率值不同程度降低。

在LaSRC大氣校正流程算法中，LC8 OLI氣溶膠的模擬應用到了近紅外和短波紅外波段歸一化指數。GF-1 WFV/MSS沒有短波紅外通道，在獲取GF-1 WFV/MSS氣溶膠時，利用歸一化植被指數進行了替代，雖然NDVI與NDVImir具有較好的線性關系，但在稀疏植被覆蓋地區差異大于植被覆蓋地區。利用前述采集的300個裸土和建筑大氣校正結果樣本數據代表稀疏植被覆蓋地區，以農田和林地代表植被覆蓋地區，以同時次對應樣本點LC8 OLI大氣校正結果為基準，比較稀疏/有植被覆蓋地區大氣校正效果。

從表3可以看出，對于裸土和建筑樣本數據GF-1 MSS大氣校正結果與LC8 OLI大氣校正結果相關系數分別為0.765～0.860和0.692～0.833，農田和林地樣本數據大氣校正結果相關系數分別為0.810～0.928和0.850～0.906。該方法在有農田和林木等植被覆蓋的中低緯度大氣校正效果較好，對稀疏植被的裸地和建筑地表大氣校正效果相對稍差。

表3 稀疏/有植被覆蓋地區大氣校正結果比較

3.3 大氣校正精度分析

Vermote等[33]已驗證LC8 OLI大氣校正結果與MODIS產品相比具有較高精度，而MODIS地表反射率產品普遍認為較為準確。本文以LC8 OLI數據大氣校正結果為基準，比較同日同區域GF-1 WFV/MSS各傳感器大氣校正結果準確性。所用數據為2016 年8月26日WFV2/MSS1/MSS2、2017年3月31日WFV1、2017年6月5日WFV4、2017年9月7日WFV3，這些數據涉及所有GF-1傳感器，且GF-1與LC8過境時傳感器采集數據有較大交叉覆蓋范圍便于二者進行比較分析。

將同日GF-1與LC8數據進行地理配準，對兩景影像分別出現的云、云陰影區域進行合并剔除，剪裁兩景影像交叉部位。根據通道波段范圍及中心波長，以GF-1 WFV/MSS通道1～4順序對應LC8 OLI通道2～5，即將2種衛星傳感器藍、綠、紅、近紅外波段兩兩對應，提取GF-1 和LC8兩景影像交叉區域所有對應像元，比較GF-1和LC8各景地表反射率算數平均值（）和二者差值（），逐像元平均絕對誤差（）、相關系數（）進行大氣校正結果的比較。

式中1、2分別為GF-1各傳感器及LC8 OLI每通道整景影像表觀反射率（TOA）或地表反射率（SR）平均值；1i、2i分別為GF-1各傳感器及LC8 OLI各景每像元大氣校正結果；為樣本容量。

所選取的每種傳感器的數據采集時間及下墊面不同，表4中GF-1及LC8各傳感器大氣校正前各通道表觀反射率平均值（GF1,TOA與LC8,TOA）相差較多。但無論各傳感器各通道平均大氣校正結果還是逐像元大氣校正結果比較，GF-1各傳感器與LC8 OLI地表反射率差異均較小。大氣校正后GF1各傳感器每通道地表反射率平均值（GF1,SR）與LC8 OLI對應通道地表反射率平均值（LC8,SR）藍通道相差0.13%～1.38%、綠通道相差為0.23%～1.84%、紅通道相差為0.71%～1.51%、近紅外通道相差為0.97%～6.64%；GF-1和LC8大氣校正結果逐像元絕對誤差，藍通道為0.81%～1.72%、綠通道為0.81%～1.94%、紅通道為0.95%～2.26%、近紅外通道為1.83%～4.39%；其中除WFV4近紅外通道校正結果相差稍多，其他誤差均低于3%。

表4 GF-1與LC8大氣校正結果比較

由于比較時進行了質檢剔除了云覆蓋數據，無損采集GF-1與LC8同日同地交叉地域覆蓋的其他所有有效像元，樣本容量都超過了千萬級別，GF-1各傳感器與LC8 OLI逐像元大氣校正結果都具有較高相關性，且均通過0.05顯著性檢驗。WFV1與OLI大氣校正結果相關系數為0.825～0.963、WFV2與OLI為0.909～0.944、WFV3與OLI為0.905～0.956、WFV4與OLI為0.951～0.972、MSS1與OLI為0.860～0.900、MSS2與OLI為0.875～0.920，總體上相關系數與波長呈正相關，這是由于波長越長對大氣穿透性越好，所需進行的校正相對較少，相應的引入誤差的可能性越小。相較于WFV系列傳感器，MSS 2個傳感器大氣校正結果與OLI相比相關性較低，是由于MSS空間分辨率8 m在像元重定位空間匹配時存在一定的偏差。從比較結果分析，對于GF-1 WFV/MSS傳感器應用LC8 OLI大氣校正算法進行大氣校正處理是適用的。

4 結 論

1）本文驗證了本地化LaSRC流程進行LC8 OLI數據大氣校正結果與USGS LC8 OLI地表反射率產品通道1～7決定系數分別為0.997、0.998、0.998、0.997、0.995、0.996、0.993，即相關性較好的基礎上，通過對逐像元幾何參數的計算、逐像元氣溶膠光學厚度的求算、MODIS產品的調用等系統建設，構建了GF-1 WFV/MSS數據大氣校正方法流程及應用平臺。

2）通過對試驗數據進行大氣校正處理，GF-1數據大氣校正后目視判讀性得到改善。研究表明，大氣校正使藍、綠通道地表反射率有所降低，大氣校正后的地表反射率能更好反映下墊面情況；陸表地類大氣校正結果隨下墊面不同表現出差異，GF-1 MSS大氣校正結果與LC8 OLI大氣校正結果相關系數在裸土和建筑樣本區域分別為0.765～0.860和0.692～0.833，農田和林地樣本區域分別為0.810～0.928和0.850～0.906。

3）應用同日同地交叉地域GF-1 WFV/MSS與LC8 OLI大氣校正結果分析表明，各通道地表反射率差值和逐像元絕對誤差較小；二者相關系數WFV1與OLI為0.825～0.963、WFV2與OLI為0.909～0.944、WFV3與OLI為0.905～0.956、WFV4與OLI為0.951～0.972、MSS1與OLI為0.860～0.900、MSS2與OLI為0.875～0.920，相關系數最低為0.825，相關關系較好。因此，對于GF-1 WFV/MSS數據應用LaSRC流程進行大氣校正處理能夠獲得較為可觀的精度。

4）與LC8 OLI相比GF-1缺少短波紅外通道，為繼承LC8 OLI大氣校正過程中氣溶膠光學厚度通過MODIS產品獲取待校正數據像元水平的做法，GF-1數據大氣校正處理時，求解氣溶膠時用NDVI線性斜率和截距代表過程變量NDVImir，實現了用待校正GF-1數據獲取其氣溶膠光學厚度的方法。近紅外波段大氣校正前后地表反射率表現恰好相反，各傳感器近紅外波段地表反射率較表觀反射率略有增加，且部分數值增加較多，并且在GF-1 WFV/MSS與LC8 OLI大氣校正結果比較時，近紅外通道校正結果也表現出偏差較大，由于LC8 OLI近紅外波段光譜響應范圍較GF-1相比偏窄，雖然利用光譜匹配系數對一些參數做了調整，難免還存在偏差，后續還需進行算法的改進及參數調整。

[1]邱鵬勛，汪小欽，茶明星，等. 基于TWDTW的時間序列GF-1 WFV農作物分類[J]. 中國農業科學，2019，52(17)：2951－2961.

Qiu Pengxun, Wang Xiaoqian, Cha Mingxing, et al. Crop identification based on TWDTW method and time series GF-1 WFV[J]. Scientia Agricultura Sinica, 2019, 52(17): 2951－2961. (in Chinese with English abstract)

[2]常布輝，王軍濤，羅玉麗，等. 河套灌區沈烏灌域GF-1/WFV遙感耕地提取[J]. 農業工程學報，2017，33(23)：188－195.

Chang Buhui, Wang Juntao, Luo Yüli, et al. Cultivated land extraction based on GF-1/WFV remote sensing in Shenwu irrigation area of Hetao Irrigation District[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(23): 188－195. (in Chinese with English abstract)

[3]Bernardo N, Watanabe F, Rodrigues F, et al. An investigation into the effectiveness of relative and absolute atmospheric correction for retrieval the TSM concentration in inland waters[J]. Modeling Earth Systems and Environment, 2016, 2(3): 1－7.

[4]趙英時. 遙感應用分析原理與方法[M]. 北京：科學出版社，2003.

[5]Richter R. A spatially adaptive fast atmospheric correction algorithm[J]. International Journal of Remote Sensing, 1996, 17(6): 1201－1214.

[6]Chavez P S. An improved dark-object subtraction technique for atmospheric scattering correction of multispectral data[J]. Remote Sensing of Environment, 1988, 24(3): 459－479.

[7]Kaufman Y J, Wald A E, Remer L A, et al. The MODIS 2.1-m channel-correlation with visible reflectance for use in remote sensing of aerosol[J]. IEEE Transactions on Geoscience and Remote Sensing, 1997, 35(5): 1286－1298.

[8]Hall F G, Strebel D E, Nickeson J E, et al. Radiometric rectification: Toward a common radiometric response among multi-date, multisensor images[J]. Remote Sensing of Environment, 1991, 35(1): 11－27.

[9]Kruse F A. Use of airborne imaging spectrometer data to map minerals associated with hydrothermally altered rocks in the northern grapevine mountains, Nevada and California[J]. Remote Sensing of Environment, 1988, 24(1): 31－51.

[10]Nieuwenhove V V, Beenhouwer J D, Carlo F D, et al. Dynamic intensity normalization using eigen flat field in X-ray imaging[J]. Optics Express, 2005, 23(21): 27975－27989.

[11]Smith G M, Milton E J. The use of the empirical line method to calibrate remotely sensed data to reflectance[J]. International Journal of Remote Sensing, 1999, 20(13): 2653－2662.

[12]Tanré D, Deroo C, Duhaut P, et al. Description of a computer code to simulate the satellite signal in the solar spectrum: the 5S code[J]. International Journal of Remote Sensing, 1990, 11(4): 659－668.

[13]Tanré D, Herman M, Deschamps P Y, et al. Atmospheric modeling for space measurements of ground reflectances, including bidirectional properties[J]. Applied Optics, 1979, 18(21): 3587－3594.

[14]Kotchenova S Y, Vermote E F, Matarrese R, et al. Validation of a vector version of the 6S radiative transfer code for atmospheric correction of satellite data. Part I: Path radiance[J]. Applied Optics, 2006, 45(26): 6762－6774.

[15]Kotchenova S Y, Vermote E F, Levy R, et al. Radiative transfer codes for atmospheric correction and aerosol retrieval: intercomparision study[J]. Applied Optics, 2008, 47(13): 2215－2226.

[16]劉佳，王利民，楊玲波，等. 基于6S模型的GF-1衛星影像大氣校正及效果[J]. 農業工程學報，2015，31(19)：159－168. Liu Jia, Wang Limin, Yang Lingbo, et al. GF-1 satelliate image atmospheric correction based on 6S model and its effcet[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(19): 159－168. (in Chinese with English abstract)

[17]孫林，于會泳，傅俏燕，等. 地表反射率產品支持的GF-1PMS氣溶膠光學厚度反演及大氣校正[J]. 遙感學報，2016，20(6)：216－228. Sun Lin, Yu Huiyong, Fu Qiaoyan, et al. Aerosol optical depth retrieval and atmospheric correction application for GF-1 PMS supported by land surface reflectance data[J]. Journal of Remote Sensing, 2016, 20(6): 216－228. (in Chinese with English abstract)

[18]孫元亨，秦其明，任華忠，等. GF-4/PMS與GF-1/WFV兩種傳感器地表反射率及NDVI一致性分析[J]. 農業工程學報，2017，33(9) ：167－173. Sun Yuanheng, Qin Qiming, Ren Huazhong, et al. Consistency analysis of surface reflectance and NDVI between GF-4/PMS and GF-1/WFV[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(9): 167－173. (in Chinese with English abstract)

[19]薛現光，周揚，胡校飛，等. 基于大氣參數擬合的GF-1多光譜影像大氣校正[J]. 海洋測繪，2018，38(5)：50－54.

Xue Xianguang, Zhou Yang, Hu Xiaofei, et al. Atmospheric correction of GF-1 multi-spectral image based on fitting of atmospheric parameters[J]. Hydrographic Surveying and Charting, 2018, 38(5): 50－54. (in Chinese with English abstract)

[20]胡勇，仲波，馬澤忠，等. 采用MODIS大氣產品的高分一號WFV數據大氣校正[J]. 遙感信息，2018，33(5)：35－40.

Hu Yong, Zhong Bo, Ma Zezhong, et al. GF-1 WFV atmospheric correction based on MODIS atmosphere produce and 6S model[J]. Remote Sensing Informaton, 2018, 33(5): 35－40. (in Chinese with English abstract)

[21]吳北嬰. 大氣輻射傳輸實用算法[M]. 北京：北京氣象出版社，1998.

[22]Berk A, Bernstein L S, Anderson G P, et al. MODTRAN cloud and multiple scattering upgrades with application to AVIRIS[J]. Remote Sensing of Environment, 1998, 65(3): 367－375.

[23]Callieco F, Dell-Acqua F. A comparison between two radiative transfer models for atmospheric correction over a wide range of wavelengths[J]. International Journal of Remote Sensing, 2011, 32(5): 1357－1370.

[24]楊貴軍，黃文江，劉三超，等. 環境減災衛星高光譜數據大氣校正模型及驗證[J] . 北京大學學報：自然科學版，2010，46(5)：821－828.Yang Guijun, Huang Wenjiang, Liu Sanchao, et al. Research on modeling and validating of atmospheric correction for HJ-1A hyperspectral imager data[J]. Acta Scientiarum Naturalium Universitatis Pekinensis, 2010, 46(5): 821－828. (in Chinese with English abstract)

[25]鄭盛，趙祥，張顥，等. HJ-1衛星CCD數據的大氣校正及其效果分析[J] . 遙感學報，2011，15(4)：709－721. Zheng Sheng, Zhao Xiang, Zhang Hao, et al. Atmospheric correction on CCD data of HJ-1 satellite and analysis of its effect[J]. Journal of Remote Sensing, 2011, 15(4): 709－721. (in Chinese with English abstract)

[26]熊世為，沈安云，李衛國，等. MODTRAN模型在HJ/CCD影像大氣校中的應用[J] . 江蘇農業學報，2016，32(2)：319－324. Xiong Shiwei, Shen Anyun, Li Weiguo, et al. Application of MODTRAN model in atmospheric correction of HJ/CCD data[J]. Jiangsu Journal of Agricultural, 2016, 32(2): 319－324. (in Chinese with English abstract)

[27]Cooley T, Anderson G P, Felde G W, et al. FLAASH, a MODTRAN4-based atmospheric correction algorithm ,its application and validation[J]. IEEE International Geosciences and Remote Sensing Symposium, 2002(3): 1414－1418.

[28]陳玲，陳理，李偉，等. 基于FLAASH模型的Worldview3大氣校正[J]. 國土資源遙感，2019，31(4)：26－31. Chen Ling, Chen Li, Li Wei, et al. Atmospheric correction of Worldview3 image based on FLAASH model[J]. Remote Sensing for Land and Resource, 2019, 31(4): 26－31. (in Chinese with English abstract)

[29]李衛國，蔣楠，王紀華. 基于薄云霧去除的ETM+影像大氣校正[J] . 農業工程學報，2013，29(增刊1)：82－88. Li Weiguo, Jiang Nan, Wang Jihua. Atmospheric correction for ETM+ image based on thin cloud removal[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2013, 29(Supp.1): 82－88. (in Chinese with English abstract)

[30]Richter R. A fast atmospheric correction algorithm applied to Landsat TM images[J]. International Journal of Remote Sensing, 1990, 11(1): 159－166.

[31]王正海，段建軍，耿欣. 基于波普匹配的Hyperion數據大氣校正方法對比研究[J] . 遙感技術與應用，2011，26(4)：432－436. Wang Zhenghai, Duan Jianjun, Geng Xin, et al. Comparative study atmospheric correction methods of Hyperion data based on spectral matching[J]. Remote Sensing for Land and Resource, 2011, 26(4): 432－436. (in Chinese with English abstract)

[32]Vermote E F, Justice C, Claverie M, et al. Preliminary analysis of the performance of the Landsat 8/OLI land surface reflectance product[J]. Remote Sensing of Environment, 2016, 185: 46－56.

[33]Vermote E F, Tanré D, Deuzé J L, et al. Second simulation of the satellite signal in the solar spectrum, 6S: An overview[J]. IEEE Geoscience and Remote Sensing Society, 1997, 35(3): 675－686.

[34]Vermote E F. Landsat8/OLI and sentinel 2 and VIIRS surface reflectance is largely based on MODIS C6(LaSRC) [EB/OL]. Workshop on Uncertainties in Remote Sensing ESA/ESRIN Frascati Italy [2017-10-25]. https://earth.esa.int/documents/ 700255/3264428/3. Vermote_ESRIN-WS-Uncertainties.pdf.

[35]Steven M D, Malthus T J, Baret F, et al. Intercalibration of vegetation indices from different sensor systems[J]. Remote Sensing of Environment, 2003, 88(4): 412－422.

[36]Tan B, Masek J G, Wolfe R, et al. Improved forest change detection with terrain illumination corrected Landsat images[J]. Remote Sensing of Environment, 2013, 136: 469－483.

[37]?ngstr?m A. The parameters of atmospheric turbidity[J]. Tellus, 1964, 16(1): 64－75.

[38]Dutton E G, Reddy P, Ryan S, et al. Features and effects of aerosol optical depth observed at Mauna Loa, Hawaii: 1982－1992[J]. Journal of Geophysical Research Atmospheres, 1994, 99(4): 8295－8306.

[39]Bodhaine B A, Norman B W, Ellsworth G D, et al. On Rayleigh optical depth calculations[J]. Journal of Atmospheric and Oceanic Technology, 1999, 16(11): 1854－1861.

Atmospheric correction method of GF-1 data based on Landsat8 product algorithm flow

Zhang Xiaoyue1,2,3, Li Linlin2,3, Wang Ying2,3, Zhang Qi2,3, Li Guochun4

(1.110166,; 2.110166,; 3.110166,; 4.,,110866,)

GF-1 satellite has the characteristics of high spatial resolution and short revisiting period, serving as the first satellite of the High-resolution Earth Observation System for National Science and Technology Major Project in China. The satellite carries two multi-spectral high-resolution cameras (panchromatic multispectral sensor, PMS) and four multi-spectral medium-resolution camera (wide field of view, WFV). GF-1 captured data play an important role in the identification of the underlying surface, and these data can be obtained free of charge from the website of the China Centre for Resources Satellite Data and Application. An important step for the application of GF-1 satellite data can be the interference removal of atmospheric molecules, aerosols, ozone, water vapor. The spectral response curves from Landsat8 (LC8) operational land imager(OLI) and GF-1 MSS/WFV were then analyzed in the visible and near the infrared bands. The results showed that the spectral range of LC8 OLI in the red- and near the infrared bands was relatively narrower than that of GF-1 MSS/WFV, whereas the spectral response function in the blue- and green bands was slightly different from that of GF-1 MSS/WFV, indicating that it is feasible to transplant LaSRC correction process to GF-1 MSS/WFV data. Since GF-1 satellite lacks the short-wave infrared band compared with LC8, the algorithm was modified to adapt to the characteristics of GF-1 channels. The atmospheric correction project was designed for the GF-1 satellite MSS/WFV data, including the algorithm analysis and code implement based on 6S atmospheric radiation simulation model and C++ programming language. Some parameters were used to estimate initial aerosol, including total atmospheric transmission, gaseous transmission, atmosphere spherical albedo and actual values of digital elevation model, atmospheric precipitation, ozone content. The loop calculation of the aerosol optical thickness(AOT) was carried out until the ratio between the red- and blue bands of GF-1 MSS/WFV data equal to the prescribed ratio of MODIS, according to the relationship between the blue- and the red surface reflectance known from MODIS. The results can be obtained the surface reflectance with the minimum residual error during different ?ngstr?m coefficients, and retrieved the aerosols in the pixel level of GF-1 MSS/WFV data. The pixel aerosol and these parameters were then substituted into 6S model to calculate the surface reflectance. Since the project was equipped a data file containing the atmospheric precipitation and ozone content at the current day, the surface reflectance could be obtained when only inputting GF-1 MSS/PMS data. Because the data of atmospheric influence gases, such as ozone and water vapor, on the same day of the data to be corrected were sometimes difficult to identify, two schemes can be provided, one is to use the data at that time, the other is to use the daily values of ozone and water vapor in past six years instead. The experimental results show that the proposed method has a good effect on the atmospheric correction in the middle and low latitudes that covered by vegetation, such as farmland and trees, but not good effect on that of the bare land and building surface that covered by sparse vegetation. Based on 6S model and LaSRC correction process, the correlation coefficient of the atmospheric correction between GF-1 MSS/WFV and LC8 OLI was from 0.825 to 0.972, indicating a high correlation of atmospheric correction results for two satellites. WFV similar spatial resolution to that of LC8 OLI was in good agreement with that of LC8 OLI atmospheric correction compared with that of MSS. The results show that it is convenient and operable for the GF-1 satellite data atmospheric correction method using the self-estimation aerosol parameters in the pixel level based on 6S model and LaSRC process. This promising atmospheric method can be very suitable for the land surface application, such as agricultural and forestry monitoring in growing season. At present, this method has been successfully implemented on Remote Sensing data processing platform Remote Sensing Desktop (RSD) in China.

satellites; remote sensing; GF-1 WFV/MSS; Landsat8 OLI; 6S; LaSRC; atmospheric correction

張曉月，李琳琳，王 瑩，張 琪，李國春. 采用Landsat8產品算法流程的高分一號數據大氣校正[J]. 農業工程學報，2020，36(1)：182－192.doi：10.11975/j.issn.1002-6819.2020.01.021 http://www.tcsae.org

Zhang Xiaoyue, Li Linlin, Wang Ying, Zhang Qi, Li Guochun. Atmospheric correction method of GF-1 data based on Landsat8 product algorithm flow[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(1): 182－192. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2020.01.021 http://www.tcsae.org

2019-07-02

2019-12-02

中國氣象局沈陽大氣環境研究所開放基金（2017SYIAE03）

張曉月，高級工程師，主要從事農業氣象和遙感數據處理及應用研究。Email：zhangxiaoyue225@163.com

10.11975/j.issn.1002-6819.2020.01.021

TP701

A

1002-6819(2020)-01-0182-11