基于輻射不變點的高光譜遙感傳感器交叉定標方法

李映潭 楊遼 王杰

(1 西華師范大學 地理科學學院,四川南充 637002)(2 中國科學院新疆生態與地理研究所,烏魯木齊 830011)

對于遙感成像測控系統,數據處理不僅要關注地物的紋理和顏色等幾何特征,還需要考慮物體的反射和輻射特性[1]。輻射定標是進行大氣校正的關鍵一環,因此獲取逐波段的輻射標定系數(增益與偏置)是光學影像預處理的必要工作。光學遙感影像輻射定標是將原始的像元亮度(DN)值轉化為輻射值的過程[2]。衛星發射前會進行大量定標實驗和飛行實驗,并在升空后進行在軌校準。這些步驟有助于建立影像DN曲線與地表輻射值曲線之間的線性關系,最終將原始的DN值精確轉換為輻射值。受圖像傳感器的非線性響應、太陽輻照度的不穩定以及傳感器屬性影響[3],導致遙感傳感器的測量值難以準確描述真實地物的輻射特性。因此,為了正確描述地表輻射特征,獲得準確的輻射值,需要對標定系數進行精確計算。高光譜遙感影像波段眾多,出現了圖譜合一的特性[4-7],因此,對其輻射標定的精度要求較高。有研究指出資源一號02D衛星高光譜數據傳感器中的120個大氣窗口波段在輻射定標方面平均存在6.18%的誤差,其中有12個波段的誤差超過了7%,誤差滿足遙感衛星應用的精度標準[8],但資源一號02D號衛星提供的標定系數僅包括增益系數,缺乏偏置系數。在實際應用中,使用默認定標參數完成影像輻射定標后,經大氣校正的影像反射光譜曲線可能出現一些偏差。

遙感傳感器交叉定標的基本原理是在待定標衛星和已定標衛星在軌道運行過程中,對同一區域進行遙感成像。通過比較兩者的成像結果并建立關系,實現對待定標衛星傳感器的輻射定標[9-12]。對于交叉定標的研究,文獻[13]最早基于水衛星(Aqua)與美國國家海洋和大氣管理局16號衛星(NOAA-16)數據提出了基于同步天底過境(SNO)的交叉標定方法,該方法嚴格限制了觀測時間與成像角度,且SNO的點多為極地區域,實用性與代表性不足;故文獻[14]后續提出來拓展的SNO方法,將適用范圍擴展到了中低緯度地區;文獻[15]結合雙向反射分布函數(BRDF)模型與6S輻射傳輸模型實現了較大時間差異與觀測天頂角差異的衛星數據交叉標定;文獻[16]通過輻亮度調整因子和反射率調整因子補償高分一號(GF-1)衛星與參考衛星之間的時間差異、天頂角的差異以及光譜響應差異,從而提高了交叉標定精度;文獻[17]以準不變定標場(PICS)交叉點,降低了標定不確定性。

受圖像傳感器的非線性響應、太陽輻照度不穩定等因素影響,遙感傳感器測量值難以準確描述地物真實的輻射特性。因此,精確獲取標定系數成為獲取地表輻射特征的關鍵。交叉定標的實施成本低、頻次高,已成為驗證光學遙感載荷輻射定標系數的常規方法,但多被用于多光譜衛星數據。為解決高光譜衛星傳感器定標參數的精確計算,文章通過交叉定標思路提出了一套完備的定標系數獲取方法。

1 方法原理

對于常規空間遙感衛星輻射定標系數的獲取與校正,需要配合地面光譜儀多點實地測量或專用定標場地,但空間衛星成像時間短,難以實現儀器點測量與遙感面測量的同步。其次,受限于衛星傳感器的空間分辨率,點測量的光譜儀由于混合像元效應,與其對應的像元光譜存在偏差。因此,為實現快速標定系數的獲取與校正,本文基于實施成本低且頻次高的交叉定標思路,引入NASA的地表礦物粉塵源調查(EMIT)高光譜數據,通過在與資源一號02D衛星經過同一區域的軌道交叉點采集數據時,利用空間、光譜與時間分辨率都較高的哨兵二號(Sentinel-2)數據尋找輻射不變點,逐波段建立關系,并得到一套新的定標系數。

選取輻射不變點是兩種數據進行交叉定標的首要工作,通過目視解譯選取輻射不變點費時費力,且目視解譯還存在主觀性的偏差,同時兩種數據分別為DN值與輻射值,無法從光譜維度獲取信息,存在信息上的缺失。基于以上因素,本文使用Sentinel-2衛星的二等A級數據(L2A)數據來選取輻射不變點。

在遙感應用中,光譜角(Spectral angle)是一種比較兩個光譜相似度的指標,通常用于遙感圖像處理和分類[18-20]。給定兩個光譜向量A和B,可以表示為包含多個波段反射率值的向量。光譜角(θ)可以通過以下公式計算

(1)

在本文中,光譜向量A與B分別代表兩種影像同一坐標位置對應的光譜向量,A·B表示向量A和向量B的點積(內積),‖A‖和‖B‖則分別表示向量A和向量B的范數(或長度)。運算結果在0～1,數值越大,表明兩個向量越相似。光譜角的引入使得Sentinel-2影像的輻射不變點選取更具有科學性,且光譜角從角度度量,不受成像時間和太陽輻照度不同引起的亮度和對比度變化的影響。通過編寫交互式數據語言(IDL)代碼對兩景影像逐點求取光譜角,得到一個角度結果矩陣,取矩陣中前100個最大值的坐標作為輻射不變點。相關代碼如圖1所示。

圖1 光譜角匹配的處理代碼

將100個輻射不變點放入EMIT與資源一號02D衛星數據進行位置檢查,共刪除了18個在云區的點,保留82個點用于計算逐波段輻射標定系數。編寫IDL代碼,基于輻射不變點坐標,以資源一號02D衛星數據的DN值為自變量,EMIT數據的輻射值為因變量,逐波段進行線性回歸。標定結果為一個逗號分隔值(CSV)表格,該表格包含了165個波段所對應的偏置、增益相關系數與標準誤差。

對處理完后的影像,以IDL代碼為基礎,結合重新標定的系數,對資源一號02D數據進行逐波段的輻射定標;將得到的輻射值影像導入基于可視化圖像環境(ENVI)軟件的快速視線大氣分析光譜超立方體算法(FLAASH)模塊編寫的批處理程序進行大氣校正處理(相關參數設置見ENVI幫助文檔),得到反射率影像;最后在刪除誤差較大的18個壞波段后,對結果進行3×3模板的光譜平滑。整體的過程框架如圖2所示。

圖2 交叉輻射標定與大氣校正過程框架

2 數據源與數據預處理

資源一號02D衛星于2019年9月12日成功發射,搭載的兩個傳感器分別提供9波段多光譜數據和166波段高光譜數據[21]。多光譜數據的光譜覆蓋范圍為0.452～1.047μm,全色波段的空間分辨率為2.5m,其他8個波段的空間分辨率為10m,成像寬幅為115km。高光譜傳感器的光譜覆蓋范圍為0.40～2.50μm,空間分辨率為30m。在可見光-近紅外波段,光譜分辨率為10nm,短波紅外波段為20nm,成像寬幅為60km。衛星的重訪周期為3天,全球覆蓋周期為55天。本文使用的是L1級產品數據,該數據級別為未地形校正的DN值影像,其中影像主文件為地理標記圖像文件格式(GeoTiff)數據文件,同時附帶可拓展標記語言(xml)元數據說明文件、有理多項式系數模型(RPC)校正參數文件、預覽圖文件、覆蓋區域矢量文件、觀測幾何角度文件以及定標系數文件。其中,影像數據文件按波段分別存為可見近紅外譜段(VN)與短波紅外譜段(SW)兩個GeoTiff文件。

EMIT是部署在“國際空間站”上的傳感器,于2022年8月開始測量,其主要任務是監測地表礦物,特別是地球干旱塵埃源覆蓋區域,在北緯52°和南緯52°之間的陽光照射區進行礦物測量[22]。EMIT Level 2A數據以網絡通用(NC)數據格式存儲,包括3個60m空間分辨率的產品:反射率影像(EMIT_L2A_RFL)、反射率不確定性影像(EMIT_L2A_RFLUNCERT)和反射掩膜影像(EMIT_L2A_MASK)。反射率產品包括285個波段,光譜覆蓋范圍為0.381～2.493μm,光譜分辨率約為7.5nm;反射率不確定性包含了逐像元、逐波段的反射率不確定估計值(見圖3)。

圖3 資源一號02D影像DN值與EMIT影像輻射值

本文使用的是一景2023年5月29日的EMIT數據與一景2021年6月6日的資源一號02D衛星影像。資源一號02D衛星影像的標準處理流程包括數據打開、輻射定標、大氣校正以及正射校正。本文使用該數據的DN值,因此無需進行輻射定標和大氣校正:①數據打開,由于ENVI不支持直接讀取該衛星數據,需要安裝國產衛星支持工具插件,并選擇.xml文件打開數據,該軟件會自動打開可見光-近紅外和短波紅外兩個圖像文件,并自動合并為166個波段的圖像文件,同時打開RPC和其他元數據文件,賦予影像基礎的位置信息;②正射校正,本文使用ENVI的RPC Orthorectification Using Reference Image工具,以同一覆蓋區域的Sentinel-2數據作為參考影像,30m的高級空間熱發射與反射輻射計(ASTER)數據作為參考數字高程模型(DEM)。處理后,得到坐標系為通用橫軸墨卡托投影(UTM)43N、空間分辨率30m的影像。最后為了統一空間分辨率,將影像重采樣到60m。

本文基于IDL程序編寫了EMIT數據預處理的批處理程序,包括以下步驟:①NC格式數據讀取,通過NC文件ID以及目標參數的ID讀取影像的相關參數并存儲為變量;②地理查找表(GLT)地理編碼,在第一步中讀取的數據包括兩個波段變量,分別逐像元記錄了經度和緯度,通過地理編碼將影像投影為經緯度;③增加元數據,最后影像的輸出格式為ENVI標準格式,通過添加頭文件信息的方式添加影像的中心波長、全波半波寬以及地圖信息;④刪除壞波段,通過對照光譜曲線,刪除誤差的波段;⑤影像重投影,將EMIT的經緯度投影轉換為Sentinel-2同一投影帶的UTM投影,確保二者地理位置完全對應,將兩種數據地理位置鏈接后,通過影像對照發現無需進行額外的影像配準。

經裁剪后,EMIT共285波段,資源一號02D衛星共166波段,為了使兩種數據的波段相對應,以資源一號02D衛星數據的光譜為參考,將EMIT數據的光譜重采樣至166波段。但兩種數據波段覆蓋范圍存在差異,光譜重采樣后的EMIT數據最后一個波段為無值,故完全對應的波段數為165。

于ESA的哥白尼數據中心下載了兩景與兩種數據時間各自相近的Sentinel-2數據[23],由于L2A數據為反射率產品,采用超級分辨率將所有波段空間分辨率增強至10m,然后將所有像元除以10000.0轉換為反射率,再重采樣至待定標數據分辨率后裁剪到研究區域,經對照二者無需再進行地理配準,可直接作為光譜角計算輸入(見圖4)。

圖4 兩個時間節點的Sentinel-2影像

3 結果與精度評價

本文使用標準誤差(Se)和相關系數(R)來評價標定結果,標準誤差是用于衡量測量值與真實值之間的偏差程度的指標,數值越小表示標定結果的精度越高;相關系數則用于衡量兩個變量之間的線性關系,數值范圍在-1～+1,越接近1表示標定結果與真實值之間的關聯性越強。165個波段評估結果及選取的輻射不變點如圖5所示。

圖5 自動選擇的樣本點以及各波段的標準誤差與相關系數

為了展示逐波段輻射標定的精度,繪制了9個常用波段的回歸結果,如圖6所示,標注了對應的相關系數、標準誤差與中心波長。

注:圖中R為相關系數;Se為標準誤差;λ為中心波長。

最終的成果影像與樣例光譜如圖7所示。在影像上均勻選取了10個點作為示例點,并展示了它們對應的輻射值曲線與反射率曲線。通過分析光譜曲線可以得知經過該方法定標的影像能夠獲取高質量的反射率光譜曲線。在影像上均勻的選擇示例點既能確保它均勻的分布在整個影像區域,又能充分代表不同地物類型和場景,這樣的選擇有助于驗證交叉定標的普適性和可靠性。光譜曲線的分析揭示了影像的輻射值與反射率之間的關系,對于提高影像數據的準確性和可信度具有重要意義,尤其是在定量遙感研究方面。

圖7 經校正后的反射率影像以及各樣點輻射值曲線與反射率曲線

通過圖6分析可知,有小部分波段存在較大誤差且相關性較低,研究表明:前幾個波段的誤差主要受到大氣的紫外吸收效應影響,而在1374nm附近與1845nm附近的波段主要受到水汽吸收的影響。[8]

700～1000nm的近紅外區域雖然相關性較高,但標準誤差較大,原因在于近紅外波段對大氣參數具有敏感性,這是大氣中的水汽和其他氣體對這個波段的輻射吸收特性。這部分數據涉及到資源一號02D衛星數據的VN與SW數據重疊部分,由圖8分析可知:以VN與SW兩種成像儀于1040nm的散點為例,其分布均滿足線性規律,但二者DN值大小存在數量級上的差別,故考慮為兩種傳感器的設計方案不同導致的該誤差,且該范圍SW成像儀的回歸結果均優于VN成像儀,故在后續處理流程中刪除了重疊部分的VN波段(即圖5中標準誤差曲線峰值部分)。

注:圖中R為相關系數;Se為標準誤差;λ為中心波長。

4 結束語

由于大氣狀況的變化以及遙感傳感器自身損耗等因素,傳感器的標定系數存在不可避免的波動,這種波動直接影響輻射定標的準確性,從而對大氣校正產生不良影響,導致反射率出現偏差,最終影響后續數據的處理與分析。因此,定標系數需要及時更新。傳統的輻射傳輸模型獲取標定系數需要選擇適當的模型,并且需要結合大量地表觀測數據與大氣參數。另外,通過光譜儀實測光譜來進行回歸分析獲取定標系數也需要大量實測數據[24],并涉及到多個地點和時間點的全面系統的數據收集,因此需要大量的人力和時間成本。本文提出的交叉輻射標定方法在執行上簡單快捷,該方法引入了NASA高精度標定的EMIT數據作為參考,并通過IDL代碼實現了大部分流程的自動化,相較于傳統方法,這一套流程框架能夠在短時間內完成數據處理,快速生成高質量的高光譜反射率產品。迄今為止,EMIT已經成功運行了一年多,并積累了大量的數據資源,豐富的數據儲備為交叉輻射標定和后續的科學研究提供了數據支撐。在今后進一步的研究中,將會考慮對不同月份的輻射標定進行詳細探討,時間維度的分析能夠更好地了解太陽輻射在不同季節和時間節點上對影像輻射標定系數的影響。通過深入研究這些時空變化,能夠更準確地把握輻射標定的動態特性,為數據的高質量應用提供更為科學的依據。這一流程框架不僅為當前的科研數據處理工作提供了新思路,也為未來EMIT數據的更全面利用和優化運營提供了有益的方法參考。在高光譜交叉輻射標定方面的深入研究進一步提升了資源一號02D衛星的數據質量,使其在遙感領域具有更廣泛的應用。