伊丕源， 李瀚波， 童鵬， 趙英俊， 張川， 田豐， 車永飛， 吳文歡

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

0 引言

在山區(qū)，由于地形的存在造成山區(qū)各部位地面目標(biāo)接收的太陽輻射能量有很大差異，也改變了太陽—地面目標(biāo)—傳感器之間的幾何位置關(guān)系，從而導(dǎo)致傳感器得到的測量值與地物真實的光譜反射率或輻射亮度之間存在差異，形成地形效應(yīng)。許多研究表明，地形效應(yīng)會影響甚至扭曲地表的反射特性，進(jìn)而會影響地表反照率的估算，并帶來較大的不確定性[1-3]。地形效應(yīng)是影響遙感定量化分析與應(yīng)用的主要障礙之一[4]，為提高遙感影像的利用價值，必須消除地形影響，即開展地形輻射校正。地形輻射校正研究是當(dāng)前遙感科學(xué)與技術(shù)發(fā)展的重要學(xué)科前沿之一[5]，國內(nèi)外學(xué)者對此開展了大量的研究，但目前的地形輻射校正方法大多是針對衛(wèi)星遙感影像，著重于消除地形角度的影響[6-13]。然而，地形因子主要包括高程因子和地形角度因子，對于衛(wèi)星遙感而言，由于衛(wèi)星軌道通常在對地距離100 km以上，所以地形高程起伏導(dǎo)致的輻射傳輸路徑變化相對較小，但對于航空高光譜遙感而言，由于其對地距離近，地形起伏造成傳感器與地物之間的輻射傳輸路徑變化較大，因此除地形角度外，地形高程變化的影響亦不可忽略。

本文主要針對地形高程變化的影響，在假定每一個高程點為水平朗伯體的前提下，計算分析了地形起伏導(dǎo)致的輻射傳輸距離變化對航空高光譜遙感成像造成的影響，并實現(xiàn)了加入高程因子的大氣輻射校正。

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

實驗數(shù)據(jù)為2012年9月12日在青海省雪鞍山地區(qū)，采用國王飛機(jī)同機(jī)搭載的CASI(compact airborne spectrographic imager)成像光譜儀和ALTM Gemini激光雷達(dá)設(shè)備獲取。設(shè)備性能參數(shù)如表1所示。

依據(jù)測區(qū)實際狀況與設(shè)備性能，設(shè)置航飛高度為6 900 m，采集的CASI航空高光譜數(shù)據(jù)波譜范圍： 380～1 050 nm(可見光—近紅外波段)，共36個波段，光譜分辨率為19 nm，預(yù)處理后插值生成的影像數(shù)據(jù)空間分辨率為0.75 m，高光譜影像像元值為輻射亮度，單位為10-3μWcm-2·sr-1·nm-1，即1 000DN=1 μWcm-2·sr-1·nm-1。航空激光雷達(dá)數(shù)據(jù)采集使用33 kHz脈沖頻率，視場角為±21°，平均點密度為0.3 個/m2，預(yù)處理后生成空間分辨率為5 m的數(shù)字高程模型(digital elevation model，DEM)。圖1(a)為實驗選擇的2012年9月12日下午15： 08采集的單航帶航空高光譜影像(采集時間600 s)，圖1(b)為其對應(yīng)的DEM，高程范圍為3 918～5 166 m，落差達(dá)1 248 m。

(b) DEM

2 實驗方法

鑒于本實驗重點在于分析高程因子變化對航空高光譜成像的影響，因此研究思路為： 假定每一個高程點為水平均一朗伯體，然后選擇相應(yīng)的高程范圍，分別解算出不同高程對應(yīng)的大氣輻射傳輸參量并對比分析。

2.1 基于MODTRAN輻射傳輸分量計算

在假定地表為水平均一朗伯體的情況下，傳感器接收的輻射亮度可表示為[14]

(1)

式中：Lpath為大氣上行輻射；S為大氣半球反照率；Fd為下行總輻射；T為地面像元與傳感器之間的大氣透過率；ρt為地表反射率。

采用MODTRAN軟件，選擇中緯度夏季(middle latitude summer)大氣模式、鄉(xiāng)村(rural，visible=23 km)氣溶膠模式； 設(shè)置傳感器垂直觀測。

選擇3 500 m，4 000 m，4 500 m，5 000 m，5 500 m共計5個地面高程值。以高程3 500 m為例，在MODTRAN4.0軟件中輸入設(shè)置的反射率和傳感器高度參數(shù)(如表2所示)，分別模擬計算傳感器在航飛高度(6 900 m)、地面高度(3 501 m)以及不同地表反射率數(shù)值(0，0.1，0.2)情況下接收的輻射亮度。然后將5組模擬計算的結(jié)果代入式(1)，解算出每個波段對應(yīng)的Lpath，S，F(xiàn)d和T。

表2 模擬輸入?yún)?shù)

參照以上參數(shù)設(shè)置，分別解算出其他高程值對應(yīng)的輻射傳輸分量。進(jìn)一步對5個高程值對應(yīng)的大氣輻射傳輸分量對比分析，以確定地形高程變化對傳感器接收信號的影響。

2.2 加入高程因子的大氣輻射校正

依據(jù)式(1)，ρt可表示為[15]

。

(2)

若經(jīng)2.1節(jié)結(jié)果分析，高程變化對大氣輻射傳輸分量有較大影響，則依據(jù)式(2)開展航空高光譜影像大氣輻射校正須考慮到各項參數(shù)受高程變化造成的影響。

實驗思路為依據(jù)一定數(shù)量高程值對應(yīng)的大氣輻射傳輸分量，通過函數(shù)擬合與插值計算，并結(jié)合高光譜像元對應(yīng)的高程值，計算出任意高程值對應(yīng)的各項大氣輻射傳輸參數(shù)，然后完成航空高光譜影像的反射率反演計算。

基于IDL語言，在ENVI+IDL編程模型下實現(xiàn)以上過程。具體實現(xiàn)步驟如下：

1)首先，將計算的各項大氣輻射傳輸分量依據(jù)標(biāo)準(zhǔn)高程值和波段存儲構(gòu)建為查找表，格式如圖2所示。將各項參數(shù)按照對應(yīng)的波段分開，逐波段將大氣參數(shù)生成矩陣，每一行對應(yīng)一個高程值，每一列對應(yīng)一項大氣參數(shù)。具體到本實驗中，包括3 500 m，4 000 m，4 500 m，5 000 m，5 500 m共5行，以及Lpath，T，S和Fd共4列。

圖2 查找表結(jié)構(gòu)示意圖

2)讀取DEM，與高光譜影像進(jìn)行坐標(biāo)范圍和空間分辨率匹配，并獲取每個高光譜像元對應(yīng)的高程值，生成高程數(shù)組。具體IDL代碼為

dem=envi_get_data(fid=fid_dem,dims=dims1,pos=0)

envi_enter_data,dem,r_fid=r_fid

resize_doit,fid=r_fid,dims=dims1,pos=0,/in_memory,interp=3,r_fid=r_fid1,rfact=rfact

elevation=envi_get_data(fid=r_fid1,dims=dims1,pos=0)

3)逐波段對各項大氣參數(shù)進(jìn)行回歸計算，解算出回歸系數(shù)，進(jìn)一步與高程數(shù)組運(yùn)算，解算出每個像元對應(yīng)的大氣輻射校正參數(shù)。具體IDL代碼為

lut=envi_get_data(fid=fid_lut,dims=dims_lut,pos=pos_lut[k])

HeightArr=[3500,4000,4500,5000,5500]

Lup=reform(lut[0,*,*],5,1)

Trans=reform(lut[1,*,*],5,1)

Sr=reform(lut[2,*,*],5,1)

Fd=reform(lut[3,*,*],5,1)

L0=regress(heightArr,Lup,const=L1)

T0=regress(heightArr,Trans,const=t1)

S0=regress(heightArr,Sr,const=s1)

F0=regress(heightArr,Fd,const=f1)

L[mask]=L1+elevation[mask]*L0

T[mask]=T1+elevation[mask]*T0

S[mask]=S1+elevation[mask]*S0

F[mask]=F1+elevation[mask]*F0

4)最后，逐波段、逐像元讀取高光譜影像輻射亮度數(shù)據(jù)，結(jié)合大氣參數(shù)，并基于式(2)進(jìn)行反射率計算。具體IDL代碼為

fork=0,num_bands-1do begin

band=envi_get_data(fid=fid_image,dims=dims_image,pos=pos_image[k])

image[mask,k]=(image[mask,k]-L[mask])/(F[mask]*T[mask]+S[mask]*(image[mask,k]-L[mask]))

endfor

image=reform(image,num_cols,num_rows,num_bands)

3 實驗結(jié)果與分析

3.1 大氣輻射傳輸分量分析

將3 500 m，4 000 m，4 500 m，5 000 m，5 500 m共5個高程值分別對應(yīng)的Lpath，T，S和Fd進(jìn)行對比分析。圖3為傳感器接收的5個不同高程的大氣上行輻射Lpath計算結(jié)果對比。

圖3 不同高程對應(yīng)的大氣上行輻射

圖3中可以看出，高程越低，與傳感器之間的距離越大，則Lpath值越大。尤其在可見光波段范圍內(nèi)，5個不同高程(對應(yīng)不同傳輸距離值)的Lpath值差異較大，這主要與該波段范圍內(nèi)大氣散射強(qiáng)烈有關(guān)。經(jīng)統(tǒng)計，可見光波段范圍內(nèi)，3 500 m與5 500 m對應(yīng)的Lpath值差異為0.23～1.23 μWcm-2·sr-1·nm-1，由于CASI影像的像元DN值單位為1 000DN=1 μWcm-2·sr-1·nm-1，則Lpath差異造成的信號DN值差異可達(dá)230～1 230之間，而且CASI高光譜影像大部分像元輻射亮度值(DN值)在2 000～8 000之間，因此這種差異影響較大。在波長大于800 nm的范圍內(nèi)，Lpath的差異逐漸變小，數(shù)值在0.05～0.14 μWcm-2·sr-1·nm-1之間，影響相對較小。由此可見，航空傳感器在獲取數(shù)據(jù)過程中，每個地面像元對應(yīng)的高程不同，其Lpath值均存在一定差異，且這種差異不可忽略。

圖4為傳感器至5個不同高程間的大氣透過率T計算結(jié)果對比。

圖4 不同高程至傳感器的大氣透過率

圖4中可見，不同高程對應(yīng)的T均不相同，高程越高，對傳感器距離值越小，則T越高，但除第32波段(938.4 nm)和第33波段(957.5 nm)外，其他波段的T差別在10-2數(shù)量級。而CASI數(shù)據(jù)的第32波段(938.4 nm)和第33波段(957.5 nm)處的T較低，主要是由于水汽在940 nm波長具有吸收作用。

圖5為不同高程對應(yīng)的大氣半球反照率S的結(jié)果對比。

圖5 不同高程對應(yīng)的大氣半球反照率

圖5中可見，首先，S的數(shù)值在短波波段范圍偏高，但數(shù)值也低于0.25，在近紅外波段范圍內(nèi)數(shù)值低于0.05。另外，高程越高，對應(yīng)的S數(shù)值越小。

圖6為不同高程接收的下行總輻射Fd的結(jié)果對比。

圖6 不同高程接收的下行總輻射

圖6中可見，高程越低，接收的Fd值越大。經(jīng)統(tǒng)計，3 500 m與5 500 m對應(yīng)的大部分波段下行輻射數(shù)值差異為1～2 μWcm-2·sr-1·nm-1，只有第32波段(938.4 nm)和第33波段(957.5 nm)，由于水汽影響大氣透過率差異較大，從而造成不同高程對應(yīng)的Fd差異更大，分別為3.1 μWcm-2·sr-1·nm-1和3.8 μWcm-2·sr-1·nm-1。因此，大部分波段下行輻射的DN值差異可達(dá)1 000～2 000，假定地物反射率為0.3，大氣透過率為0.95，則此部分對傳感器信號造成的DN值差異約在285～570范圍，其影響亦不可忽略。

綜合以上分析，不同高程對應(yīng)的Lpath，T，S和Fd均存在差異，且不可忽略。因此要實現(xiàn)精確的山地航空高光譜影像大氣輻射校正，必須考慮加入高程因子的影響。

3.2 反射率計算結(jié)果對比

本文所選實驗區(qū)為高海拔無人區(qū)，因此缺乏地面實測光譜數(shù)據(jù)為參考驗證。為分析加入高程因子校正后的效果，將解算的反射率結(jié)果與FLAASH大氣校正處理后的結(jié)果進(jìn)行對比分析。FLAASH是遙感圖像反射率反演的常用方法，基于MODTRAN輻射傳輸模型開發(fā)[16-17]。在FLAASH大氣校正過程中，同樣選擇中緯度夏季(middle latitude summer)大氣模式、鄉(xiāng)村(rural，visible=23 km)氣溶膠模式，高程參數(shù)選擇測區(qū)平均高程為4 500 m。

選擇實驗區(qū)中變質(zhì)巖地表、苔蘚植被、冰雪以及巖漿巖地表共4類地物(如圖7所示)，對反演后的反射率進(jìn)行對比。

(a) 所選特征地物分布

(b) 變質(zhì)巖地表(1號區(qū)域，高程4 920 m) (c) 苔蘚植物(2號區(qū)域，高程4 390 m)

(d) 冰雪(3號區(qū)域，高程4 950 m) (e) 巖漿巖地表(4號區(qū)域，高程4 260 m)

對比圖7中所選4類地物的FLAASH校正后與加入高程因子校正后的反射率，結(jié)果如圖8所示。

(a) 變質(zhì)巖FLAASH校正結(jié)果 (b) 變質(zhì)巖加入高程因子校正后結(jié)果

(c) 苔蘚FLAASH校正結(jié)果 (d) 苔蘚加入高程因子校正后結(jié)果

(e) 冰雪FLAASH校正結(jié)果 (f) 冰雪加入高程因子校正后結(jié)果

(g) 巖漿巖FLAASH校正結(jié)果 (h) 巖漿巖加入高程因子校正后結(jié)果

從圖8可以看出，F(xiàn)LAASH大氣校正結(jié)果與加入高程因子校正后的結(jié)果相比，其反射率譜形基本一致，主要差別在于短波波段，F(xiàn)LAASH校正計算的結(jié)果甚至出現(xiàn)負(fù)值，無疑是錯誤的。本文認(rèn)為如第3.1節(jié)所述，不同高程對應(yīng)的Lpath，T，S和Fd之間均存在差異，且波長越短差異越大，而FLAASH校正在假定地表為均一朗伯體的前提下，參考單一高程值進(jìn)行校正，從而導(dǎo)致反射率計算結(jié)果存在誤差，且在短波波段尤為明顯。因此，對于山區(qū)航空高光譜影像而言必須加入高程因子的影響，才能實現(xiàn)更為精確的大氣輻射校正。

4 結(jié)論

針對地形高程因子，在假定每一個高程點為水平朗伯體的前提下，計算分析了地形起伏導(dǎo)致的輻射傳輸距離變化對航空高光譜遙感成像的影響，并實現(xiàn)了加入高程因子的大氣輻射校正。主要結(jié)論如下：

1)基于MODTRAN的計算結(jié)果表明，地形高程的變化對航空高光譜影像的成像過程存在重要影響。

2)加入高程因子的大氣輻射校正結(jié)果與FLAASH大氣校正結(jié)果相比，兩者在反射率曲線的譜形方面總體一致，但反射率數(shù)值存在差異，尤其是在短波波段FLAASH大氣校正結(jié)果中甚至出現(xiàn)負(fù)值，無疑是錯誤的。

3)無論是基于MODTRAN計算的大氣輻射校正參數(shù)(大氣上行輻射、大氣透過率)，還是FLAASH大氣校正結(jié)果，都可以發(fā)現(xiàn)高程因子變化的影響在短波波段更為顯著，這主要與短波波段的大氣散射輻射強(qiáng)烈有關(guān)。

高程因子的變化對山地航空高光譜影像成像過程的影響不可忽略，要實現(xiàn)精確的大氣輻射校正就必須消除其影響。未來還需要從地形高程和地形角度2個方面共同開展研究，才能進(jìn)一步實現(xiàn)高精度的航空影像地形輻射校正。