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

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

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

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

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

1 高光譜數(shù)據(jù)反演

1.1 高光譜成像儀

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

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

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

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

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

1.2 高光譜數(shù)據(jù)反演

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

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

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

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

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

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

(1)

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

(2)

式中：L為最大光程差。

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

2 數(shù)據(jù)質量提升

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

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

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

2.1 干涉數(shù)據(jù)不一致性校正

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

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

(3)

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

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

(4)

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

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

(5)

將探測器不一致性校正系數(shù)和光學系統(tǒng)不一致性系數(shù)合并，得到干涉數(shù)據(jù)相對輻射校正系數(shù)C。

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

(6)

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

2.2 復原后高光譜數(shù)據(jù)不一致性校正

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

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

(7)

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

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

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

(8)

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

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

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

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

(9)

對再次正整數(shù)化數(shù)據(jù)進行直方圖統(tǒng)計，獲取各譜段灰度分布直方圖的概率密度函數(shù)為

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

(10)

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

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

(11)

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

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

(12)

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

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

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

(13)

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

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

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

3 試驗結果及分析

3.1 干涉數(shù)據(jù)校正結果

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

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

3.2 高光譜數(shù)據(jù)校正結果

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

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

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

(14)

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

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

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

圖9 環(huán)境減災二號A衛(wèi)星高光譜數(shù)據(jù)與實測數(shù)據(jù)對比Fig.9 Comparison between HJ-2A satellite hyperspectral data and measured data

圖10 環(huán)境減災二號B衛(wèi)星高光譜數(shù)據(jù)與實測數(shù)據(jù)對比Fig.10 Comparison between HJ-2B satellite hyperspectral data and measured data

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

圖11 高亮區(qū)域Fig.11 Highlighted region

圖12 環(huán)境減災二號A衛(wèi)星VNIR通道高光譜數(shù)據(jù)與實測數(shù)據(jù)對比Fig.12 Comparison between HJ-2A satellite VNIR channel hyperspectral data and measured data

圖13 環(huán)境減災二號B衛(wèi)星VNIR通道高光譜數(shù)據(jù)與實測數(shù)據(jù)對比Fig.13 Comparison between HJ-2B satellite VNIR channel hyperspectral data and measured data

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

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

4 結束語

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