趙世湖，尹 丹

（1.國(guó)家測(cè)繪地理信息局衛(wèi)星測(cè)繪應(yīng)用中心，北京 100830；2.北京大學(xué)物理學(xué)院，北京 100871）

光學(xué)遙感數(shù)字相機(jī)具有作業(yè)效率高、影像動(dòng)態(tài)范圍大、幾何精度穩(wěn)定等優(yōu)點(diǎn)，同時(shí)帶來了多中心投影、視場(chǎng)拼接等新問題［1，2］。傳統(tǒng)的光學(xué)傳遞函數(shù)（Optical Transfer Function，OTF）模型及理論將光學(xué)遙感成像系統(tǒng)等效為具有空間移不變特性的空間頻率濾波器，描述了光學(xué)成像系統(tǒng)的物與像空間頻譜關(guān)系［3，4］。在膠片式航空攝影中，航空膠片解像力極高，可達(dá)200線對(duì)／mm以上，可視為對(duì)連續(xù)地物空間信號(hào)的連續(xù)記錄。因此，基于空間移不變特性的光學(xué)傳遞函數(shù)在一定程度上能夠反映膠片航攝相機(jī)系統(tǒng)對(duì)遙感地物中各個(gè)空間頻率分量的響應(yīng)特性［5］。隨著數(shù)字成像傳感器（如Charge Coupled Device，CCD）的應(yīng)用，基于CCD數(shù)字傳感器的模擬／數(shù)字（Analog／Discrete，A／D）信號(hào)離散采樣的空間移變效應(yīng)凸顯，即輸入的模擬影像信號(hào)在時(shí)間和空間的移位將導(dǎo)致數(shù)字離散采樣輸出的變化［6，7］。數(shù)字遙感成像系統(tǒng)已不再具有空間移不變特性，而恰恰是空間移變系統(tǒng)［2］。因此，在空間移變條件下，基于傳統(tǒng)傳遞函數(shù)理論的光學(xué)系統(tǒng)點(diǎn)擴(kuò)散函數(shù)與物函數(shù)卷積研究數(shù)字成像過程，是一種不嚴(yán)格的近似方法［8］。本文建立了CCD數(shù)字傳感器空間移變離散采樣過程的數(shù)學(xué)函數(shù)表征，分析了數(shù)字傳感器離散采樣物理過程所引起物像頻譜的空間移變響應(yīng)，并通過數(shù)字離散采樣模擬仿真程序驗(yàn)證了空間移變特性及其影響。本研究將對(duì)光學(xué)數(shù)字成像系統(tǒng)優(yōu)化設(shè)計(jì)與質(zhì)量評(píng)價(jià)、遙感影像信息高精度提取具有參考意義。

1 傳遞函數(shù)在光學(xué)遙感數(shù)字成像中的局限性

調(diào)制傳遞函數(shù)（Modulation Transfer Function，MTF）的定義是影像的調(diào)制度與景物的調(diào)制度的比值，稱為模量傳遞函數(shù)值，其反映了成像系統(tǒng)的空間頻率響應(yīng)特性，是評(píng)價(jià)成像系統(tǒng)質(zhì)量的重要指標(biāo)。光學(xué)傳遞函數(shù)理論建立的重要假設(shè)前提是系統(tǒng)具有空間移變特性，而CCD數(shù)字成像傳感器是一個(gè)采樣離散化的過程，是空間移變系統(tǒng)。在空間移不變條件不成立的情況下，輸入的連續(xù)物面函數(shù)f（x，y）在時(shí)間和空間的移位將導(dǎo)致數(shù)字離散采樣輸出的變化，即F（x，y）不等于F′（x＋x′，y＋y′），如圖1所示。光學(xué)遙感數(shù)字成像系統(tǒng)凸顯空間移變特性方面很多，如動(dòng)態(tài)成像過程導(dǎo)致曝光過程中的像移、基于多CCD外視場(chǎng)拼接的航攝相機(jī)系統(tǒng)空間采樣不一致性以及航攝相機(jī)的內(nèi)參數(shù)變化等［9］。

圖1 數(shù)字傳感器采樣成像系統(tǒng)空間移變特性Fig.1 Digital sensor′s space variant character

2 光學(xué)遙感數(shù)字成像空間離散采樣數(shù)學(xué)分析

在數(shù)字遙感成像中，物面函數(shù)o（x，y，z）經(jīng)過光學(xué)鏡頭投影后在二維像面得到同樣連續(xù)的像面函數(shù)g（x，y），g（x，y）經(jīng)離散采樣后變?yōu)椴蓸雍瘮?shù)gs（x，y），如圖2所示。光學(xué)鏡頭系統(tǒng)PSF具有空間移不變特性，可表示為h（x，y），CCD傳感器具有空間移變特性，兩者綜合組成光學(xué)數(shù)字遙感成像系統(tǒng)的空間移變點(diǎn)擴(kuò)散函數(shù)，可表示為h（x，y，Δx，Δy）。h（x，y，Δx，Δy）不僅與空間點(diǎn)位置（x，y）有關(guān)，而且與圖像像元的亞像素位置點(diǎn)（Δx，Δy）有關(guān)。

圖2 數(shù)字成像空間信息轉(zhuǎn)換過程Fig.2 Digital sensor′s information transform

設(shè)任意輸入物面函數(shù)的采樣網(wǎng)格為（X，Y），其中X、Y都是整數(shù)，微小偏移量（Δx，Δy），且小于單位采樣網(wǎng)格大小，即0＜Δx＜X，0＜Δy＜Y，則在物面二維空間中的任一點(diǎn)光源函數(shù)可表示為：

點(diǎn)光源經(jīng)過光學(xué)鏡頭h（x，y）轉(zhuǎn)換后，得到在二維連續(xù)像面空間中的連續(xù)像面函數(shù)g（x，y；Δx，Δy）可表示為：

連續(xù)像面函數(shù)g（x，y；Δx，Δy）經(jīng)過CCD傳感器空間離散采樣后，得到抽樣函數(shù)gs（x，y；Δx，Δy）。利用二維δ函數(shù)comb（x／X，y／Y）進(jìn)行連續(xù)像面函數(shù)的離散采樣，即二維δ函數(shù)與像面函數(shù)的乘積，離散采樣函數(shù)gs（x，y；Δx，Δy）就可以表示為：

根據(jù)δ脈沖函數(shù)的篩選與線性相移性質(zhì)，二維離散采樣函數(shù)可變換為：

根據(jù)卷積定理的傅里葉變換性質(zhì)，二維采樣成像像面函數(shù)gs（x，y）的頻譜Gs（fx，fv）可表示為：

由此可見，在考慮空間移變特性的情況下，空間域上對(duì)連續(xù)像面函數(shù)g（x，y；Δx，Δy）的離散采樣，導(dǎo)致輸出函數(shù)gs（x，y；Δx，Δy）的頻譜Gs的兩個(gè)變化：1）周期性復(fù)現(xiàn)在頻率平面上的（n／X，m／Y）點(diǎn)為中心的位置；2）頻域中相位的線性移動(dòng)與縮放。

3 Monte Carlo光子追跡算法與遙感成像仿真

3.1 光學(xué)遙感成像的概率統(tǒng)計(jì)模型

基于概率統(tǒng)計(jì)的光學(xué)系統(tǒng)成像模型如圖3所示，從物面上Δx處發(fā)出的光子數(shù)量N與該面元上的光強(qiáng)φ（xi）Δxi呈正比，即：

式中：K為一定數(shù)量光子數(shù)與對(duì)應(yīng)光強(qiáng)的比例常數(shù)，φ（xi）為在Δxi處單位面積上的光強(qiáng)。從xi附近的Δx區(qū)域內(nèi)發(fā)出的光子概率為：

設(shè)∑φ（xi）＝1，則p（xi）＝φ（xi）。同理，在像面上，存在p（yi）＝φ（yi）。

圖3 成像過程的概率統(tǒng)計(jì)表達(dá)Fig.3 Probability statistics expression of imaging process

3.2 Monte Carlo光子追軌成像仿真算法

Monte Carlo（MC）光子追軌遙感成像算法是基于上述光學(xué)系統(tǒng)成像概率統(tǒng)計(jì)模型，通過抽樣和統(tǒng)計(jì)求解像面光強(qiáng)分布未知特征量。算法流程如圖4所示。

（1）光子隨機(jī)產(chǎn)生。產(chǎn)生初始光子的位置和方向是MC模擬的第一步。物面光子產(chǎn)生具有朗伯性，向半球空間隨機(jī)產(chǎn)生光子。考慮到提高模擬仿真效率，本研究采用產(chǎn)生光子包（可視為若干光子的集合體）的方法改進(jìn)MC方法的運(yùn)行效率，同時(shí)賦予光子包一定的權(quán)重值ω。

（2）移動(dòng)光子包。光子包產(chǎn)生之后，就要傳播微小的步長(zhǎng)距離△s，該距離采用不固定傳輸步長(zhǎng)方法計(jì)算。不固定的步長(zhǎng)值△s為：

圖4 MC光子追軌成像仿真算法流程Fig.4 Imaging simulation scheme of Monte Carlo Photon Tracing

式中：ξ1為在［0，1］之間隨機(jī)分布的數(shù)值。

（3）光子傳播過程中的位置與方向。光子傳播的位置與方向由5個(gè)參量（3個(gè)位置坐標(biāo)，2個(gè)方向角）構(gòu)成，分別采用三維空間笛卡兒坐標(biāo)系的3個(gè)位置坐標(biāo)和3個(gè)方向余弦表示，即：

式中：x、y和z代表光子當(dāng)前傳播起始位置，x′、y′和z′代表光子傳播步長(zhǎng)△s后的終止位置，ux、uy和uz代表當(dāng)前光子傳播方向與X、Y和Z軸夾角的余弦。

（4）光子吸收。只有進(jìn)入光學(xué)成像鏡頭視場(chǎng)的光子包經(jīng)過鏡頭散射、吸收后才能夠達(dá)到像面。光子包在鏡頭中和傳輸介質(zhì)中可能發(fā)生多次碰撞，一部分能量被吸收，一部分能量被散射，其散射和吸收比可通過單次散射反照率表示，即：

式中：μs為粒子的散射系數(shù)，μa為粒子的吸收系數(shù)，a描述了光子碰撞后被散射的概率。由此，碰撞后光子的剩余權(quán)重值ω′可表示為ω′＝a×ω。

（5）光子消亡。按照光子碰撞散射、吸收理論計(jì)算，光子剩余權(quán)重值ω′永不可能為0。但光子總是存在消亡的過程，因此，經(jīng)過多次散射、吸收作用后，當(dāng)光子剩余權(quán)重值小于設(shè)定的最小閾值時(shí)，則判定光子消亡。

（6）光子散射方向。光子碰撞后的運(yùn)動(dòng)方向由歸一化的相位函數(shù)P描述，計(jì)算光子散射后新的傳播方向。設(shè)該相位函數(shù)P定義為Henyey－Greenstein相位函數(shù)，即：

式中：g為各向異性因子，ξ2為［0，1］之間隨機(jī)數(shù)。光子每次碰撞后，新的傳播方向可表示為u′x、u′y和u′z。

式中：（θ，φ）為碰撞角度。

4 空間移變點(diǎn)擴(kuò)散函數(shù)模擬仿真結(jié)果

4.1 空間移變遙感成像仿真過程

首先，設(shè)CCD數(shù)字成像傳感器在物面空間的采樣網(wǎng)格（X，Y）的間隔為0.06mm，模擬生成黑白刃邊景物e（x，y），其中，模擬景物發(fā)生亞像元空間移變偏移，偏移量分別為0.00mm（無(wú)像素偏移）、0.02 mm（1／3個(gè)亞像元偏移）、0.04mm（2／3個(gè)亞像元偏移）和0.06mm（1個(gè)像元偏移），分別對(duì)應(yīng)表1中的C1、C2、C3和C4項(xiàng)；其次，基于前述光學(xué)遙感成像概率統(tǒng)計(jì)模型與Monte Carlo光子追跡方法，實(shí)現(xiàn)光學(xué)遙感數(shù)字采樣成像模擬，獲得黑白刃邊景物的光學(xué)遙感退化影像ex（x，y）；最后，通過分析黑白刃邊數(shù)字影像ex（x，y），提取移變點(diǎn)擴(kuò)散函數(shù)h（x，y）進(jìn)行定量對(duì)比分析。

在不同空間移變量條件下光學(xué)遙感數(shù)字成像黑白刃邊模擬結(jié)果如圖5所示。

圖5 空間移變CCD靶面光強(qiáng)分布仿真Fig.5 Space variant imaging photon distribution on CCD

通過不同空間移變條件下的黑白刃邊遙感影像函數(shù)ex（x）的一階差分提取不同亞像元空間移變偏移量下的系統(tǒng)線擴(kuò)散響應(yīng)，并采用一維高斯函數(shù)擬合得到系統(tǒng)的空間移變線擴(kuò)散函數(shù)lx（x），即：

由空間移變線擴(kuò)散函數(shù)lx（x），計(jì)算移變點(diǎn)擴(kuò)散函數(shù)hx（x，y）。lx（x）是系統(tǒng)的空間移變點(diǎn)擴(kuò)散函數(shù)hx（x，y）在線光源方向的積分，即：

對(duì)于大部分光學(xué)系統(tǒng)而言，點(diǎn)擴(kuò)散函數(shù)h（x，y）趨于Gauss型，即：

其中：a為高斯函數(shù)的幅值系數(shù)；b為高斯函數(shù)的對(duì)稱軸位置因子；c為高斯函數(shù)的形狀特征因子，它反映了點(diǎn)擴(kuò)散函數(shù)的寬度，c越小說明光學(xué)系統(tǒng)成像質(zhì)量越好；C是h（x，y）的圓形區(qū)域。

平均點(diǎn)擴(kuò)散函數(shù)ha（x，y）對(duì)不同空間移變量黑白刃邊數(shù)字遙感影像進(jìn)行加權(quán)平均后，利用點(diǎn)擴(kuò)散函數(shù)提取方法求得。根據(jù)上述點(diǎn)擴(kuò)散函數(shù)計(jì)算方法，不同空間移變量條件下的空間移變點(diǎn)擴(kuò)散函數(shù)hx（x，y）如圖6所示。

圖6 空間移變點(diǎn)擴(kuò)散函數(shù)高斯擬合結(jié)果Fig.6 Gauss fitting result of space variant PSF

4.2 空間移變點(diǎn)擴(kuò)散函數(shù)比較分析

空間移變點(diǎn)擴(kuò)散函數(shù)hx（x，y）的定量對(duì)比分析，主要從兩方面進(jìn)行考察：

（1）ha（x，y）與hx（x，y）的形狀特征因子c誤差Error＿c，主要分析空間移變效應(yīng)對(duì)光學(xué)遙感成像系統(tǒng)成像質(zhì)量評(píng)價(jià)的影響。

式中：cx為空間點(diǎn)擴(kuò)散函數(shù)的形狀因子，ca為空間平均點(diǎn)擴(kuò)散函數(shù)的形狀因子。

（2）由ha（x，y）與hx（x，y）的對(duì)稱軸位置因子b可反算得到空間移變的估計(jì)值，并與空間移變量的理論值相比較，得到誤差Error＿b，主要分析空間移變效應(yīng)對(duì)點(diǎn)、線特征信息定位精度的影響。

式中：bx為空間移變點(diǎn)擴(kuò)散函數(shù)的對(duì)稱軸位置因子；b1為空間移變量為0.00mm的對(duì)稱軸位置因子；P為采樣網(wǎng)格的尺寸，本仿真中P＝0.06mm；m為與空間移變點(diǎn)擴(kuò)散函數(shù)對(duì)應(yīng)的理論空間移變量，例如0.02mm、0.04mm和0.06mm。

由表1可知，在0.00mm、0.02mm、0.04mm和0.06mm空間偏移量情況下，平均點(diǎn)擴(kuò)散函數(shù)ha（x，y）與移變點(diǎn)擴(kuò)散函數(shù)hx（x，y）的形狀特征因子誤差Error＿c分別為－15.55%、22.82%、－15.40%和－16.62%；位置因子誤差Error＿b分別為0%、5.3%、6.5%和0.9%。

表1 光學(xué)成像空間移變效應(yīng)仿真結(jié)果Table 1 Space variant simulation results

5 結(jié)論

本文深入分析了光學(xué)遙感數(shù)字成像在空間離散采樣條件下的空間移變特性，并通過模擬仿真系統(tǒng)進(jìn)行驗(yàn)證，對(duì)嚴(yán)密遙感成像模型建立與遙感信息提取具有參考意義。

其一，明確指出航空膠片影像與數(shù)字成像傳感器的根本區(qū)別之一——是否具有空間移變特性，是經(jīng)典光學(xué)傳遞函數(shù)理論應(yīng)用于遙感數(shù)字成像系統(tǒng)的局限性所在。數(shù)字成像傳感器對(duì)連續(xù)像面光強(qiáng)的離散采樣涉及空間信息的采樣與重建一系列遙感數(shù)據(jù)處理過程。

其二，根據(jù)數(shù)字成像的物理過程，引入物面函數(shù)的亞像元空間位移，建立了面陣CCD模／數(shù)轉(zhuǎn)換的數(shù)學(xué)方程，進(jìn)行了移變特性的數(shù)學(xué)分析，指出了連續(xù)像面函數(shù)g（x，y；Δx，Δy）經(jīng)離散采樣后函數(shù)gs（x，y；Δx，Δy）的頻譜Gs變化。

其三，在經(jīng)典MTF測(cè)量基礎(chǔ)上，利用能夠反映數(shù)字采樣空間移變特性的移變點(diǎn)擴(kuò)散傳遞函數(shù)提取方法，通過光子追跡遙感成像仿真算法，獲得在不同空間移變量條件下的數(shù)字影像，定量分析并驗(yàn)證了數(shù)字成像空間移變特性對(duì)系統(tǒng)點(diǎn)擴(kuò)散函數(shù)提取精度的影響。

［1］ 劉先林，段福州，宮輝力.航空攝影科技發(fā)展成就與未來展望［J］.前沿科學(xué)，2007（3）：10－14.

［2］ 董龍，李濤.CCD成像模型及仿真系統(tǒng)研究［J］.系統(tǒng)仿真學(xué)報(bào)，2009，21（4）：4520－4525.

［3］ 宋敏，鄶新凱，鄭亞茹.輸入信號(hào)與CCD像元位相匹配對(duì)調(diào)制傳遞函數(shù)的影響［J］.光學(xué)技術(shù)，1999，137（3）：84－86.

［4］ LEGER D D.MTF measurement using spotlight［J］.International Geoscience and Remote Sensing Symposium.IEEE Transactions on Geoscience and Remote Sensing，1994，4：2010－2012.

［5］ LEE H－C.Review of image－blur models in a photographic system using the principles of optics［J］.Optical Engineering，1990，29（5）：405－421.

［6］ 徐伯慶，李鑫欣.成像系統(tǒng)的空變模型及其空域逆濾波［J］.儀器儀表學(xué)報(bào)，2006，27（6）：2162－2168.

［7］ COSTELLO T P，MIKHAEL W B.Efficient restoration of space variant blurs from physical optics by sectioning with modified wiener filtering［J］.Digital Signal Processing，2003，13：1－22.

［8］ SITTER D N.Method for the measurement of the MTF of sampled imaging systems from bar－target patterns［J］.Applied Optics，1995，34（4）：746－751.

［9］ YAN L，ZHAO S H，LIAN Z H.A convertible aerial camera system integrated with film and CCD image sensor［C］.International Symposium on Photoelectronic Detection and Imaging，ISPDI2007，Proceeding of SPIE.2007.p662414.