環(huán)境減災(zāi)二號A/B衛(wèi)星高光譜數(shù)據(jù)反演與質(zhì)量提升

馮向朋 陳鐵橋,2 張耿 王爽 劉學(xué)斌

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

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

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

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

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

1.1 高光譜成像儀

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

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

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

圖2 高光譜成像儀時空聯(lián)合調(diào)制獲取干涉信息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ù)的反演過程實際上是目標(biāo)圖譜數(shù)據(jù)立方體的重構(gòu)過程，過程如圖3所示[4]。對于理想干涉圖，僅需要進(jìn)行實線框中描述的點干涉圖提取和傅立葉變換就可以反演出高光譜數(shù)據(jù)。虛線框所代表的流程是在數(shù)據(jù)并非理想時的預(yù)處理和修正過程。

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

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

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

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

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

(1)

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

(2)

式中：L為最大光程差。

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

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

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

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

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

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

干涉數(shù)據(jù)校正聯(lián)合使用實驗室定標(biāo)法和在軌統(tǒng)計法進(jìn)行校正。探測器響應(yīng)不一致性校正系數(shù)通過實驗室定標(biāo)的方式獲得：去掉前置光學(xué)系統(tǒng)，僅留下探測器對太陽模擬器不同亮度、不同增益數(shù)據(jù)進(jìn)行采集，然后計算得到探測器響應(yīng)不一致性定標(biāo)系數(shù)C1。光學(xué)系統(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ù)和光學(xué)系統(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 復(fù)原后高光譜數(shù)據(jù)不一致性校正

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

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

(7)

式中：Dorg(k,l)為正整數(shù)化的第k個譜段第l行值；Dfloat(k,l)為復(fù)原后第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ù)值在當(dāng)前樣本中的個數(shù)；M(k,l)為第k個譜段第l行所有數(shù)據(jù)在當(dāng)前統(tǒng)計樣本中的總數(shù)。

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

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

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

(9)

對再次正整數(shù)化數(shù)據(jù)進(jìn)行直方圖統(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行的累計直方圖。

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

(12)

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

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

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

(13)

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

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

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

3 試驗結(jié)果及分析

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

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

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

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

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

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

光譜曲線準(zhǔn)確度計算公式為

(14)

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

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

圖8 敦煌定標(biāo)場Fig.8 Dunhuang calibration field

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

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

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

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

圖12 環(huán)境減災(zāi)二號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)境減災(zāi)二號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)計結(jié)果。其中：VNIR通道數(shù)據(jù)光譜準(zhǔn)確度最低為93.72%，最高可達(dá)到97.26%；SWIR通道數(shù)據(jù)最低為92.06%，最高可達(dá)到98.38%。

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

4 結(jié)束語

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