環境減災二號A/B衛星高光譜數據反演與質量提升

馮向朋 陳鐵橋,2 張耿 王爽 劉學斌

(1 中國科學院西安光學精密機械研究所，西安 710119)(2 中國科學院大學，北京 100049)

環境減災二號A/B衛星高光譜成像儀是面向“十二五”減災與環境應用需求，研制的一套星載寬覆蓋、高分辨率高光譜成像儀民用載荷。其接替超期服役12年的環境減災一號A衛星超光譜成像儀，通過雙星同軌組網，在全球目標區獲取高光譜遙感圖像數據，用于支持我國環境監測、防災減災等業務工作，同時為國土資源、水利、農業、林業、地震等多個領域提供衛星數據資源支撐和業務化應用服務。

面向上述領域應用研究，為定量分析提供可靠的數據來源，需要對高光譜遙感圖像數據進行相對輻射校正，實現數據質量提升，使得最終生成的數據具備良好的空間信噪比和光譜準確度。相對輻射校正，尤其是遙感數據的相對輻射校正，都是以校正傳感器響應度差異為手段來達到消除條帶效應和提高數據輻射質量的目的[1]。傳統的高光譜遙感圖像數據校正，主要分為實驗室定標法和在軌統計法2種。實驗室定標法采用均勻光源完成對光學系統的不一致標定；在軌統計法基于大量在軌數據統計結果與規律，完成相對輻射校正。

環境減災二號A/B衛星高光譜成像儀采用大孔徑靜態干涉高光譜成像技術，使用時空聯合調制的方式，直接獲取到的數據為不同地物不同光程差下的原始數據。如果使用傳統的實驗室定標法，均勻光源經過其光學系統后，被探測器捕獲的是帶有干涉條紋的數據，不能直接生成相對輻射定標系數。而且，其數據產品并不是傳感器輸出的原始干涉數據，是經過干涉數據反演得到的光譜數據。因此，其條帶噪聲的去除不僅依賴于對傳感器獲取數據進行響應不一致性校正，還需要校正經過復雜數學計算和能量分離之后光譜數據的不一致性，這樣才能為高光譜遙感圖像數據定量分析提供可靠的數據來源。復雜計算和能量分離帶來的不一致性，不能使用實驗室定標法完成，只能使用在軌統計法進行。但是，如果對浮點型數據進行諸如直方圖統計等在軌數據分析，現有計算機是不可能完成的。為了提升環境減災二號A/B衛星高光譜數據質量，本文綜合使用實驗室定標法和在軌統計法，對干涉數據和復原后光譜數據進行不一致性校正，從而有效地去除了圖像條帶，并保持了良好的光譜準確度。

1 高光譜數據反演

1.1 高光譜成像儀

環境減災二號A/B衛星高光譜成像儀基于大孔徑靜態干涉光譜成像技術原理，主要包括擺鏡、干涉儀、傅氏鏡和探測器[2]。其基本組成如圖1所示。

圖1 高光譜成像儀基本組成Fig.1 Basic component of hyperspectral imager

高光譜成像儀通過時空聯合調制的方式，在同一時刻獲得全部目標點各自不同光程差位置上的干涉信息，再利用衛星的移動和時間維的采樣共同獲取完整的干涉信息[3]。其具體過程如圖2所示。

圖2 高光譜成像儀時空聯合調制獲取干涉信息Fig.2 Interferometric information obtained from hyperspectral imager by combined space-time modulation

在圖2中，第1幀圖像包含了t1時刻第1行地物探測器第1個光程差位置上的干涉信息、第2行地物探測器第2個光程差位置上的干涉信息，一直到第m行地物探測器第m個光程差位置上的干涉信息。其中：m為探測器干涉方向像元數目，也就是探測器所有光程差位置的個數。以此類推，第i幀(i=1，2，3，4，…)圖像包含了ti時刻第i行地物探測器第1個光程差位置上的干涉信息、第i+1行地物探測器第2個光程差位置上的干涉信息，一直到第i+m-1行地物探測器第m個光程差位置上的干涉信息。本文稱這種數據為大孔徑靜態干涉光譜成像(LASIS)數據。

1.2 高光譜數據反演

高光譜數據的反演過程實際上是目標圖譜數據立方體的重構過程，過程如圖3所示[4]。對于理想干涉圖，僅需要進行實線框中描述的點干涉圖提取和傅立葉變換就可以反演出高光譜數據。虛線框所代表的流程是在數據并非理想時的預處理和修正過程。

圖3 高光譜數據反演一般過程Fig.3 General process of hyperspectral data recovering

由高光譜成像原理可知，要獲取同一地物所有光程差下的干涉數據，需要將原始的LASIS數據進行抽取，即要獲取第i行地物的完整干涉數據，就必須抽取第i幀數據的第1個光程差位置、第i+1幀數據的第2個光程差位置，一直到第i+m+1行的第m個光程差位置，就可得到1行地物完整的干涉數據，具體如圖4[5]所示。

圖4 地物目標點完整干涉信息獲取方法Fig.4 Acquisition method of integral interference information of ground object target point

高光譜成像儀基本原理建立在光譜與干涉信息之間具有傅立葉變換關系的基礎上，儀器通過對干涉信息的采集與變換間接獲取光譜信息。假設某一點干涉數據為I(x)，其中，x為沿飛行方向某列探測單元的位置；其復原后點光譜為B(v)={v1,v2,…,vr}，其中，v1為起始波長，vr表示終止波長。I(x)和B(v)相互之間的變換關系可以通過式(1)和式(2)來表示,其中，d為干涉儀剪切量，f為傅氏鏡的焦距。

某列探測單元位置x處的干涉強度可表示為[6]

(1)

提取出某一地物目標點的完整干涉圖后，經過傅立葉逆變換可以得到該目標點復原后點光譜[6]為

(2)

式中：L為最大光程差。

在干涉儀出現伊始，高光譜數據反演只是通過如式(2)所示簡單的傅立葉變換實現[7]。隨著干涉光譜技術的發展，相關的數據處理技術不斷出現，使得高光譜數據處理的精度不斷提高，數據處理流程也不斷完善。與高光譜數據處理技術相關的各種誤差修正方法及數據處理方法相繼被提出，從不同的角度提高了數據處理的精度，完善了高光譜數據處理技術，逐步形成一套通用的光譜數據處理流程(如圖3所示)。其計算復雜度高，算法參數多樣，可適應不同場景與應用下的光譜復原需求。

2 數據質量提升

在高光譜數據反演過程中，暗電流去除、直流分量去除、壞像元修正、切趾、相位修正、絕對輻射校正環節都有較為固定的方法。例如：切趾過程中較為常用的三角窗、海明窗、海寧窗、Blackman、Norton-Beer等切趾函數[8-10]；相位修正過程中經典的Forman法[11]和Mertz法[12]；去直流分量過程中的差分法[13]、擬合法[14]。

相對輻射校正作為校正探測器響應不一致性、消除條帶噪聲的直接手段，在干涉數據復原過程中有著重要的作用。但是，由于干涉儀的存在，傳統的干涉數據相對輻射校正多采用實驗室定標法，去掉前置光學系統(擺鏡、干涉儀、傅氏鏡)，直接獲取探測器不一致性校正系數。這種方式忽略了擺鏡每點的反射率不一致性、干涉儀半透半反性能問題和傅氏鏡干涉效率問題。

為了提升復原后高光譜數據質量，本文對干涉數據和復原后高光譜數據進行不一致性校正，即利用實驗室定標法和統計方式對在軌原始干涉數據直接進行不一致性校正，從而提升復原后高光譜圖像的空間一致性。復原后高光譜數據不一致性校正，是對復原后各個譜段數據進行直方圖統計，然后進行直方圖匹配，得到直方圖查找表進行校正。

2.1 干涉數據不一致性校正

干涉數據校正聯合使用實驗室定標法和在軌統計法進行校正。探測器響應不一致性校正系數通過實驗室定標的方式獲得：去掉前置光學系統，僅留下探測器對太陽模擬器不同亮度、不同增益數據進行采集，然后計算得到探測器響應不一致性定標系數C1。光學系統的不一致性通過在軌統計法獲得，即：統計不同時相下數據，然后通過計算得到每個光程差的相對輻射校正系數C2。具體計算方法如下。

獲取多軌原始干涉數據，累加之后求平均值，得到平均干涉數據，如式(3)所示。

(3)

式中：IAVG(p,q)為第p個光程差第q行的平均值；Ia(p,q)為第a幀數據第p個光程差第q行的值；N為干涉數據幀數。

每個光程差下的累計平均值為

(4)

式中：IAVG(p)為第p個光程差下的累計平均值；W為干涉數據幀的列數。

以單行光程差平均值除以行積累平均值，得到相對輻射校正系數C2。其中，C2(p,q)為相對輻射校正中第p個光程差第q列的系數。

(5)

將探測器不一致性校正系數和光學系統不一致性系數合并，得到干涉數據相對輻射校正系數C。

C(p,q)=C1(p,q)C2(p,q)

(6)

式中：C(p,q)為第p個光程差第q列的系數。

2.2 復原后高光譜數據不一致性校正

統計誤差及切趾、相位修正、逆傅立葉變換等處理過程中產生的計算誤差，導致干涉數據相對輻射校正之后復原數據依然存在條帶。而且，這些復雜計算過程產生的不一致性難以通過單一的線性關系進行校正。因此，本文通過直方圖匹配的統計學方法[15]對復原后光譜數據進行校正。具體步驟如下。

光譜復原后，因經過復雜計算，數據為浮點FMax，將所有譜段的數據都正整數化到[1,Z]區間，對區間進行直方圖統計。歸一化過程為

(7)

式中：Dorg(k,l)為正整數化的第k個譜段第l行值；Dfloat(k,l)為復原后第k個譜段第l行原始浮點值；round表示臨近取整函數。

統計各個譜段灰度分布直方圖概率密度，見式(8)。

PDNorg(k,l)=nDNorg(k,l)/M(k,l)

(8)

式中：PDNorg(k,l)為正整數化后第k個譜段第l行數值的概率密度；nDNorg(k,l)為正整數化后第k個譜段第l行數值在當前樣本中的個數；M(k,l)為第k個譜段第l行所有數據在當前統計樣本中的總數。

找到能量響應最高的譜段，對概率密度分布函數進行95%的截斷，去掉概率分布的拖尾部分，取截斷之后的最大值。截斷之后的灰度值分布個數如圖5所示。

圖5 95%截斷后能量響應最高譜段的灰度值分布Fig.5 Digital number distribution of 95% truncated for the highest energy response spectral band

假設此時所有譜段中最大值不超過Zsta(圖5中最大值不超過2000)，那么要重新將復原后浮點型數據正整數化到[1,Z]區間內，即

(9)

對再次正整數化數據進行直方圖統計，獲取各譜段灰度分布直方圖的概率密度函數為

PD(k,l)=nD(k,l)/R

(10)

式中：PD(k,l)為第k個譜段第l行地物灰度值為D的概率密度函數；nD(k,l)為第k個譜段第l行地物灰度值為D的個數；R為統計地物的行數。

求取各像元的累計直方圖的概率密度函數為

(11)

式中：SD(k,l)為灰度值從1到D在第k個譜段第l行的累計直方圖。

分別對各譜段的累計直方圖進行平均計算，得到不同譜段下的平均累計直方圖，見式(12)。

(12)

式中：SD(k)為譜段k從灰度值1到D的累計直方圖概率密度函數。

分別對各譜段各行的累計直方圖和對應譜段的平均累計直方圖進行匹配計算，找出各譜段各行累計直方圖和對應譜段平均累計直方圖的灰度映射關系，見式(13)。

D(k,l)=gSg(k)≤SD(k,l)≤Sg+1(k)

(13)

式中：Sg(k)為譜段k在灰度值為g時的平均累計直方圖概率密度函數。

式(13)表示直方圖匹配關系，即：當Sg(k)≤SD(k,l)≤Sg+1(k)時，D(k,l)的值匹配為g。

根據式(13)建立灰度值查找表，在上述統計過程中，本文選取了大量的非飽和干涉數據進行反演和計算，包括裸地、植被、水體、山地、雪地、建筑物等。查找表建立完成后，進行校正。首先使用式(9)進行正整數化，然后對照直方圖查找表進行灰度值映射。如果被查找灰度值不在[1，Z]區間，將其設置為Z。

3 試驗結果及分析

3.1 干涉數據校正結果

進行干涉數據校正之后，直觀結果如圖6所示，干涉數據列方向條帶與分塊現象消除明顯。

圖6 干涉數據相對輻射校正前后對比Fig.6 Comparison of interference data before and after radiometric correction

3.2 高光譜數據校正結果

如圖7所示，高光譜數據相對輻射校正后，其數據質量大大提升；而且，條帶噪聲得到了有效的去除，保持了很好的光譜準確度。

直方圖匹配校正不是基于線性關系的校正，可能會對光譜準確度產生一定的影響。因此，本文使用敦煌定標場和高亮地區實測數據進行分析。分析結果表明：直方圖匹配不會對光譜曲線準確度造成過大影響。

光譜曲線準確度計算公式為

(14)

圖7 光譜數據相對輻射校正前后圖像質量對比Fig.7 Image quality comparison of spectral data before and after radiometric correction

選取敦煌定標場定標區域對光譜準確度進行測試，將L1A級高光譜數據進行絕對輻射校正，得到入瞳輻亮度數據，將定標區域實測圖像處理質量標準數據反算到入瞳輻亮度，對比2組數據進行計算，得到光譜曲線準確度。測試區域如圖8所示，選取的圖像區域如圖8中紅色矩形框所示，其中，可見光近紅外(VNIR)通道數據取10×10像元，短波紅外(SWIR)通道數據取5×5像元(紅外去除水汽吸收譜段)。環境減災二號A/B衛星敦煌定標場定標數據光譜準確度曲線如圖9和圖10所示。從圖9和圖10可以看出：環境減災二號A/B衛星的高光譜成像數據光譜準確度已經分別達到97%以上和92%以上。

圖8 敦煌定標場Fig.8 Dunhuang calibration field

圖9 環境減災二號A衛星高光譜數據與實測數據對比Fig.9 Comparison between HJ-2A satellite hyperspectral data and measured data

圖10 環境減災二號B衛星高光譜數據與實測數據對比Fig.10 Comparison between HJ-2B satellite hyperspectral data and measured data

選取敦煌定標場北部高亮區域(見圖11)不同時相高光譜數據，進行絕對輻射校正得到入瞳輻亮度數據，與實際測量數據進行對比。環境減災二號A/B衛星光譜準確度結果如圖12和圖13所示，與實際測量數據的光譜相似度分別達到了93.72%和94.56%。

圖11 高亮區域Fig.11 Highlighted region

圖12 環境減災二號A衛星VNIR通道高光譜數據與實測數據對比Fig.12 Comparison between HJ-2A satellite VNIR channel hyperspectral data and measured data

圖13 環境減災二號B衛星VNIR通道高光譜數據與實測數據對比Fig.13 Comparison between HJ-2B satellite VNIR channel hyperspectral data and measured data

表1為上述試驗數據的統計結果。其中：VNIR通道數據光譜準確度最低為93.72%，最高可達到97.26%；SWIR通道數據最低為92.06%，最高可達到98.38%。

表1 光譜曲線準確度分析結果(絕對光譜準確度)Table 1 Results of spectral curve accuracy analysis (absolute spectral accuracy)

4 結束語

干涉光譜成像數據條帶噪聲的去除是一個系統性問題，不能僅僅依靠單一的相對輻射校正手段，需要綜合使用實驗室定標法和在軌數據統計法。本文所述數據質量提升方法，通過降低原始干涉數據響應不一致性噪聲和反演后光譜數據中的響應不一致性噪聲和計算噪聲，具備良好的相對輻射校正效果，并保持了92%以上光譜準確度。在軌數據特征統計過程中，僅僅剔除了少量無效數據(云、飽和數據等)生成干涉數據相對輻射校正系數與復原后光譜數據直方圖查找表，有效地解決了載荷光學系統星地環境不一致帶來的差異性。綜合使用不同的相對輻射校正方法，將原始干涉數據、實驗室定標系數校正后干涉數據、干涉數據在軌狀態，反演后光譜數據在軌統計結果有機結合，實現更好的干涉光譜成像數據質量提升效果。環境減災二號A/B衛星高光譜數據質量方法，給出了從機理上校正不同光程差干涉數據的理論，利用直方圖統計法和歸一化方法將反演后浮點數量化到整數范圍，并切除拖尾，使得基于復雜校正原理的反演后浮點高光譜數據校正成為可能。為了完善本文所述復原后相對輻射校正機理，進一步提升光譜準確度，后續將改進復原后直方圖匹配方法為線性修正模式。