光学遥感数字成像系统空间移变效应仿真与分析

赵世湖，尹 丹

（1.国家测绘地理信息局卫星测绘应用中心，北京 100830；2.北京大学物理学院，北京 100871）

光学遥感数字相机具有作业效率高、影像动态范围大、几何精度稳定等优点，同时带来了多中心投影、视场拼接等新问题［1，2］。传统的光学传递函数（Optical Transfer Function，OTF）模型及理论将光学遥感成像系统等效为具有空间移不变特性的空间频率滤波器，描述了光学成像系统的物与像空间频谱关系［3，4］。在胶片式航空摄影中，航空胶片解像力极高，可达200线对／mm以上，可视为对连续地物空间信号的连续记录。因此，基于空间移不变特性的光学传递函数在一定程度上能够反映胶片航摄相机系统对遥感地物中各个空间频率分量的响应特性［5］。随着数字成像传感器（如Charge Coupled Device，CCD）的应用，基于CCD数字传感器的模拟／数字（Analog／Discrete，A／D）信号离散采样的空间移变效应凸显，即输入的模拟影像信号在时间和空间的移位将导致数字离散采样输出的变化［6，7］。数字遥感成像系统已不再具有空间移不变特性，而恰恰是空间移变系统［2］。因此，在空间移变条件下，基于传统传递函数理论的光学系统点扩散函数与物函数卷积研究数字成像过程，是一种不严格的近似方法［8］。本文建立了CCD数字传感器空间移变离散采样过程的数学函数表征，分析了数字传感器离散采样物理过程所引起物像频谱的空间移变响应，并通过数字离散采样模拟仿真程序验证了空间移变特性及其影响。本研究将对光学数字成像系统优化设计与质量评价、遥感影像信息高精度提取具有参考意义。

1 传递函数在光学遥感数字成像中的局限性

调制传递函数（Modulation Transfer Function，MTF）的定义是影像的调制度与景物的调制度的比值，称为模量传递函数值，其反映了成像系统的空间频率响应特性，是评价成像系统质量的重要指标。光学传递函数理论建立的重要假设前提是系统具有空间移变特性，而CCD数字成像传感器是一个采样离散化的过程，是空间移变系统。在空间移不变条件不成立的情况下，输入的连续物面函数f（x，y）在时间和空间的移位将导致数字离散采样输出的变化，即F（x，y）不等于F′（x＋x′，y＋y′），如图1所示。光学遥感数字成像系统凸显空间移变特性方面很多，如动态成像过程导致曝光过程中的像移、基于多CCD外视场拼接的航摄相机系统空间采样不一致性以及航摄相机的内参数变化等［9］。

图1 数字传感器采样成像系统空间移变特性Fig.1 Digital sensor′s space variant character

2 光学遥感数字成像空间离散采样数学分析

在数字遥感成像中，物面函数o（x，y，z）经过光学镜头投影后在二维像面得到同样连续的像面函数g（x，y），g（x，y）经离散采样后变为采样函数gs（x，y），如图2所示。光学镜头系统PSF具有空间移不变特性，可表示为h（x，y），CCD传感器具有空间移变特性，两者综合组成光学数字遥感成像系统的空间移变点扩散函数，可表示为h（x，y，Δx，Δy）。h（x，y，Δx，Δy）不仅与空间点位置（x，y）有关，而且与图像像元的亚像素位置点（Δx，Δy）有关。

图2 数字成像空间信息转换过程Fig.2 Digital sensor′s information transform

设任意输入物面函数的采样网格为（X，Y），其中X、Y都是整数，微小偏移量（Δx，Δy），且小于单位采样网格大小，即0＜Δx＜X，0＜Δy＜Y，则在物面二维空间中的任一点光源函数可表示为：

点光源经过光学镜头h（x，y）转换后，得到在二维连续像面空间中的连续像面函数g（x，y；Δx，Δy）可表示为：

连续像面函数g（x，y；Δx，Δy）经过CCD传感器空间离散采样后，得到抽样函数gs（x，y；Δx，Δy）。利用二维δ函数comb（x／X，y／Y）进行连续像面函数的离散采样，即二维δ函数与像面函数的乘积，离散采样函数gs（x，y；Δx，Δy）就可以表示为：

根据δ脉冲函数的筛选与线性相移性质，二维离散采样函数可变换为：

根据卷积定理的傅里叶变换性质，二维采样成像像面函数gs（x，y）的频谱Gs（fx，fv）可表示为：

由此可见，在考虑空间移变特性的情况下，空间域上对连续像面函数g（x，y；Δx，Δy）的离散采样，导致输出函数gs（x，y；Δx，Δy）的频谱Gs的两个变化：1）周期性复现在频率平面上的（n／X，m／Y）点为中心的位置；2）频域中相位的线性移动与缩放。

3 Monte Carlo光子追迹算法与遥感成像仿真

3.1 光学遥感成像的概率统计模型

基于概率统计的光学系统成像模型如图3所示，从物面上Δx处发出的光子数量N与该面元上的光强φ（xi）Δxi呈正比，即：

式中：K为一定数量光子数与对应光强的比例常数，φ（xi）为在Δxi处单位面积上的光强。从xi附近的Δx区域内发出的光子概率为：

设∑φ（xi）＝1，则p（xi）＝φ（xi）。同理，在像面上，存在p（yi）＝φ（yi）。

图3 成像过程的概率统计表达Fig.3 Probability statistics expression of imaging process

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

Monte Carlo（MC）光子追轨遥感成像算法是基于上述光学系统成像概率统计模型，通过抽样和统计求解像面光强分布未知特征量。算法流程如图4所示。

（1）光子随机产生。产生初始光子的位置和方向是MC模拟的第一步。物面光子产生具有朗伯性，向半球空间随机产生光子。考虑到提高模拟仿真效率，本研究采用产生光子包（可视为若干光子的集合体）的方法改进MC方法的运行效率，同时赋予光子包一定的权重值ω。

（2）移动光子包。光子包产生之后，就要传播微小的步长距离△s，该距离采用不固定传输步长方法计算。不固定的步长值△s为：

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

式中：ξ1为在［0，1］之间随机分布的数值。

（3）光子传播过程中的位置与方向。光子传播的位置与方向由5个参量（3个位置坐标，2个方向角）构成，分别采用三维空间笛卡儿坐标系的3个位置坐标和3个方向余弦表示，即：

式中：x、y和z代表光子当前传播起始位置，x′、y′和z′代表光子传播步长△s后的终止位置，ux、uy和uz代表当前光子传播方向与X、Y和Z轴夹角的余弦。

（4）光子吸收。只有进入光学成像镜头视场的光子包经过镜头散射、吸收后才能够达到像面。光子包在镜头中和传输介质中可能发生多次碰撞，一部分能量被吸收，一部分能量被散射，其散射和吸收比可通过单次散射反照率表示，即：

式中：μs为粒子的散射系数，μa为粒子的吸收系数，a描述了光子碰撞后被散射的概率。由此，碰撞后光子的剩余权重值ω′可表示为ω′＝a×ω。

（5）光子消亡。按照光子碰撞散射、吸收理论计算，光子剩余权重值ω′永不可能为0。但光子总是存在消亡的过程，因此，经过多次散射、吸收作用后，当光子剩余权重值小于设定的最小阈值时，则判定光子消亡。

（6）光子散射方向。光子碰撞后的运动方向由归一化的相位函数P描述，计算光子散射后新的传播方向。设该相位函数P定义为Henyey－Greenstein相位函数，即：

式中：g为各向异性因子，ξ2为［0，1］之间随机数。光子每次碰撞后，新的传播方向可表示为u′x、u′y和u′z。

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

4 空间移变点扩散函数模拟仿真结果

4.1 空间移变遥感成像仿真过程

首先，设CCD数字成像传感器在物面空间的采样网格（X，Y）的间隔为0.06mm，模拟生成黑白刃边景物e（x，y），其中，模拟景物发生亚像元空间移变偏移，偏移量分别为0.00mm（无像素偏移）、0.02 mm（1／3个亚像元偏移）、0.04mm（2／3个亚像元偏移）和0.06mm（1个像元偏移），分别对应表1中的C1、C2、C3和C4项；其次，基于前述光学遥感成像概率统计模型与Monte Carlo光子追迹方法，实现光学遥感数字采样成像模拟，获得黑白刃边景物的光学遥感退化影像ex（x，y）；最后，通过分析黑白刃边数字影像ex（x，y），提取移变点扩散函数h（x，y）进行定量对比分析。

在不同空间移变量条件下光学遥感数字成像黑白刃边模拟结果如图5所示。

图5 空间移变CCD靶面光强分布仿真Fig.5 Space variant imaging photon distribution on CCD

通过不同空间移变条件下的黑白刃边遥感影像函数ex（x）的一阶差分提取不同亚像元空间移变偏移量下的系统线扩散响应，并采用一维高斯函数拟合得到系统的空间移变线扩散函数lx（x），即：

由空间移变线扩散函数lx（x），计算移变点扩散函数hx（x，y）。lx（x）是系统的空间移变点扩散函数hx（x，y）在线光源方向的积分，即：

对于大部分光学系统而言，点扩散函数h（x，y）趋于Gauss型，即：

其中：a为高斯函数的幅值系数；b为高斯函数的对称轴位置因子；c为高斯函数的形状特征因子，它反映了点扩散函数的宽度，c越小说明光学系统成像质量越好；C是h（x，y）的圆形区域。

平均点扩散函数ha（x，y）对不同空间移变量黑白刃边数字遥感影像进行加权平均后，利用点扩散函数提取方法求得。根据上述点扩散函数计算方法，不同空间移变量条件下的空间移变点扩散函数hx（x，y）如图6所示。

图6 空间移变点扩散函数高斯拟合结果Fig.6 Gauss fitting result of space variant PSF

4.2 空间移变点扩散函数比较分析

空间移变点扩散函数hx（x，y）的定量对比分析，主要从两方面进行考察：

（1）ha（x，y）与hx（x，y）的形状特征因子c误差Error＿c，主要分析空间移变效应对光学遥感成像系统成像质量评价的影响。

式中：cx为空间点扩散函数的形状因子，ca为空间平均点扩散函数的形状因子。

（2）由ha（x，y）与hx（x，y）的对称轴位置因子b可反算得到空间移变的估计值，并与空间移变量的理论值相比较，得到误差Error＿b，主要分析空间移变效应对点、线特征信息定位精度的影响。

式中：bx为空间移变点扩散函数的对称轴位置因子；b1为空间移变量为0.00mm的对称轴位置因子；P为采样网格的尺寸，本仿真中P＝0.06mm；m为与空间移变点扩散函数对应的理论空间移变量，例如0.02mm、0.04mm和0.06mm。

由表1可知，在0.00mm、0.02mm、0.04mm和0.06mm空间偏移量情况下，平均点扩散函数ha（x，y）与移变点扩散函数hx（x，y）的形状特征因子误差Error＿c分别为－15.55%、22.82%、－15.40%和－16.62%；位置因子误差Error＿b分别为0%、5.3%、6.5%和0.9%。

表1 光学成像空间移变效应仿真结果Table 1 Space variant simulation results

5 结论

本文深入分析了光学遥感数字成像在空间离散采样条件下的空间移变特性，并通过模拟仿真系统进行验证，对严密遥感成像模型建立与遥感信息提取具有参考意义。

其一，明确指出航空胶片影像与数字成像传感器的根本区别之一——是否具有空间移变特性，是经典光学传递函数理论应用于遥感数字成像系统的局限性所在。数字成像传感器对连续像面光强的离散采样涉及空间信息的采样与重建一系列遥感数据处理过程。

其二，根据数字成像的物理过程，引入物面函数的亚像元空间位移，建立了面阵CCD模／数转换的数学方程，进行了移变特性的数学分析，指出了连续像面函数g（x，y；Δx，Δy）经离散采样后函数gs（x，y；Δx，Δy）的频谱Gs变化。

其三，在经典MTF测量基础上，利用能够反映数字采样空间移变特性的移变点扩散传递函数提取方法，通过光子追迹遥感成像仿真算法，获得在不同空间移变量条件下的数字影像，定量分析并验证了数字成像空间移变特性对系统点扩散函数提取精度的影响。

［1］ 刘先林，段福州，宫辉力.航空摄影科技发展成就与未来展望［J］.前沿科学，2007（3）：10－14.

［2］ 董龙，李涛.CCD成像模型及仿真系统研究［J］.系统仿真学报，2009，21（4）：4520－4525.

［3］ 宋敏，郐新凯，郑亚茹.输入信号与CCD像元位相匹配对调制传递函数的影响［J］.光学技术，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］ 徐伯庆，李鑫欣.成像系统的空变模型及其空域逆滤波［J］.仪器仪表学报，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.