孔卓,楊海濤,鄭逢杰,李揚(yáng),齊濟(jì)

(1.航天工程大學(xué) 研究生院,北京 101407;2.航天工程大學(xué) 太空安全研究中心,北京 101407;3.航天工程大學(xué) 航天信息學(xué)院,北京 101407)

0 引言

高分五號(Gaofen-5,GF-5)是我國高分專項(xiàng)中的一顆重要衛(wèi)星,也是世界首顆能夠?qū)崿F(xiàn)對大氣和陸地綜合觀測的全譜段高光譜衛(wèi)星[1],首星于2018年5月9日發(fā)射,02星于2021年9月7日發(fā)射。GF-5共搭載有6臺不同的傳感器,其中的可見光短波紅外高光譜傳感器(advanced hyperspectral imager,AHSI)光譜范圍400～2 500 nm,空間分辨率30 m,可見光波段和短波紅外波段光譜分辨率分別為5 nm和10 nm,是承擔(dān)對地高光譜分辨率觀測的重要手段[2]。

作為一種高光譜衛(wèi)星,GF-5在成像過程中接收到的地表輻射會受到大氣的影響,致使圖像不能正確反映真實(shí)地物的性質(zhì),這嚴(yán)重影響了對GF-5圖像的判讀與應(yīng)用。因此,在實(shí)現(xiàn)對GF-5圖像的應(yīng)用前,應(yīng)首先對其進(jìn)行大氣校正,最大限度還原地物真實(shí)光譜信息。

實(shí)現(xiàn)快速高效、自動化批量化的大氣校正能夠有效促進(jìn)GF-5圖像的應(yīng)用,但總體看,目前對GF-5大氣校正的研究相對較少。現(xiàn)有研究中,劉梓欽等[3]使用光譜超立方體快速視線大氣分析模型(fast line-of-sight atmospheric analysis of spectral hypercubes,FLAASH)和6S模型對GF-5圖像進(jìn)行了大氣校正,并與實(shí)測光譜數(shù)據(jù)進(jìn)行了對比,發(fā)現(xiàn)兩種方法均能有效實(shí)現(xiàn)大氣校正,且6S效果略優(yōu)于FLAASH;Feng[4]采用中等光譜分辨率大氣透射模型(the moderate spectral resolution atmospheric transmittance algorithm and computer model,MODTRAN)完成了GF-5圖像的大氣校正,與FLAASH校正結(jié)果及實(shí)測數(shù)據(jù)相比,MODTRAN大氣校正精度有明顯提高。這些研究已實(shí)現(xiàn)并驗(yàn)證了采用輻射傳輸模型實(shí)現(xiàn)GF-5圖像大氣校正的有效性,但由于輻射傳輸模型需輸入各種大氣參數(shù)(如氣溶膠、水汽含量等),當(dāng)缺失相應(yīng)參數(shù)時(shí),將導(dǎo)致大氣校正無法完成或效果下降,現(xiàn)有研究并未對此進(jìn)行考慮,因此應(yīng)加強(qiáng)對大氣校正同步獲取大氣參數(shù)的研究。此外,現(xiàn)有研究均是直接應(yīng)用輻射傳輸模型實(shí)現(xiàn)GF-5圖像的大氣校正,每次大氣校正均需輸入相關(guān)參數(shù),這難以實(shí)現(xiàn)圖像自動化批量化的大氣校正處理,因此也具有一定的局限性。

本文結(jié)合現(xiàn)有GF-5及其他高光譜圖像大氣校正的研究成果,提出了一種基于MODTRAN查找表的GF-5圖像大氣校正方法,該方法基于MODTRAN輻射傳輸模型先后構(gòu)造能見度查找表、水汽含量查找表和反射率查找表,在大氣校正前實(shí)現(xiàn)了大氣參數(shù)的獲取,同時(shí)查找表也有效代替了輻射傳輸模型的運(yùn)行,大氣校正過程無需輸入其他參數(shù),從而實(shí)現(xiàn)GF-5圖像自動化、批量化的大氣校正,為GF-5圖像大氣校正的研究提供參考。

1 大氣校正原理和方法

1.1 實(shí)驗(yàn)數(shù)據(jù)及計(jì)算過程

實(shí)驗(yàn)數(shù)據(jù)為1景2019年10月15日拍攝于北京懷柔地區(qū)的GF-5 AHSI圖像(圖1),圖像中心經(jīng)度為116.58°E,中心緯度為40.68°N,圖像無云,地表覆蓋類型主要為植被、裸土、水體、城鎮(zhèn)等,地表類型較為復(fù)雜,地表平均海拔875 m。

圖1 本文所用GF-5圖像

GF-5 AHSI圖像共有330個光譜通道,其中可見光近紅外波段(visible-near infrared,VNIR)為150個光譜通道,短波紅外波段(shortwave infrared,SWIR)為180個光譜通道。同時(shí),VNIR后6個通道和SWIR前4個通道波長重疊,可將其進(jìn)行合并,本文保留SWIR通道,合并后通道數(shù)為324個。

完成波段合并后,本文方法的主要步驟如圖2所示。首先對GF-5圖像完成輻射定標(biāo),之后調(diào)用MODTRAN模型分別構(gòu)造能見度查找表、水汽查找表和反射率查找表,逐像元完成能見度和水汽含量的反演,最后將相關(guān)參數(shù)帶入反射率查找表,完成GF-5圖像的大氣校正。

圖2 基于MODTRAN查找表的大氣校正步驟圖

本實(shí)驗(yàn)采用Python編程,并使用Python調(diào)用MODTRAN完成模型運(yùn)行,計(jì)算過程中分別在反射率查找表的構(gòu)建、水汽含量、能見度的逐像元反演、圖像的逐像元大氣校正過程中采用并行運(yùn)算進(jìn)行加速。對比實(shí)驗(yàn)采用FLAASH輻射傳輸模型進(jìn)行實(shí)現(xiàn),FLAASH是一種在MODTRAN基礎(chǔ)上開發(fā)的輻射傳輸模型,它利用圖像數(shù)據(jù)與大氣特性,以物理方法精確推算地表反射率,并在計(jì)算過程中考慮了鄰近效應(yīng)的影響,同時(shí)能夠校正傾斜觀測導(dǎo)致的誤差,計(jì)算精度得到了很大提升,已廣泛應(yīng)用于高光譜和多光譜圖像的大氣校正中[5]。

1.2 輻射定標(biāo)

遙感圖像記錄的原始數(shù)據(jù)(digital number,DN值)不能直接應(yīng)用于后續(xù)圖像的處理中,首先需通過輻射定標(biāo)將其轉(zhuǎn)化為表觀輻亮度或表觀反射率[6],對于GF-5,輻射定標(biāo)如式(1)所示。

LTOA=DN×gain

(1)

式中:LTOA為表觀輻亮度;gain為增益值,可從圖像數(shù)據(jù)中讀取。

根據(jù)表觀輻亮度可求得相應(yīng)波段的表觀反射率,如式(2)所示。

(2)

1.3 MODTRAN大氣校正

MODTRAN是由美國光譜科技公司、空間物理實(shí)驗(yàn)室聯(lián)合開發(fā)的一種輻射傳輸模型,它采用FORTRAN語言編寫,可以計(jì)算波數(shù)為0～50 000 cm-1的大氣透射率和輻射,在高光譜圖像的大氣校正中有著廣泛的應(yīng)用[7]。

采用MODTRAN實(shí)現(xiàn)大氣校正所采用的輻射傳輸方程[8]如式(3)所示。

(3)

式中:L為表觀輻亮度;Lp為程輻射;Gb、Gt為校正參數(shù),分別表示大氣的直射和散射影響;S為大氣反照率,可以用來解釋大氣與地表間多重散射的影響;ρ為地表反射率。式中L已知,通過運(yùn)行MODTRAN模型可得到其余4個參數(shù),由此可計(jì)算得到地表反射率ρ。

由于MODTRAN無交互式頁面,對其的運(yùn)行需通過編寫相應(yīng)的執(zhí)行文件(tape5文件),并運(yùn)行“modtran.exe”實(shí)現(xiàn),通過運(yùn)行可以得到多個文件,其中的tape7文件中存儲有一系列與輻射相關(guān)的參數(shù),大氣校正中主要需用到GRND RFLT、DRCT RFLT、TOTAL RAD、PTH THRML、SOL SCAT這幾個參數(shù),為便于后續(xù)計(jì)算,可將其合并簡化表示,用TOTR表示TOTAL RAD,用GRFL表示GRND RFLT,用DRFL表示DRCT RFLT,用PATH表示PTH THRML+SOL SCAT。

于是可求得大氣校正所需參數(shù),如式(4)～式(7)所示。

Lp=TOTR0=PATH0

(4)

(5)

Gb=(PATH1.0-PATH0)·(1-S)

(6)

Gt=(DRFL1.0-DRFL0)·(1-S)

(7)

式中:下標(biāo)中的0、0.5、1.0表示運(yùn)行MODTRAN時(shí)的反射率的大小。

1.4 大氣參數(shù)反演

MODTRAN的運(yùn)行需輸入各種大氣參數(shù),其中氣溶膠和水汽由于時(shí)空異變性大,且對大氣校正的結(jié)果影響顯著,是最關(guān)鍵和最難確定的參數(shù)。因此在實(shí)現(xiàn)大氣校正前應(yīng)首先確定其數(shù)值,本文分別采用V5.2法和APDA法反演獲得研究區(qū)域的氣溶膠和水汽含量。

1)氣溶膠反演。大氣校正中常用550 nm處的氣溶膠光學(xué)厚度或能見度表示氣溶膠對輻射傳輸?shù)挠绊?在MODTRAN中,二者的關(guān)系如式(8)所示。

(8)

式中:VIS為能見度;EXT550為550 nm氣溶膠光學(xué)厚度。

目前常見的氣溶膠反演方法主要有偏振反演[9]、暗目標(biāo)法[10]、結(jié)構(gòu)函數(shù)法[11]、深藍(lán)算法[12]等,其中暗目標(biāo)法應(yīng)用最為廣泛。近年來,在暗目標(biāo)法的基礎(chǔ)上開發(fā)出了應(yīng)用于MODIS氣溶膠反演的Version 5.2(V5.2)算法[13],相對于暗目標(biāo)法,V5.2算法精度更高,效果更好,由于GF-5 AHSI圖像包含有V5.2算法所需的相應(yīng)波段,因此可采取該方法實(shí)現(xiàn)能見度的反演。

2)水汽反演。在近紅外波段,常見的基于像素的水汽反演方法有連續(xù)內(nèi)插波段比值法(continuous interpolation band ratio,CIBR)、大氣預(yù)處理微分吸收法(atmospheric preconditioning differential absorption,APDA)、二階導(dǎo)數(shù)算法(second derivative algorithm,SODA)、光譜擬合技術(shù)等。CIBR方法在高反射率(>0.4)的表面上效果最好,但APDA也可以在較暗的曲面(反射率>0.1)上獲得良好的性能。因此,本文選擇APDA方法進(jìn)行水汽反演。

APDA方法最少需要兩個波段,一個為位于大氣窗口區(qū)域(不受水汽影響)的參考波段,一個為位于大氣水汽吸收區(qū)域的測量(計(jì)算)波段。由于GF-5波段數(shù)充足,可選擇3個通道(1個水汽吸收通道、2個參考通道)實(shí)現(xiàn)水汽反演。根據(jù)王崇倡等[14]的建議,選擇GF-5 AHSI圖像通道合并后的第159通道(SWIR第15波段)作為水汽吸收通道,149(SWIR第5波段)和169(SWIR第25波段)通道作為參考通道,因此可求得RAPDA[15]。此外,RAPDA可通過擬合公式計(jì)算[16]。于是,設(shè)置好初始水汽含量,即可迭代求得最終水汽含量,如圖3所示。

圖3 APDA法求解水汽含量流程圖

1.5 查找表構(gòu)建

1)能見度查找表。根據(jù)V5.2法反演能見度原理,能見度查找表中所需參數(shù)包括:不同能見度、470 nm和660 nm波段地表反射率。

由于需要計(jì)算地表反射率,需對470 nm和660 nm波段設(shè)置0、0.5和1.0反射率及不同能見度得到。對于能見度的取值,根據(jù)Gauanter等[17]的建議,可取10～200 km,步長統(tǒng)一為5 km。

2)水汽查找表。根據(jù)APDA法反演水汽原理,水汽查找表中所需參數(shù)包括:不同水汽含量、水汽吸收波段及參考波段程輻射、RAPDA。

因此,可設(shè)置反射率為0、0.5,能見度為圖像平均能見度,在不同水汽含量下運(yùn)行MODTRAN模型得到程輻射,之后計(jì)算得到相應(yīng)RAPDA。

根據(jù)Gauanter等的建議,水汽可取0.3～5.0 cm,因此本文在計(jì)算過程中設(shè)置水汽取值步長0.3～3.0 cm為0.1 cm,3.0～5.0 cm為0.3 cm。

3)反射率查找表。根據(jù)胡順石等[18]總結(jié)的構(gòu)建MODTRAN查找表的原則,同時(shí)考慮GF-5中已有數(shù)據(jù),本文構(gòu)建的GF-5反射率查找表共設(shè)置水汽含量、能見度、觀測天頂角、太陽天頂角、太陽方位角、地面海拔6個變化參數(shù),相應(yīng)取值及步長如表1所示。對于查找表中沒有的數(shù)值,可采用反距離插值法得到。

表1 反射率查找表參數(shù)設(shè)置

與水汽和能見度查找表相比,反射率查找表所變化的參數(shù)更多,計(jì)算時(shí)間也更長。為提升查找表構(gòu)建速度,可在MODTRAN運(yùn)行過程中采用并行運(yùn)算進(jìn)行加速。本文在查找表構(gòu)建過程中采用16線程并行運(yùn)算,計(jì)算時(shí)間縮短為原計(jì)算時(shí)間的約1/16。

上述3個查找表只需一次建立即可滿足不同觀測條件和地理幾何條件下GF-5圖像的大氣校正,而對于一副圖像的校正,由于GF-5數(shù)據(jù)中包含有圖像的觀測天頂角、太陽天頂角、太陽方位角數(shù)據(jù),同時(shí)GF-5數(shù)據(jù)中包含有4個角點(diǎn)和中心點(diǎn)的經(jīng)緯度信息,經(jīng)插值可以得到每個像元的經(jīng)緯度值,通過DEM文件即可求得每個像元的海拔值,因此可固定觀測天頂角、太陽天頂角、太陽方位角取值,對能見度、水汽含量及海拔逐像元取值,從而實(shí)現(xiàn)該圖像的逐像元大氣校正。

1.6 評價(jià)指標(biāo)

采用主觀評價(jià)法和客觀評價(jià)法相結(jié)合的方法實(shí)現(xiàn)對本文方法大氣校正效果的質(zhì)量評價(jià)。

1)主觀評價(jià)法即由判讀人員目視實(shí)現(xiàn)對大氣校正質(zhì)量的比較。由于這種方法帶有較大的主觀性,本文僅將其作為一種輔助評價(jià)方法。

2)客觀評價(jià)法。本文選取了均方根誤差(root mean square error,RMSE)、決定系數(shù)(coefficient of determination,R2)、歸一化植被指數(shù)(normalized difference vegetation index,NDVI)這3種廣泛使用的指標(biāo)完成對本文大氣校正方法效果的評價(jià)。

均方根誤差是真實(shí)值和預(yù)測值之間差異的標(biāo)準(zhǔn)差,用于說明兩種數(shù)據(jù)間的離散程度,值越小,表示兩組數(shù)據(jù)間的差異越小。

決定系數(shù)用于說明真實(shí)值與預(yù)測值間的相似程度,取值范圍為0～1,值越接近于1,說明兩組數(shù)據(jù)間的擬合效果越好。

歸一化植被指數(shù)定義為近紅外通道與可見光通道反射率之差與之和的商,可用于反映土地覆蓋植被狀況,預(yù)測值越接近于真實(shí)值,說明大氣校正效果越好。

2 大氣校正結(jié)果與評價(jià)

2.1 能見度、水汽含量反演效果評價(jià)

1)目視效果對比。對于能見度,由于FLAASH對暗目標(biāo)的選取為分塊處理,并且不能生成能見度圖,因此此節(jié)僅給出采用本文方法得到的暗目標(biāo)圖和能見度分布灰度圖,如圖4所示,其中圖中黑色像元為暗目標(biāo)。

圖4 能見度反演

根據(jù)Kaufman等[19]、Teillet等[20]總結(jié)的暗目標(biāo)選取原則,采用暗目標(biāo)法時(shí)一般可取濃密植被、清潔水體或陰影區(qū)域作為暗目標(biāo)。結(jié)合圖1和圖4(a)可以看出,本文方法所選取的暗目標(biāo)與暗目標(biāo)的選取原則相一致,同時(shí)由于考慮了條帶噪聲,圖中出現(xiàn)了多條非暗目標(biāo)條帶。

采用本文方法及FLAASH得到的水汽分布如圖5所示,圖中像元灰度越低,水汽含量越低。

圖5 水汽含量分布灰度圖

由圖5可以看出,本文方法與FLAASH得到的水汽含量分布較為一致。在水汽含量的反演過程中,除圖像右下方的水體外,由于部分像元數(shù)據(jù)存在異常,使得反演得到的水汽含量為負(fù)值或非常大,FLAASH將其進(jìn)行了剔除,因此圖5(b)中存在一些黑色的像元,即被剔除水汽含量的像元,而本文則是對存在問題的像元進(jìn)行插值,從而使每個像元都具有水汽含量,因此圖5(a)中沒有黑色像元。

2)反演結(jié)果對比。采用本文方法及FLAASH反演得到的能見度和水汽含量如表2所示。

表2 本文方法及FLAASH反演能見度、水汽含量值

由表3可知,本文方法計(jì)算得到的水汽和能見度與FLAASH計(jì)算得到的能見度和水汽含量有一定的差別,但數(shù)值相差較小,說明了本算法在計(jì)算能見度和水汽含量時(shí)具有可行性。計(jì)算得到的最大值和最小值相差較大,分析原因主要是由于云的影響和部分像元存在壞數(shù)據(jù),以及不同地表情況的影響,導(dǎo)致這些像元的能見度和水汽含量與整體相差較大,由此也說明了本算法通過逐像元的能見度和水汽含量的反演,能夠更好地反映出不同地物受到大氣的影響,從而更加精確地實(shí)現(xiàn)大氣校正。

2.2 大氣校正效果評價(jià)

1)目視效果對比。經(jīng)本文方法和FLAASH大氣校正后,圖像的目視效果如圖6所示?？梢钥闯?經(jīng)大氣校正后,與原始圖像相比,圖像的清晰度得到了明顯改善,對比度也得到了提高。

圖6 大氣校正前后對比圖

2)反射率對比。為評價(jià)本文方法對實(shí)現(xiàn)典型地物大氣校正的效果,分別比較了植被、裸土、水體3種典型地物經(jīng)本文方法和FLAASH大氣校正后的光譜曲線,并計(jì)算了相應(yīng)的RMSE、R2和NDVI,如圖7、表3所示。其中圖7中的光譜曲線由于1 400 nm和1 900 nm附近為強(qiáng)大氣吸收區(qū)域,光譜信號較弱,將其數(shù)值進(jìn)行了去除。

由圖7可以看出,經(jīng)過本文方法大氣校正后,3種典型地物均能夠較好地符合其光譜特性,且與經(jīng)FLAASH校正后的光譜曲線具有良好的一致性,說明本文方法能較好地消除大氣的影響,恢復(fù)地物原有的光譜屬性,得到正常情況下地物的反射率。同時(shí)整幅圖像和3種不同典型地物的RMSE和R2值均保持在較高的水平,說明本文方法的大氣校正效果能夠達(dá)到與FLAASH相同的水平,對實(shí)現(xiàn)不同地物的大氣校正具有一定的普適性,本文算法對實(shí)現(xiàn)GF-5 AHSI圖像的大氣校正具有一定的適用性。

圖7 不同地物大氣校正反射率對比圖

表3 本文方法與FLAASH曲線評價(jià)指標(biāo)

同時(shí)對比大氣校正前后NDVI比值的變化可以發(fā)現(xiàn),兩種方法均有效校正了大氣導(dǎo)致的NDVI比值的降低,本文方法與FLAASH方法NDVI比值變化一致,誤差在±0.02之內(nèi),同時(shí)對植被NDVI比值的校正幅度最大,由0.51增加到了校正后的0.74,說明本文方法能夠有效還原地物原有的NDVI值。

3)分波段對比。為進(jìn)一步評價(jià)本文方法大氣校正的效果,本文對GF-5各個波段的大氣校正效果進(jìn)行了比較,如圖8所示?？梢钥闯?本文方法在絕大多數(shù)波段上與FLAASH保持了較好的一致性,但在1 000 nm、1 400 nm和1 900 nm附近波段誤差較大,分析原因在1 000 nm附近波段主要由于采用MDOTRAN完成大氣校正后數(shù)據(jù)本身波動就較大,FLAASH對此進(jìn)行了光譜平滑處理,消除了這些波動,本文方法未考慮平滑處理,因此導(dǎo)致誤差較大,而1 400 nm和1 900 nm附近波段則是由于位于強(qiáng)水汽吸收波段附近,受到水汽影響較大,因此大氣校正效果較難控制。

3 結(jié)束語

本文基于MODTRAN大氣輻射傳輸模型構(gòu)建查找表,逐像元完成能見度和水汽含量的反演,并實(shí)現(xiàn)了對GF-5 AHSI高光譜圖像的逐像元大氣校正,以經(jīng)FLAASH大氣校正后的圖像作為參考,對本文所提出方法的結(jié)果進(jìn)行了驗(yàn)證。經(jīng)對比可以發(fā)現(xiàn),所提出的方法能夠有效實(shí)現(xiàn)逐像元的能見度和水汽含量反演,反演效果較好。同時(shí),該方法有效消除了大氣對輻射傳輸產(chǎn)生的畸變效應(yīng),較好地恢復(fù)了地物原有的光譜屬性,改善了圖像質(zhì)量。與FLAASH大氣校正相比,該方法能夠?qū)崿F(xiàn)圖像的自動化、批量化大氣校正,為提高GF-5圖像的業(yè)務(wù)化大氣校正提供了一定參考。

由于沒有真實(shí)地物實(shí)測光譜數(shù)據(jù),本文對所提方法大氣校正精度的驗(yàn)證存在一定的困難,同時(shí)本文在參數(shù)反演過程中沒有考慮到水汽含量和氣溶膠對輻射傳輸?shù)膮f(xié)同影響,這也給大氣校正帶來了一定的誤差。因此,下一步需針對這兩個方面繼續(xù)進(jìn)行深入研究,進(jìn)一步提升對GF-5 AHSI高光譜圖像的大氣校正效果。