伊丕源， 劉原麟， 武 鼎， 趙英俊， 劉洪成

(核工業(yè)北京地質(zhì)研究院， 遙感信息與圖像分析技術(shù)國家級重點實驗室， 北京 100029)

遙感影像采集受到太陽輻射、大氣條件、地形等多種因素的影響，導(dǎo)致多時相影像中的同類地物光譜特征存在差異，進而會影響到定量產(chǎn)品的一致性[1-2]，因此須要進行輻射歸一化處理。輻射歸一化包括絕對輻射校正和相對輻射歸一化[3]。絕對輻射校正是將傳感器獲得像元灰度值轉(zhuǎn)換成地表反射率，通常是通過大氣輻射校正實現(xiàn)，主要包括統(tǒng)計學模型和基于大氣輻射傳輸理論的模型。此外，為控制和減少不同衛(wèi)星影像之間，由于成像條件、傳感器響應(yīng)等因素造成的差異，許多學者對遙感影像的相對輻射歸一化開展了研究，并提出了多種校正方法[4-5]，基本原理是針對多時相遙感圖像的各個波段分別建立對應(yīng)的波段之間的校正方程，進行輻射歸一化計算處理，從而減少各種影響因素產(chǎn)生的噪聲。目前的相對輻射歸一化校正處理研究主要是針對多時相、多傳感器的衛(wèi)星遙感影像，然而與衛(wèi)星影像不同，航空高光譜遙感(尤其是推掃式傳感器)每條航帶的數(shù)據(jù)采集需要一定的時間過程，且作業(yè)時間內(nèi)太陽輻射、大氣條件處于持續(xù)變化之中，而這些變化必定對傳感器接收的輻亮度信息和反射率反演產(chǎn)生一定的影響。但針對航空光譜遙感的相對輻射歸一化校正方法研究較少。

高精度的反射率光譜反演是基于高光譜圖像開展應(yīng)用研究的前提和基礎(chǔ)[6-7]。首先基于常用的統(tǒng)計學模型開展了大氣輻射校正，進一步分析了時間因子對反射率反演的影響。然后，提出了顧及時間因子的相對輻射歸一化校正方法，從而減少數(shù)據(jù)采集過程中時間變化對同類地物反射率校正結(jié)果的影響。

1 實驗數(shù)據(jù)

實驗區(qū)位于新疆卡姆斯特地區(qū)，2018年10月14日采用塞斯納飛機搭載航空成像光譜儀(compact airborne spectrographic imager, CASI)完成了高光譜影像采集。航飛高度為海拔2 600 m，采集的CASI航空高光譜數(shù)據(jù)波譜范圍：380～1 050 nm(可見光與近紅外波段)，48個波段，空間分辨率為1 m。影像像元值為DN(digital number)值，與輻亮度單位之間的關(guān)系為1 000 DN=1 μW/(cm2·sr·nm)。航線敷設(shè)如圖1所示，每條航線長度為80 km，航帶旁向重疊率約15%，航帶1的采集作業(yè)時間為13：06—13：21，航帶2的采集作業(yè)時間為13：28—13：43。

2 實驗方法

2.1 大氣輻射校正方法

高光譜影像絕對輻射校正是將傳感器獲得的輻亮度值轉(zhuǎn)換為反射率值。目前反射率光譜反演主要有統(tǒng)計學模型和基于大氣輻射傳輸理論的模型等[8]?；诮y(tǒng)計學模型的方法，如經(jīng)驗線性法，需要兩個以上光譜均一、有一定面積大小的目標。實際數(shù)據(jù)采集過程中，在每條航帶中布設(shè)明暗地物是難以實施的，通常是在某一航帶中布設(shè)明暗地物，并以其校正系數(shù)對其他航帶進行大氣輻射校正。基于大氣輻射傳輸?shù)哪Ｐ?，常用的軟件包?S、Modtran、FLAASH等[9-11]。統(tǒng)計學模型選擇常用的經(jīng)驗線性法，并基于大氣輻射傳輸模型對不同時間的大氣下行輻射等參數(shù)進行計算與對比分析。

2.1.1 經(jīng)驗線性法

假定圖像像元DN值與反射率R存在線性關(guān)系為

DN=aR+b

(1)

式(1)中:a、b為校正系數(shù)。通常需要在圖像中兩個以上光譜均一、有一定面積的定標地物，分別為暗目標和亮目標，利用線性回歸[式(1)]建立起兩個目標與對應(yīng)的圖像像元平均值間的關(guān)系，即解算出校正系數(shù)a、b[12]。

實驗過程中，在航帶1的北端、航線2的南端都鋪設(shè)了定標黑白布作為明暗地物，并同步測量了其光譜反射率，然后分別以其為參考，計算對應(yīng)的校正系數(shù)。首先，以航帶1的校正系數(shù)對航帶1、航帶2進行校正，然后以航帶2的校正系數(shù)對航帶2自身進行校正；然后選擇兩航帶重疊區(qū)域的同名地物，對其反射率計算結(jié)果進行對比。

2.1.2 大氣輻射傳輸法

大氣輻射傳輸模型表征了大氣對遙感傳感器接收到的地表輻射信號的影響。對于可見光-短波紅外遙感而言，在假定地表為水平均一朗伯體的情況下，依據(jù)大氣輻射傳輸方程，傳感器接收的輻射信號L可表示為[13-14]

(2)

式(2)中:Lpath為路徑輻射；S為大氣半球反照率；Fd為下行總輻射；Tv為地面像元與傳感器之間的大氣傳輸透過率；ρ為地表反射率。實驗中，選用Modtran軟件，基于兩條航帶各自的獲取時間，計算出航空高光譜影像采集過程中不同時間因子對應(yīng)的大氣輻射傳輸參數(shù)，輔助分析時間變化的影響。

2.2 相對輻射歸一化校正方法

對于多時相衛(wèi)星影像的相對輻射歸一化校正，主要有基于全局影像統(tǒng)計信息的輻射歸一化方法和基于樣本集的輻射歸一化方法[15]。

基于全局影像統(tǒng)計信息的輻射歸一化方法，如直方圖匹配、平均值-標準偏差歸一化法等，容易造成原始光譜特征的扭曲，不適用于高光譜影像處理[16]。

基于樣本集的輻射歸一化方法是通過選擇不變的像元作為控制集樣本進行輻射歸一化，如偽不變特征法(pseudo-invariant feature method, PIF)[17]、時相不變集群法(temporally invariant cluster, TIC)[18]、全景回歸法(wall-to-wall regression, WWR)[19]、基于加權(quán)不變點歸一化法(weighted invariant pixels，WIP)[20]等。以上方法都是基于選擇的樣本集構(gòu)建線性回歸方程，實現(xiàn)輻射歸一化計算，主要區(qū)別在于選擇樣本集的方法不同。常用的偽不變特征法是通過選取多景遙感影像中的特征地物點，以參考影像的特征點光譜值為自變量，其他影像的特征點光譜值因變量，通過線性回歸求解校正系數(shù)[21]。

主要采用基于樣本集的歸一化的思想，在航空高光譜影像的旁向重疊區(qū)域內(nèi)選取同名地物像點作為樣本集，具體步驟如下。

(1)選取兩條航帶S1、S2重疊區(qū)域中的同名地物像元作為樣本集，并分別記錄S1、S2航帶中同名地物像元的反射率數(shù)值、成像時間，即P1={V11，V12,…,V1n}，P2={V21，V22,…,V2n}。其中，n為同名地物像元編號，V1n、V2n分別表示第n個同名地物點在航帶1和航帶2中的光譜值。

(2)假定單航帶影像的采集時間長度為T，開始采集時間為Tstart，設(shè)置時間差閾值tthr(0

(3)

依據(jù)記錄的同名像元成像時間，將樣本集P1、P2分為k個樣本子集{P11,P12,…,P1k}，{P21,P22,…,P2k}，每個子集的成像時間計為tk,可表示為

(4)

(3)假定選取航帶S1為參考航帶(可選擇明暗地物所在的航帶)，則以P1k中的同名地物像元光譜值為參考值，P2k中的同名地物像元光譜值為因變量，通過線性回歸求解出線性系數(shù){(a1,b1)，(a2,b2), …, (ak,bk)}，其中a為斜率、b為截距。

(4)以每個子集的成像時間tk為自變量，線性系數(shù)為因變量，通過線性回歸計算出航帶S2的對應(yīng)的線性校正系數(shù){c,d}。

(5)以航帶S2的像元成像時間為自變量，結(jié)合線性系數(shù){c,d}，對航帶S2進行校正計算。

3 實驗結(jié)果分析

首先對大氣輻射校正參數(shù)進行對比，然后在航帶中選擇4種典型地物，對比其反射率光譜校正的結(jié)果，分析確認時間變化的影響。進一步對相對輻射歸一化方法的校正結(jié)果，與常用的偽不變特征點方法校正結(jié)果進行對比。

3.1 大氣輻射校正參數(shù)對比

實際數(shù)據(jù)采集過程中，在每條航帶中布設(shè)明暗地物是難以實施的，通常是在個別航帶中布設(shè)明暗地物，并以其校正系數(shù)對其他航帶進行大氣輻射校正?；诖髿廨椛鋫鬏斈Ｐ偷男Ｕ彩且院綆Р杉^程中的某一時間點為參考進行大氣校正參數(shù)解算。然而，航空高光譜遙感(尤其是推掃式傳感器)每條航帶的數(shù)據(jù)采集需要一定時間，且作業(yè)時間內(nèi)太陽輻射、大氣條件處于持續(xù)變化之中，而這些變化必定對傳感器接收的輻亮度信息和反射率反演產(chǎn)生一定的影響。

3.1.1 線性回歸系數(shù)對比

實驗中在航線1、航帶2的首端都布設(shè)了黑白布靶標，航空高光譜影像中記錄的黑白布輻亮度DN值如圖2所示。隨著時間變化，反射率強的白布對應(yīng)的輻亮度DN值變化較為明顯，可見光譜段的DN值差異在220～530[即0.22～0.53 μW/(cm2·sr·nm)]。

圖2 黑白布輻亮度DN值

分別基于首端、末端黑白布輻亮度值和實測的黑白布反射率進行經(jīng)驗線性回歸計算。表1列出了部分波段對應(yīng)的線性回歸系數(shù)。黑白布的反射率是一定的，而不同時間對應(yīng)的黑白布輻亮度存在差異，因此計算的校正系數(shù)也不同。

表1 黑白布1、黑白布2對應(yīng)的校正系數(shù)

3.1.2 大氣輻射傳輸參數(shù)對比

基于Modtran軟件模擬計算獲取了單航帶采集時間(15 min，13：06—13：21)，以及兩條航帶整體采集時間內(nèi)(37 min，13：06—13：43)對應(yīng)的大氣輻射傳輸分量，并進行對比。

總體上，大氣上行輻射、地面至傳感間的大氣透過率、大氣半球反照率變化非常小。大氣上行輻射在可見光譜段稍大，單航帶時間內(nèi)在1.04×10-3～8.88×10-3μW/(cm2·sr·nm)，兩航帶時間內(nèi)在1.43×10-3～11.78×10-3μW/(cm2·sr·nm)。大氣透過率的變化，只有水汽吸收波段(820、940 nm)在10-3～10-4量級，其他波段在10-5量級。大氣半球反照率十分穩(wěn)定，數(shù)值變化在10-5～10-6量級。

地面接收的下行總輻射存在一定變化，但不同波段的變化也不同，如圖3所示?？梢姽庾V段(400～760 nm)差值單航帶時間(15 min)內(nèi)在0.30～0.49 μW/(cm2·sr·nm)(即DN值為300～490)，兩航帶時間(37 min)內(nèi)為0.50～0.85 μW/(cm2·sr·nm)(即DN值為500～850)；940～960 nm波段受水汽吸收影響，變化較?。黄渌V段的單航帶時間(15 min)內(nèi)變化在0.15～0.30 μW/(cm2·sr·nm)(即DN值為150～300)，兩航帶時間(37 min)內(nèi)在0.27～0.51 μW/(cm2·sr·nm)(即DN值為270～510)。這一點與圖2中白布DN值的差異相對應(yīng)。

圖3 不同時間對應(yīng)的下行總輻射

3.2 地物反射率光譜對比分析

3.2.1 絕對輻射校正結(jié)果分析

選取實驗區(qū)5種典型地物的反射率反演結(jié)果進行對比，分析時間變化對大氣輻射傳輸模型、經(jīng)驗線性模型校正結(jié)果的影響。表2為航帶重疊區(qū)同名地物影像采集時間。如圖4～圖8所示，分別為黃土、沙土、石英正長巖、鹽堿的反射率光譜曲線對比。

表2 航帶重疊區(qū)同名地物影像采集時間

選取航帶重疊區(qū)域中不同分布位置的4種典型地物，對比其反射率校正結(jié)果。

圖4中所選的黃土位于航帶的北端。兩條航帶中的采集時間差為35 min。經(jīng)驗線性法校正結(jié)果中，航帶1中黃土的反射率計算結(jié)果與航帶2的結(jié)果數(shù)值差異較大，差值在0.02～0.96。航帶2中分別采用1號黑白布和2號黑白布的校正結(jié)果也存在一定差異，但整體一致性較好。

圖4 黃土反射率計算結(jié)果

圖5中所選的沙土地物，兩條航帶中的采集時間差為29 min。經(jīng)驗線性法校正結(jié)果中，航帶1中黃土的反射率計算結(jié)果與航帶2的反射率結(jié)果在可見光譜段差值較小，在0.01～0.02，近紅外譜段差值較大，在0.03～0.09。航帶2中分別采用1號黑白布和2號黑白布的校正結(jié)果也存在一定差異，但整體一致性較好。

圖5 沙土反射率計算結(jié)果

圖6中所選的石英正長巖，兩條航帶中的采集時間差為15 min。15 min時間內(nèi)，太陽輻照度等的變化并不大。經(jīng)驗線性法的校正結(jié)果主要是近紅外波段略有差異。此外，在382.5、396.9、411.3 nm 3個波段出現(xiàn)負值，且數(shù)值不一致。主要是由于高程變化引起的，1號航帶中黑白布鋪設(shè)地點高程為海拔980 m，2號航帶中黑白布鋪設(shè)地點高程為海拔759 m，石英正長巖地物點的高程為海拔904 m。航空高光譜數(shù)據(jù)采集過程中，飛行高度一定，地面高程變化造成傳感器對地距離的變化，而不同傳輸距離對應(yīng)的輻射傳輸分量在短波波段差異較大[22]。圖4、圖5中黃土、沙土地物的高程分別為968、1 012 m，與1號黑白布、測區(qū)高程中值接近，受高程變化影響較小。

圖6 石英正長巖反射率計算結(jié)果

圖7中測區(qū)南端的鹽堿地物，兩條航帶中的采集時間差為9 min，且自身輻亮度(反射率)值較高，校正結(jié)果的數(shù)值一致性較好。

圖7 鹽堿反射率計算結(jié)果

3.2.2 相對輻射校正結(jié)果對比

實驗中，設(shè)置時間差閾值為5 min，自數(shù)據(jù)采集時間開始每隔5 min時間內(nèi)的同名地物作為一個樣本子集，共構(gòu)建3個樣本子集，以航帶1為參考航帶，分別計算出每個樣本子集的回歸系數(shù)，進一步結(jié)合成像時間計算航帶2的回歸系數(shù)。

此外，基于選定的同名地物像元，采用偽不變特征法，同樣以航帶1為參考航帶，計算回歸系數(shù)進行校正，與本文方法進行對比，結(jié)果如圖8～圖11所示。

圖8 黃土反射率校正前后對比

圖9 沙土反射率校正前后對比

圖10 石英正長巖反射率校正前后對比

圖11 鹽堿反射率校正前后對比

以均方根誤差作為歸一化校正結(jié)果與參考反射率數(shù)值的相似性度量，均方根誤差的值越小，則表明校正前后的光譜數(shù)值越接近。具體統(tǒng)計結(jié)果如表3所示。通過對均方根誤差的統(tǒng)計對比，所提出的方法校正結(jié)果優(yōu)于常用的偽不變特征點法。

表3 均方根誤差統(tǒng)計

4 結(jié)論

分析了航空高光譜影像數(shù)據(jù)采集過程中，太陽輻射、大氣傳輸?shù)入S著時間變化，對地物反射率反演帶來的影響。得出如下主要結(jié)論。

(1)隨著時間變化，不同時刻地物接收的下行總輻射存在一定差異，在37 min時間內(nèi)可達0.27～0.51 μW/(cm2·sr·nm)，不可忽略。對于經(jīng)驗線性校正法而言，不同時刻的明暗參考地物所計算的校正系數(shù)不同，會對大氣輻射校正結(jié)果產(chǎn)生一定的影響。

(2)選取兩條航帶重疊區(qū)域4類同名地物的反射率光譜對比，基于經(jīng)驗線性法的反演結(jié)果，在反射率譜形方面基本一致，但是反射率數(shù)值存在一定偏差。同名地物采集的時間間隔越長(間隔30 min以上)，校正結(jié)果的數(shù)值偏差越大。

(3)提出的相對輻射歸一化方法以航帶間同名地物像元作為不變特征點，并加入了時間因子影響，能夠有效地去除航帶數(shù)據(jù)采集過程中太陽輻射、大氣變化等因素的影響，提高航帶間同名地物像元的光譜一致性。