商子健，金建文，許振峰,李 鶴

(1.沈陽(yáng)市勘察測(cè)繪研究院有限公司,遼寧 沈陽(yáng) 110004；2.遼寧工程技術(shù)大學(xué) 測(cè)繪與地理科學(xué)學(xué)院，遼寧 阜新 123000)

0 引言

高分一號(hào)(GF-1)衛(wèi)星是中國(guó)高分辨率對(duì)地觀測(cè)系統(tǒng)的第一顆衛(wèi)星，于2013年4月26日在酒泉衛(wèi)星發(fā)射中心發(fā)射，具有空間分辨率高、重訪周期短、覆蓋寬等特點(diǎn)。GF-1衛(wèi)星搭載有2臺(tái)多光譜相機(jī)，多光譜空間分辨率為8 m，有藍(lán)、綠、紅和近紅外4個(gè)波段，具有周期為4天的重訪能力[1]。為我國(guó)國(guó)土、農(nóng)業(yè)、環(huán)保等部門提供高效、連續(xù)、穩(wěn)定的空間觀測(cè)服務(wù)，在自然災(zāi)害評(píng)估、生態(tài)環(huán)境監(jiān)測(cè)、城市和交通管理等領(lǐng)域發(fā)揮重要作用。

遙感成像時(shí)，大氣對(duì)遙感器入瞳信號(hào)的貢獻(xiàn)可超過(guò)80%[2],使其獲取的遙感信號(hào)不能反映地物的真實(shí)信息。因此，為了獲取地物的真實(shí)反射率數(shù)據(jù)，需要對(duì)遙感影像進(jìn)行大氣校正，消除大氣分子、氣溶膠等大氣組分造成的吸收和瑞利散射影像[3]。目前，大氣校正的方法主要有基于圖像本身的暗像元法、基于地面的線性經(jīng)驗(yàn)?zāi)Ｐ头ā⒒诖髿馔叫Ｕ齼x的大氣校正和基于輻射傳輸模型法等4種方法[4-6]?；趫D像本身的暗像元法要求影像具有反演大氣參數(shù)的短波紅外波段[7-8]?；诘孛娴木€性經(jīng)驗(yàn)?zāi)Ｐ头ㄐ枰赖匚镎鎸?shí)的光譜信息[9]，這2種方法在缺少輔助數(shù)據(jù)情況下，難以開(kāi)展GF-1影像的大氣校正?；诖髿馔叫Ｕ齼x的大氣校正要求衛(wèi)星上搭載同步的大氣校正，而我國(guó)GF-1衛(wèi)星缺少大氣校正儀，因此基于大氣同步校正儀的大氣校正方法也不適應(yīng)GF-1影像大氣校正?；谳椛鋫鬏斈Ｐ偷拇髿庑Ｕ▽儆谖锢砟Ｐ头椒ǎ褂媚Ｐ吞峁┑哪M大氣參數(shù)可以進(jìn)行GF-1影像的大氣校正。6S輻射傳輸模型配置了典型大氣模式、氣溶膠類型模式和常用的傳感器的光譜響應(yīng)函數(shù)，實(shí)用性較強(qiáng)，廣泛應(yīng)用于光學(xué)影像的大氣校正。因此，本文使用6S輻射傳輸模型進(jìn)行GF-1影像進(jìn)行大氣校正研究。

針對(duì)輔助數(shù)據(jù)不足和GF-1衛(wèi)星缺少反演大氣參數(shù)的短波紅外波段，以及無(wú)法使用基于大氣同步校正儀的GF-1大氣校正等問(wèn)題，本文基于6S輻射傳輸模型，設(shè)計(jì)并實(shí)現(xiàn)了適合GF-1衛(wèi)星數(shù)據(jù)大氣校正算法與程序，并用天津城區(qū)和山區(qū)兩景影像做算法驗(yàn)證。由于沒(méi)有實(shí)測(cè)光譜數(shù)據(jù)，本文選用商業(yè)軟件ENVI提供的FLAASH大氣校正結(jié)果驗(yàn)證了算法的大氣校正效果。

1 研究方法

1.1 算法介紹

本文基于6S模型獲取大氣頂輻亮度轉(zhuǎn)化為地表反射率的相關(guān)參數(shù)，根據(jù)轉(zhuǎn)化參數(shù)計(jì)算每景影像的地表反射率。算法以GF-1星8 m分辨率多光譜數(shù)據(jù)，其波段設(shè)置如表1所示。

表1 GF-1多光譜影像波段設(shè)置Tab.1 Band setting of GF-1 multispectral image 單位：nm

元數(shù)據(jù)及傳感器公開(kāi)參數(shù)為輸入?yún)?shù)，無(wú)需其他輔助數(shù)據(jù)。

基于6S模型GF-1影像大氣校正示意如圖1所示，主要包括輻射定標(biāo)、6S參數(shù)設(shè)置和6S大氣校正3部分。其中，DEM是全球0.01°空間分辨率高程數(shù)據(jù)；GF-1輻射定標(biāo)數(shù)據(jù)從中國(guó)資源衛(wèi)星應(yīng)用中心官網(wǎng)下載；2.5 nm光譜重采樣函數(shù)由中國(guó)資源衛(wèi)星應(yīng)用中心提供的GF-1光譜響應(yīng)函數(shù)重采樣得到。

圖1 6S模型GF-1影像大氣校正示意Fig.1 Atmospheric correction process of 6S model GF-1 image

1.2 輻射定標(biāo)

在用6S模型對(duì)影像大氣校正前，采用式(1)對(duì)影像進(jìn)行絕對(duì)輻射定標(biāo)：

Lλ=Gain×DN，

(1)

式中，Lλ為波段的表觀輻射亮度(W·m-2·sr-1·μm-1)；Gain為增益；DN為原始影像像元的灰度值。式(1)中的增益系數(shù)在影像的頭文件中可以獲得。

1.3 表觀反射率計(jì)算

表觀反射率也稱大氣頂層反射率，表示衛(wèi)星傳感器接收的光譜輻亮度與大氣頂層太陽(yáng)輻亮度的比值[10]。利用式(2)可將表觀輻亮度轉(zhuǎn)換為表觀反射率：

(2)

式中，d為日地距離校正因子；ESUNλ為中心波長(zhǎng)為λ的大氣上界太陽(yáng)光譜輻照度(W·m-2·μm-1)，其值可根據(jù)GF-1多光譜影像的光譜響應(yīng)函數(shù)與大氣上界太陽(yáng)光譜輻照度，通過(guò)式(3)積分得到：

(3)

式中，Eλ為大氣上界太陽(yáng)光譜輻照度；fλ為光譜響應(yīng)函數(shù)；λ1,λ2為光譜響應(yīng)函數(shù)的起始和終止波長(zhǎng)。

GF-1多光譜影像各波段的大氣上界太陽(yáng)光譜輻照度值如表2所示，θ為太陽(yáng)天頂角,由GF-1頭文件得到。

表2 GF-1多光譜影像大氣上界太陽(yáng)光譜輻照度Tab.2 Upper atmospheric solar spectral irradiance of GF-1 multispectral image 單位：W·m-2·μm-1

1.4 6S參數(shù)設(shè)置及大氣校正

本文對(duì)6S原始模型添加了GF-1衛(wèi)星支持，并編譯為可執(zhí)行文件，用Python語(yǔ)言編寫代碼，自動(dòng)調(diào)用6S程序逐波段計(jì)算大氣校正系數(shù)，得到校正后的地表反射率。整個(gè)校正過(guò)程中，程序自動(dòng)讀取GF-1原始影像及元數(shù)據(jù)，用戶只需確定影像的能見(jiàn)度或氣溶膠厚度，即可進(jìn)行地表反射率自動(dòng)化計(jì)算。主要步驟如下：

① 6S模型參數(shù)設(shè)置：利用幾何參數(shù)(太陽(yáng)天頂角和方位角、衛(wèi)星天頂角和方位角、傳感器高度、成像時(shí)間、全球0.01°空間分辨率的高程數(shù)據(jù)、2.5 nm光譜重采樣響應(yīng)函數(shù))和大氣參數(shù)(大氣模式、氣溶膠類型)配置6S輸入文件；

② 調(diào)用6S程序，得到GF-1影像各波段的大氣校正系數(shù)xa、xb、xc；

③ 利用6S模型提供的式(4)和式(5)計(jì)算得到校正后地表反射率[11]：

y=xa·Lλ-xb，

(4)

(5)

2 實(shí)驗(yàn)結(jié)果與分析

調(diào)用6S程序后，所得到GF-1影像各波段大氣校正系數(shù)的統(tǒng)計(jì)結(jié)果如表3所示。

表3 GF-1衛(wèi)星多光譜影像6S大氣校正輸出參數(shù)Tab.3 Output parameters of 6S atmospheric correction for GF-1 multispectral imagery

2.1 地物光譜對(duì)比分析

本文選取天津2017年2月27日城區(qū)和2016年4月9日山區(qū)的GF-1星8 m分辨率的多光譜數(shù)據(jù)進(jìn)行算法實(shí)驗(yàn)，并用ENVI軟件中FLAASH大氣校正結(jié)果作為對(duì)比。圖2，圖3分別是天津城區(qū)和山區(qū)大氣校正前后真彩色合成效果圖。對(duì)比發(fā)現(xiàn)，經(jīng)過(guò)本文算法大氣校正后，圖像的清晰度、對(duì)比度明顯提升。

(a)城區(qū)

(a)城區(qū)

大氣校正后，目視選擇城區(qū)水體、水泥地和山區(qū)植被、裸地等每種典型地物的20個(gè)純像元，計(jì)算其平均反射率，以定量比較大氣校正前后的結(jié)果[12-13]。典型地物大氣校正前后反射率對(duì)比如表4所示。

表4 典型地物大氣校正前后反射率對(duì)比Tab.4 Reflectance comparison of typical surface features before and after atmospheric correction

由表4可以看出，本文算法在藍(lán)、綠、紅波段大氣校正后，除了裸地在紅波段校正后反射率反而增大，其余地物校正后反射率都小于校正前；在近紅外波段，除了水體校正后的反射率減小，其余地物在校正后的反射率都有所增加。所以，本文算法較好的去除了大氣的影響，使各地物的反射率趨于正常值，消除了在藍(lán)、綠、紅波段氣溶膠散射對(duì)地表反射的增強(qiáng)，有效校正了在近紅外波段水汽對(duì)地表反射率的削弱。

水體大氣校正、水泥地大氣校正、植被大氣校正和裸地大氣校正前后反射率變化分別如圖4～圖7所示。

(a)校正前

(a)校正前

(a)校正前

(a)校正前

從圖4～圖7可以看出，經(jīng)過(guò)本文算法大氣校正后，4種典型地物的光譜曲線均符合其光譜特性，并與FLAASH大氣校正模型具有良好的一致性，表明本文算法在城區(qū)和山區(qū)影像大氣校正中具有一定的普適性。

2.2 大氣校正對(duì)NDVI的影響分析

植被指數(shù)是植被定量遙感中常用的參數(shù)，通常是由2個(gè)或2個(gè)以上的波段計(jì)算得到。歸一化植被指數(shù)(NDVI)能夠有效檢驗(yàn)大氣校正的效果[14-15]。本文選用NDVI作為指標(biāo)，對(duì)大氣校正前后影像的NDVI值變化進(jìn)行分析，其計(jì)算式如下：

(6)

式中，ρnir，ρred分別為近紅外和紅外波段的反射率。

選取天津城區(qū)水體、水泥地和山區(qū)植被、裸地等典型地物提取植被指數(shù)信息[16]。經(jīng)式(6)計(jì)算得到典型地物大氣校正前后NDVI的變化趨勢(shì)，如圖8所示。

圖8 大氣校正前后典型地物NDVI的變化Fig.8 Changes of NDVI of typical ground objects before and after atmospheric correction

由圖8可以看出，本文算法很好地校正了大氣對(duì)地物NDVI的降低，且與FLAASH模型大氣校正后地物的NDVI值變化趨勢(shì)一致，說(shuō)明本文大氣校正算法在城區(qū)和山區(qū)都能較好的去除大氣影像，還原地物的NDVI值。其中，植被NDVI增幅最大，從校正前的0.079 1增加到校正后的0.252 1(見(jiàn)表4)。大氣校正前后NDVI統(tǒng)計(jì)如表5所示。

表5 大氣校正前后NDVI統(tǒng)計(jì)表Tab.5 NDVI statistics before and after atmospheric correction

由表5可以看出，本文大氣校正后的NDVI的均值和標(biāo)準(zhǔn)差都大于校正前，即校正后影像目視效果更佳，更有利于進(jìn)一步進(jìn)行植被遙感的定量化研究。

2.3 大氣校正效率分析

在內(nèi)存64 GB、CPU主頻2.10 GHz、磁盤剩余空間200 GB的運(yùn)行環(huán)境下，F(xiàn)LAASH在參數(shù)配置完畢后運(yùn)行總耗時(shí)252 s，而本文方法總耗時(shí)118 s，單景影像的處理效率提高了53%。因此，本文算法較FLAASH相比，提高了大氣校正效率。

3 結(jié)束語(yǔ)

本文基于6S輻射傳輸模型，設(shè)計(jì)并實(shí)現(xiàn)了適合GF-1 衛(wèi)星數(shù)據(jù)大氣校正算法與程序，并以FLAASH大氣校正結(jié)果作為參考，對(duì)校正結(jié)果進(jìn)行了驗(yàn)證，得出結(jié)論：① 通過(guò)典型地物校正前后的反射率和NDVI比較表明，本文算法較好的還原了水體、水泥地、植被和裸地等地物類型的光譜曲線，增加了地表的NDVI值，并與FLAASH校正結(jié)果趨勢(shì)具有較好的一致性；② 本文大氣校正算法基于6S輻射傳輸模型，校正過(guò)程自動(dòng)化，效率高，為GF-1影像的業(yè)務(wù)化大氣校正提供了參考；③ 由于沒(méi)有地物真實(shí)反射率數(shù)據(jù)，無(wú)法對(duì)校正的結(jié)果進(jìn)一步定量分析，下一步需要增加地面實(shí)測(cè)數(shù)據(jù)對(duì)校正結(jié)果進(jìn)一步分析；其次，本文假設(shè)地表均一，同時(shí)選取550 nm氣溶膠光學(xué)厚度，這可能導(dǎo)致更大的大氣校正誤差，下一步研究需要考慮非均一地表帶來(lái)的臨近像元效應(yīng)的影響以及氣溶膠光學(xué)厚度的準(zhǔn)確獲取，來(lái)進(jìn)一步提升大氣校正效果。