王宇軒 黃紅蓮 劉曉 孫曉兵 林軍
(1 中國科學院合肥物質(zhì)科學研究院安徽光學精密機械研究所,合肥 230031) (2 中國資源衛(wèi)星應用中心,北京 100094)
環(huán)境減災二號A/B衛(wèi)星于2020年9月27日成功發(fā)射,目前已正常在軌運行1年多。其裝有16 m相機和大氣校正儀。其中,大氣校正儀通過同步獲取大氣參數(shù),對16 m相機進行大氣校正[1]。
大氣中的氣溶膠顆粒物組分多樣、光學性質(zhì)復雜,且在時域和空域具有高度變化的特征,是大氣中的主要不確定因素[2],單個譜段1 h內(nèi)氣溶膠光學厚度變化量可達到0.15以上。目前,進行大氣校正使用的輔助數(shù)據(jù)多來自中分辨率成像光譜儀(MODIS)產(chǎn)品,需要對MODIS原始產(chǎn)品進行時間匹配和空間匹配、分辨率重采樣及其他處理,這就造成了大氣參數(shù)的離散性,最終導致校正效果不理想。為了得到高精度校正圖像,16 m相機和大氣校正儀共同搭載于環(huán)境減災二號A/B衛(wèi)星進行同時同地觀測,獲取用于大氣校正的時間匹配和空間匹配同步大氣參數(shù)。
因為大氣校正儀可以獲取有效大氣參數(shù),因此本文選擇6S模型作為大氣輻射傳輸模型,實現(xiàn)環(huán)境減災二號A/B衛(wèi)星大氣校正程序[3]。本文通過大氣校正儀數(shù)據(jù)預處理程序?qū)Νh(huán)境減災二號A/B衛(wèi)星大氣校正儀在軌獲取數(shù)據(jù)進行預處理和水汽、氣溶膠參數(shù)同步反演,再基于大氣參數(shù)查找表逐像元對圖像進行校正處理。最后,與敦煌地區(qū)和合肥地區(qū)星地同步實測地物光譜數(shù)據(jù)以及MODIS產(chǎn)品進行對比,驗證了大氣校正技術的精確性。
為排除衛(wèi)星遙感器觀測的非目標物帶來的干擾信息,需要對遙感圖像進行大氣校正,這也是定量遙感的一個必備環(huán)節(jié)[4]。衛(wèi)星遙感器對地面進行觀測時,衛(wèi)星入瞳處接收到的輻亮度可分為3個部分:①目標地表輻射;②大氣程輻射;③目標物周圍環(huán)境反射的太陽輻射引起的地表背景輻射[5]。其中:目標像元對衛(wèi)星遙感器入瞳信號貢獻一般為80%,因此使衛(wèi)星遙感器接收到的目標地物信息真實性較低,需要進行大氣校正,以消除大氣程輻射和地表背景輻射的干擾,從而獲取地物的真實反射率。本節(jié)主要介紹16 m相機和大氣校正儀以及大氣校正的流程和原理。
環(huán)境減災二號A/B衛(wèi)星均裝載了16 m相機。2臺相機設計原理相同,運行軌道完全重合,經(jīng)過同一點時間相隔1天;采用線陣推掃式,共設置5個譜段(如表1所示);采用四鏡頭大視場拼接,1幅圖像像元數(shù)為12 000×12 000×5。大氣校正儀的具體參數(shù)如表2所示,它共設置9個譜段,包括可見光近紅外(VNIR)譜段,針對偏振輻射具備多光譜的探測能力[6]。利用大氣校正儀與16 m相機相互對應譜段的光譜范圍設置和同時同地觀測的特征,再經(jīng)過遙感圖像和輔助數(shù)據(jù)的像元經(jīng)緯度匹配,就可以逐像元獲取同步高精度的水汽和氣溶膠等大氣參數(shù),用于遙感圖像的大氣校正。這一部分主要實現(xiàn)的數(shù)據(jù)處理目標有:地面數(shù)據(jù)預處理,即L0級數(shù)據(jù)到L1級數(shù)據(jù)的生成;反演同步大氣參數(shù),即L1級數(shù)據(jù)到L2級數(shù)據(jù)的生成。
表1 16 m相機參數(shù)Table 1 Parameters of 16m camera
表2 大氣校正儀參數(shù)Table 2 Parameters of PSAC
對于環(huán)境減災二號A/B衛(wèi)星的大氣校正,首先將有云圖像剔除,再使用大氣校正儀同步數(shù)據(jù)根據(jù)經(jīng)緯度信息與L1級圖像進行范圍精匹配。將大氣校正儀同步反演得到的水汽含量和氣溶膠L2級產(chǎn)品參數(shù),分別針對大氣吸收和大氣散射影響作為輸入?yún)?shù)進行大氣校正中的輻射校正部分和目標區(qū)域的鄰近效應校正部分,最終實現(xiàn)圖像的大氣同步校正。處理過程如圖1所示。
圖1 大氣校正過程Fig.1 Process of atmospheric correction
對于本文選擇的6S輻射傳輸模型,首先將地面假設為表面均一、各向同性的朗伯體,則衛(wèi)星遙感器接收到的表觀輻亮度可表示為[7]
(1)
式中:Lp為程輻射;τv為接收到的大氣透射率,即直接透射率和漫透射率之和;ρ為像元反射率;ρ1為地面綜合反射率;Eg為ρ=0時水平地表上的總輻射通量;S為大氣半球反照率。
根據(jù)式(1)可解出地面綜合反射率為
(2)
大氣校正整個過程可分為2步:①根據(jù)大氣校正儀同步反演的水汽含量、氣溶膠參數(shù),以及16 m相機的譜段設置、光譜響應函數(shù)、觀測幾何,使用大氣輻射傳輸模型計算大氣中的透過率、程輻射和L1級產(chǎn)品中所包含的太陽天頂角、太陽方位角、衛(wèi)星天頂角、衛(wèi)星方位角等參數(shù),針對每個譜段構建大氣校正查找表,針對每個像元進行查找,得到輻射程、總透過率、地面輻照度等參數(shù)。②根據(jù)式(2),對16 m相機的L2級輻亮度圖像和①中查詢到的參數(shù),從入瞳處的輻亮度中扣除大氣程輻射,并且去除大氣衰減影響,反演出含鄰近效應的地面綜合反射率。
在成像過程中,衛(wèi)星入瞳處接收到的輻亮度不僅有目標像元的貢獻,更包含周圍自然環(huán)境反射太陽輻射帶來的干擾,這種現(xiàn)象被稱為鄰近效應。在衛(wèi)星定量化遙感過程中,鄰近效應會影響衛(wèi)星對地面目標的遙感能力[8]。目標像元貢獻的這部分輻射信息不包含地面目標物的光譜信息,因此要提高衛(wèi)星圖像質(zhì)量及定量化遙感的精度,就必須去除掉周圍像元的干擾,進行鄰近效應校正。當目標像元與周圍像元反射率較為均一、差別較小時,目標像元反射的太陽輻射和周圍像元反射的太陽輻射近似,鄰近效應的影響相對較小;當目標地表周圍反射率組合較為復雜,目標像元反射的太陽輻射和周圍像元反射的太陽輻射就會有較大差異,鄰近效應的影響也會增大。研究表明:圖像鄰近效應的影響因素還有很多,如觀測波長、衛(wèi)星分辨率、大氣能見度、氣溶膠類型及大氣分子的散射相函數(shù)[9]。鄰近效應中周圍干擾像元的影響可以看作目標表面的輻射場和大氣點擴散函數(shù)的卷積。因此,背景反射率可由周圍像元反射率的點擴散函數(shù)近似求得[10],如式(3)所示。
ρm=?f(r(x,y))·ρ(x,y)dxdy
(3)
式中:ρm為背景反射率;f(r(x,y))為大氣點擴散函數(shù);ρ(x,y)為坐標為(x,y)的點的反射率;x為水平方向距離;y為豎直方向距離。
考慮到求解大氣點擴散函數(shù)的復雜性,一般選擇基于經(jīng)驗的近似方法代替大氣點擴散函數(shù)的精確求解,以離散形式求背景反射率[11],如式(4)所示。
(4)
式中:p(x,y)為鄰近效應中的權重函數(shù),它描述了與目標像元的距離為r(由x,y求出)的點對目標像元的貢獻率;k和l分別為橫向和縱向的像元數(shù)。
為了更好地進行鄰近效應校正,一般采用漫射與直射透過率之比作為鄰近效應校正系數(shù)[12]。本文根據(jù)遙感器空間分辨率越高鄰近效應的影響越強的理論,增加校正系數(shù)以提高校正效果,如式(5)所示。
ρt=ρ1+q(ρ1+ρm)
(5)
式中:ρt為地面真實反射率;q為校正系數(shù)。
為了驗證環(huán)境減災二號A/B衛(wèi)星快速監(jiān)控環(huán)境變化、快速獲取災害圖像的能力,本文選擇火山噴發(fā)和浮塵天氣2個典型災害場景進行對比驗證。數(shù)據(jù)分別來源于2021年2月23日環(huán)境減災二號B衛(wèi)星16 m相機CCD3于意大利埃特納附近(15°9′E,37°6′N)和2021年4月15日環(huán)境減災二號A衛(wèi)星16 m相機CCD3相機于中國北京(115°6′E,40°4′N)拍攝的2組圖像。圖2為大氣校正前后假彩色合成后的典型地物場景對比,由藍譜段(B1)、綠譜段(B2)、紅譜段(B3)3個譜段合成。需要校正的圖像尺寸為192 km×192 km,整個校正過程平均時間為27 min??梢钥闯觯?組對比圖像相較于校正前變得清晰且銳化,對比度增加,信息熵增加,校正效果顯著,能較好地去除衛(wèi)星入瞳處非目標物的干擾,具體參數(shù)信息如表3所示。
圖2 多光譜圖像大氣校正前后對比Fig.2 Comparison of multispectral images before and after atmospheric correction
表3 校正結(jié)果Table 3 Correction results
為驗證校正效果,本文利用MODIS同時同地生產(chǎn)的反射率產(chǎn)品與校正后的圖像進行對比驗證。MODIS地表反射率產(chǎn)品較為常用的有2種,分別是土(Terra)衛(wèi)星生產(chǎn)的MOD09GA和水(Aqua)衛(wèi)星生產(chǎn)的MYD09GA反射率產(chǎn)品[13],由于需要圖像各個譜段對應才能進行對比驗證,選擇將同期過境的MOD09GA產(chǎn)品和MYD09GA產(chǎn)品中的4個譜段組合后與環(huán)境減災二號衛(wèi)星16 m相機前4個譜段實現(xiàn)一一對應。因為MOD09GA和MYD09GA反射率產(chǎn)品空間分辨率為500 m,數(shù)據(jù)格式為HDF5,所以需要進行數(shù)據(jù)預處理。在預處理流程中,首先進行WGS84投影坐標轉(zhuǎn)換,再對圖像重采樣處理,最后進行譜段合成[14]。在預處理后的圖像中,對戈壁灘和植被這2種典型地物的區(qū)域分別隨機選取1個像元,提取2種產(chǎn)品每個譜段在所選位置的反射率值,環(huán)境減災二號衛(wèi)星16 m相機產(chǎn)品和MODIS產(chǎn)品反射率對比如圖3所示,數(shù)據(jù)來源分別為2021年4月4日環(huán)境減災二號A衛(wèi)星16 m相機CCD3拍攝于合肥(118°1′E,31°5′N)和2020年11月27日環(huán)境減災二號A衛(wèi)星16 m相機CCD1拍攝于敦煌(93°8′E,41°1′N)。
從圖3中可以看出:經(jīng)過大氣校正后,環(huán)境減災二號衛(wèi)星16 m相機所選的2種典型地物絕大部分譜段整體與MODIS光譜曲線基本一致,植被由于葉綠素產(chǎn)生的550 nm處的小反射峰和其兩側(cè)460 nm,670 nm的2個吸收帶都能得到較好的體現(xiàn);戈壁灘除B1譜段存在誤差外,其余譜段誤差都相對較小。
圖3 16 m相機產(chǎn)品和MODIS產(chǎn)品反射率對比Fig.3 Reflectance comparison of 16m camera and MODIS products
如圖4所示,分別選取敦煌戈壁灘、合肥小麥田和水體(5像元×4個譜段)作為典型地物對比反射率,可以看出相關系數(shù)達到0.98以上。
圖4 3種典型地物的16 m相機產(chǎn)品與MODIS產(chǎn)品反射率對比Fig.4 Reflectance comparison of three typical features of 16m camera and MODIS products
第1組地面實測數(shù)據(jù)來源于2021年1月27日對敦煌戈壁灘(94°24′E,40°17′N)的測量,氣溶膠模式為沙漠氣溶膠模式。由于敦煌地區(qū)冬季干旱少雨,測量時大氣模式為中緯度冬季大氣,實測時間與衛(wèi)星圖像拍攝時間帶來的反射率誤差可以忽略不計[15];第2組地面實測數(shù)據(jù)來源于2021年3月25日對合肥附近小麥田(117°50′E,31°26′N)的測量。
如表4和表5所示,經(jīng)過校正后,戈壁灘和小麥田反射率更接近地面實測,整體來看2組光譜變化趨勢具有較好一致性,每個譜段反射率誤差值相比于校正前都有明顯減小。敦煌戈壁灘校正后誤差均小于18%;合肥小麥田因植物某些譜段反射率較小,造成相對誤差較大,但誤差仍在可接受范圍內(nèi)。校正后的光譜曲線與地面實測光譜曲線相關性較高,且相關系數(shù)均高于0.95。
表4 合肥小麥田大氣校正前后反射率與地面實測數(shù)據(jù)Table 4 Reflectance of Hefei wheat field before and after atmospheric correction and measurement data
表5 敦煌戈壁灘大氣校正前后反射率與地面實測數(shù)據(jù)Table 5 Reflectance of Dunhuang Gobi before and after atmospheric correction and measurement data
本文對于環(huán)境減災二號衛(wèi)星16 m相機結(jié)合大氣校正儀,基于6S輻射傳輸模型實現(xiàn)了大氣校正,并選擇合肥地區(qū)和敦煌地區(qū)衛(wèi)星圖像根據(jù)MODIS產(chǎn)品和地面實測數(shù)據(jù)對校正后的圖像進行校正效果研究,主要結(jié)論如下。
(1)大氣校正儀能提供較高精度的大氣參數(shù),對提高大氣校正精度具有重要意義。
(2)基于6S模型和大氣校正儀獲取高精度同步大氣參數(shù)完成了大氣校正,并對16 m相機進行逐像元校正。
(3)將校正后的圖像分別與MODIS產(chǎn)品和地面實測地表反射率進行對比,光譜曲線具有較好的一致性,大氣校正效果明顯。