交叉輻射對高分辨率遙感成像質(zhì)量的影響

齊懷川 王合順

（1 中國資源衛(wèi)星應(yīng)用中心，北京 100094）

（2 中國科學(xué)院遙感應(yīng)用研究所，北京 100101）

1 引言

介于遙感衛(wèi)星、地表目標(biāo)之間的大氣環(huán)境是遙感成像鏈路中的一個(gè)環(huán)節(jié)，是影響遙感衛(wèi)星數(shù)據(jù)質(zhì)量的不可忽略的因素[1]。進(jìn)入遙感器的光輻射除了來自目標(biāo)的輻射外，還包括一部分未到達(dá)目標(biāo)被大氣直接散射而進(jìn)入視場的光輻射，此部分稱之為程輻射；另一部分由于大氣散射導(dǎo)致目標(biāo)周圍背景進(jìn)入遙感器視場的一部分臨近像元輻射，此部分稱之為交叉輻射[2]。通常情況下，程輻射是影響遙感衛(wèi)星質(zhì)量、遙感定量反演的主要因素[3]，是大氣校正需要重點(diǎn)解決的問題，交叉輻射對質(zhì)量的影響可以忽略。但是隨著遙感衛(wèi)星空間分辨率的不斷提高，交叉輻射的影響逐漸增大，能夠計(jì)算不同大氣條件高分辨率觀測情況下的交叉輻射量是完善大氣模型的基礎(chǔ)，可以用于大氣效應(yīng)處理算法的開發(fā)改善高分辨率遙感衛(wèi)星的數(shù)據(jù)質(zhì)量[4-5]。

2 交叉輻射計(jì)算方法

在定義場景目標(biāo)背景表面為朗伯面的情況下，入瞳處觀測地面表觀輻射亮度分布 Ls(x,y,λ)可表示為[6-7]

式中 h (x, y,θv)為交叉輻射造成的點(diǎn)擴(kuò)散函數(shù)；*為卷積計(jì)算；為照射到表面的輻照度；T (θv)為觀測方向余弦；ρ(x,y)為目標(biāo)場景反射率分布；Lp()λ為程輻射貢獻(xiàn)量。本文將利用蒙特卡洛隨機(jī)統(tǒng)計(jì)方法，統(tǒng)計(jì)計(jì)算得到 h (x,y,θv)，用并以此表征交叉輻射對遙感圖像質(zhì)量的影響。根據(jù)光子傳輸?shù)幕ヒ仔?，本?jié)通過蒙特卡洛算法模擬光子從遙感器和目標(biāo)之間的傳輸，統(tǒng)計(jì)光子的空間分布最終得到近似點(diǎn)擴(kuò)散函數(shù)。

2.1 蒙特卡洛算法

作為一種隨機(jī)抽樣的試驗(yàn)統(tǒng)計(jì)方法，蒙特卡洛算法的基本思路是建立一個(gè)待解決問題的隨機(jī)過程模型，使過程所求參數(shù)是該問題的解或是描述該問題解的重要參量[8]。最終通過若干次抽樣試驗(yàn)統(tǒng)計(jì)近似求出該問題的解。求解的精確程度同隨機(jī)過程描述的完整性和抽樣次數(shù)有關(guān)。結(jié)合本文的蒙特卡洛算法可歸結(jié)為4個(gè)步驟[9]：

1）構(gòu)造描述問題的隨機(jī)過程，光子傳遞作為一個(gè)隨機(jī)過程，主要問題是描述該過程當(dāng)中的如是否碰撞、碰撞粒子性質(zhì)、碰撞取向等環(huán)節(jié)。

2）實(shí)現(xiàn)從確定概率分布當(dāng)中抽樣，本文將光子傳輸主要過程抽樣定義為[0,1]之間的均勻分布，特殊分布主要是在光子發(fā)生散射后其方向的獲取事件當(dāng)中，這些事件的概率分布同散射相位函數(shù)有關(guān)。

3）實(shí)現(xiàn)過程的模擬后，建立過程隨機(jī)量作為求解問題的估計(jì)量。若這問題的數(shù)學(xué)期望正好是所求問題的解則這種估計(jì)量就是問題的無偏估計(jì)。

4）統(tǒng)計(jì)獲得數(shù)據(jù)，得到模擬事件的最終結(jié)果或結(jié)論。

2.2 利用蒙特卡洛算法描述光子傳輸問題

根據(jù)上述提到的蒙特卡洛算法將光子在大氣中的傳輸描述為若干隨機(jī)事件的組合，基本過程描述如下：光子從目標(biāo)（中心像元）反射直接進(jìn)入接收視場；光子從目標(biāo)反射經(jīng)一次反射進(jìn)入接收視場；來自目標(biāo)外的輻射經(jīng)散射進(jìn)入接收視場?；诿商乜逅惴ǖ墓庾觽鬏斢幸韵聨c(diǎn)假設(shè)[10]：

1）光子散射是獨(dú)立的；2）光子的偏振性不考慮；3）光子傳輸過程是互易的；4）不考慮湍流的影響；5）忽略光子與地氣間的多次反彈。

光子運(yùn)動(dòng)狀態(tài)描述如圖1所示。下面對光子移動(dòng)過程進(jìn)行詳細(xì)闡述，蒙特卡洛算法光子傳輸模型的建立首先需要設(shè)定光子的狀態(tài),這里引入狀態(tài)參量S，其主要包括光子的位置P（x,y,z）,運(yùn)動(dòng)狀態(tài)M（u,v），其中u為散射角；v為散射平面內(nèi)的方位角；l為光子走過的大氣光學(xué)厚度，本文稱之為光學(xué)距離，可表示為

式中 ka與ks分別為吸收和散射系數(shù)；W為光子能量狀態(tài)初始值，可設(shè)為1。設(shè)S（P，M，l，W）為光子的狀態(tài)信息，光子從發(fā)射到最后吸收或到達(dá)目標(biāo)表面可由一系列狀態(tài)參量（S1，S2，S3，…）組成利用矩陣，可表示為

光子的移動(dòng)過程按下面描述為：

1）按照S0[P0(0,0,z0),M0(0,0),l0,1]的初始狀態(tài)發(fā)射光子。

2）獲得 Si+1的狀態(tài)信息。按照蒙特卡洛算法首先定義l'=ln(1 -r)，r為[0,1]之間的均勻分布，由此得到 l '。如果 l '＜li，則表示第i+1 狀態(tài)發(fā)生碰撞，于是有l(wèi)i+1=li-l'得到此次位置處的光學(xué)距離，根據(jù)式（2）得到zi+1，并由xi+1=(zi－zi+1)tan(ui)cos(vi)、yi+1=(zi－zi+1)tan(ui)sin(vi),得到完整的Pi+1(x,y,z)，能量閾值控制可根據(jù)Wi+1=WiW 得到，其中W為散射發(fā)生的反照度。

3）對于散射角u，通過判斷k＜[k分子/（k分子+k氣溶膠）]的關(guān)系確定發(fā)生碰撞的光子是與大氣分子還是與大氣氣溶膠粒子發(fā)生碰撞。其中k分子為大氣分子散射系數(shù)；k氣溶膠為氣溶膠粒子散射系數(shù)；k 服從[0,1]之間的均勻分布。

4）重復(fù)步驟2）～3），直至光子到達(dá)地面；如果Wi+1小于設(shè)定的閾值，記錄并重新發(fā)射光子。

光子傳輸過程中形象描述如圖2所示，圖中描述了光子某次發(fā)生碰撞后的情形，程序流程圖如圖3所示。

圖1 光子運(yùn)動(dòng)狀態(tài)描述Fig.1 The description of photon movement

3 計(jì)算結(jié)果與分析

將光子傳輸?shù)碾S機(jī)過程通過編程實(shí)現(xiàn)，利用計(jì)算機(jī)循環(huán)控制開展試驗(yàn)并統(tǒng)計(jì)試驗(yàn)數(shù)據(jù)，點(diǎn)擴(kuò)散函數(shù)矩陣h(x,y)如圖3所示，出于計(jì)算量的考慮選擇5×5 單元表示，h(3,3)為中心像元，計(jì)算不同大氣、不同分辨率、不同觀測角度情況下交叉輻射導(dǎo)致的點(diǎn)擴(kuò)散函數(shù)即得中心像元h(3,3)和周圍各像元輻射的貢獻(xiàn)量。

圖2 模擬流程圖Fig.2 Simulation flow chart

3.1 交叉輻射貢獻(xiàn)量同大氣能見度的關(guān)系

能見度同交叉輻射關(guān)系參數(shù)設(shè)定如下：大氣能見度(vis)分別為25、15、5km；傳感器高度(h)為785 000m；空間分辨率(GSD)為2m；觀測天頂角(vza)為0°；波長(λ)為0.5μm。由這些參數(shù)結(jié)合點(diǎn)擴(kuò)散函數(shù)，通過交叉輻射計(jì)算可得的h(x,y)見表1～3。

大氣能見度降低，表示氣溶膠厚度增加，大氣散射必然加強(qiáng)。由表4 可知，隨著能見度的降低，中心像元的貢獻(xiàn)量逐漸減小，交叉輻射的影響逐漸增大，這和定性分析的結(jié)論一致。在GSD=2m 時(shí)，當(dāng)能見度為vis=5km，中心像元貢獻(xiàn)量為87.6%，根據(jù)此時(shí)的點(diǎn)擴(kuò)散函數(shù)h(x,y)獲得的在奈奎斯特頻率處（2m）MTF 的值約為0.862 4，較vis=30km 時(shí)下降了0.087 5。

圖3 點(diǎn)擴(kuò)散函數(shù)矩陣Fig.3 Point spread function matrix

表1 vis=25km 時(shí)h(x,y)的分布情況Tab.1 Distribution of h(x,y)when vis=25km

表2 vis=15km 時(shí)h(x,y)的分布情況Tab.2 Distribution of h(x,y)when vis=15km

表3 vis=5km 時(shí)h(x,y)的分布情況Tab.3 Distribution of h(x,y) when vis=5km

表4 不同能見度下中心像元的共享量Tab.4 Distribution of central pixel with various visibility

3.2 交叉輻射貢獻(xiàn)量同空間分辨率的關(guān)系

空間分辨率同交叉輻射關(guān)系參數(shù)設(shè)定如下：vis為15km；h為785 000m；GSD 分別為100、30、2m；觀測天頂角vza為0°；波長λ為0.5μm，可得的h(x,y)見表5～7。

表5 GSD=100m 時(shí)h(x,y)的分布情況Tab.5 Distribution of h(x,y) when GSD=100m

表6 GSD=30m 時(shí)h(x,y)的分布情況Tab.6 Distribution of h(x,y) when GSD=30m

表7 GSD=2m 時(shí)h(x,y)的分布情況Tab.7 Distribution of h(x,y)when GSD=2m

不同空間分辨率情況下，如表8所示，由于散射引起的交叉輻射逐漸凸顯，隨著分辨率的增加，中心像元貢獻(xiàn)量逐漸減小，MTF 值逐漸下降，但趨勢較緩，僅下降約0.02。當(dāng)GSD=2m，在vis=15km 時(shí)MTF為0.932 4 較GSD=100m 時(shí)下降了0.019 5。

表8 不同分辨率下中心像元共享量Tab.8 Distribution of central pixel with various GSDS

3.3 交叉輻射貢獻(xiàn)量同觀測天頂角度的關(guān)系

觀測天頂角同交叉輻射關(guān)系參數(shù)設(shè)定如下：大氣能見度vis為15km；傳感器高度h為785 000m；空間分辨率GSD為2m；觀測天頂角vza 分別為0°、30°、60°；波長λ為0.5μm，可得的h(x,y)如表9～11所示。

表9 vza=0°時(shí)h(x,y)的分布情況Tab.9 Distribution of h(x,y)when vza=0°

表10 vza=30°時(shí)h(x,y)的分布情況Tab.10 Distribution of h(x,y)when vza=30°

表11 vza=60°時(shí)h(x,y)的分布情況Tab.11 Distribution of h(x,y)when vza=60°

隨著觀測角度的增加，散射引起的交叉輻射逐漸增大，如表12所示，當(dāng)vza=60°時(shí)，點(diǎn)擴(kuò)散函數(shù)h(x,y)中心像元共享量相對0°時(shí)下降了22.3%，MTF 值為0.665，較0°減小了0.267 4，降幅明顯?？梢婋S著觀測角度的增加，交叉輻射效應(yīng)導(dǎo)致傳函降低的情況不能忽視。

表12 不同觀測角度下交叉輻射影響Tab.12 Distribution of central pixel with various vza

4 結(jié)束語

本文介紹了大氣對遙感衛(wèi)星圖像數(shù)據(jù)得影響，給出了一種基蒙特卡洛算法計(jì)算交叉輻射貢獻(xiàn)量的方法，通過交叉輻射貢獻(xiàn)量表征其影響下的點(diǎn)擴(kuò)散函數(shù)，分析得出交叉輻射同大氣能見度、空間分辨率、觀測角度的關(guān)系，大氣能見度和觀測角度對成像質(zhì)量的影響不能忽視，隨著空間分辨率的升高，交叉輻射影響增大，但較上述兩種因素來說較小。因此，在大氣質(zhì)量下降的今天，當(dāng)我國高分辨率大角度觀測時(shí)，大氣交叉輻射對成像性能的影響不能忽視，大氣校正對提高成像輻射質(zhì)量十分必要。

References)

[1]齊懷川.高分辨率可見光空間光學(xué)遙感成像鏈路分析研究[D].北京：中國空間技術(shù)研究院,2010.QI Huaichuan.Study on the Image Chain for High Resolution Visible Light Space Optical Remote Sensing[D].Beijing：China Academy of Space Technology,2010.(in Chinese)

[2]王小飛, 毛志華, 陳建裕.海岸帶區(qū)域SPOT 衛(wèi)星數(shù)據(jù)大氣校正處理[J].海洋學(xué)研究,2011,29(1)：68-72.WANG Xiaofei,MAO Zhihua,CHEN Jianyu.Atmospheric Correction of the SPOT Satellite Data of the Coastal Zones[J].Journal of Marine Sciences,2011,29(1)：68-72.(in Chinese)

[3]何紅艷,楊居奎,齊文雯.大氣對遙感衛(wèi)星圖像品質(zhì)的影響分析[J].航天返回與遙感,2011,32(2)：42-47.HE Hongyan,YANG Jukui,QI Wenwen.Analysis of Atmosphere’s Influence on Image Quality of Remote Sensing Satellite[J].Spacecraft Recovery & Remote Sensing,2011,32(2)：42-47.(in Chinese)

[4]劉兆軍, 周峰, 阮寧娟, 等.一種光學(xué)遙感成像系統(tǒng)優(yōu)化設(shè)計(jì)新方法研究[J].航天返回與遙感,2011,32(2)：34-41.LIU Zhaojun,ZHOU Feng,RUAN Ningjuan,et al.Study on Trade-off Design Method for Space Optical Remote Sensing Imaging System[J].Spacecraft Recovery & Remote Sensing,2011,32(2)：34-41.(in Chinese)

[5]陳世平.關(guān)于航天遙感的若干問題[J].航天返回與遙感,2011,32(3)：1-8.CHEN Shiping.Some Issues about Space Remote Sensing[J].Spacecraft Recovery & Remote Sensing,2011,32(3)：1-8.(in Chinese)

[6]Steve A Cota,Jabin T Bell.PICASSO-An End-to-End Image Simulation Tool for Space and Airborne Imaging Systems[J].Journal of Applied Remote Sensing,2010,4(1)：ArticleID 043535.

[7]徐希孺.遙感物理[M].北京：北京大學(xué)出版社,2003.XU Xiru.Physics for Remote Sensing[M].Beijing：Peking University Press,2003.(in Chinese)

[8]徐希孺, 王平榮.用蒙特-卡洛方法計(jì)算大氣點(diǎn)擴(kuò)散函數(shù)[J].遙感學(xué)報(bào),1999,3(4)：268-278.XU Xiru,WANG Pingrong.Computing Atmospheric Point Spread Function by Monte-Carlo Method[J].Journal of Remote Sensing,1999,3(4)：268-278.(in Chinese)

[9]Robert A.Remote Sensing-Models and Methods for Image Processing[M].Arizona America：Academic Press,2007.

[10]李宏順,劉偉.用逆向蒙特卡羅法分析衛(wèi)星遙感中的鄰近效應(yīng)[J].華中科技大學(xué)學(xué)報(bào),2004,32(11)：31-33.LI Hongshun,LIU Wei.Analysis of the Adjacency Effect in Satellite Remote Sensing by Using Backward Monte Carlo Method[J].J .Huazhong Univ.of Sci.& Tech,2004,32(11)：31-33.(in Chinese)