饶 鹏, 王成良
1.中国科学院上海技术物理研究所,上海200083
2.中国科学院红外探测与成像技术重点实验室,上海200083
3.空军预警学院空天预警实验室,武汉430019
随着航天技术的不断发展,气象卫星、资源卫星、海洋卫星以及侦察监视卫星等的广泛应用与发展极大地促进了对地观测光学载荷的发展,也使光学遥感技术得以迅猛发展,从定性逐渐走向定量化应用,形成了较多的定量数值产品.然而,影响定量化进一步提高的主要因素之一就是杂散光[1].杂散光是指光学系统中除了目标光线外,扩散于探测器表面上的其他非目标光线辐射能量,也就是到达光学系统像面的非成像光束.随着光学遥感定量化应用的深入,杂散光一直是这一领域的研究热点和难点[2-3],其中杂散光分析是研究基础,已成为该领域的热门研究方向.杂散光分析方法在20世纪70年代到80年代得到了长足发展[4],主要有光线追迹法、近轴近似法、光流密度追迹法、区域法和蒙特卡罗方法.蒙特卡罗方法[5]以适用范围最广、精度最高、建模便捷而成为当前应用最广泛的杂散光分析方法,其分析结果可用点源透过率(point source transmittance,PST)和杂光系数这两类指标来描述.点源透过率是按角度特性来表述系统的抗杂散光能力,对视场外光源幅度进行归一化,形状为理想点源;杂光系数则对视场内外光源幅度进行归一化,形状为面源.
当前的杂散光分析方法以及分析结果的表达主要存在两方面的不足:
1)杂散光二维分布表达不足.遥感图像上杂散光分布的差异是像素级的[6-7],于是要求像面上的杂散光分布达到像素级分辨率;而以往的杂散光分析方法则把光学像面当作一个整体处理,得到的分析结果是像面上平均的统计结果,不能分辨相邻像素单元的差别.
2)分析结果与遥感图像不能比对.遥感图像是在实际的运行环境中获得的,能真实地记录杂散光在像面上的分布,直观地反映杂散光的特征,是分析杂散光机理和来源的重要途径.然而,以往的杂散光分析方法缺少成像过程和光照环境等环节,其分析结果不能与遥感图像直接比对,造成了杂散光机理不清晰以及杂散光来源不明确等难题.
针对上述不足,为了具备杂散光仿真图像与遥感图像的可比性,本文提出一种基于成像过程仿真的杂散光分析方法,并通过试验验证了该方法的可行性和实用性.
基于成像过程仿真的杂散光分析主要依据有效载荷的扫描成像原理、焦平面探测器排列和杂散光二维分布进行空间采样和视场拼接,按不同杂散光机理形成图层化的杂散光图像.点源透过率PST函数是杂散光分析最常用的指标,是评价光学系统杂散光抑制能力的主要指标,它有两种定义形式:
1)归一化到轴上点源的点源透过率函数(point source transmittance normalized to the power from an on-axis point source,PSTN)
式中,Pd(φ)表示从离轴角φ的点源发出的辐射落在像面上的辐射通量,Pd(0)表示从位于轴上的同一点源发出的辐射落在像面的辐射通量.
2)点源垂直照度透过率函数(point sourcenormal irradiance transmittance,PSNIT)
式中,Ed(φ)表示由离轴角为φ的点源引起的像面辐照度,Ei表示垂直于该点源的系统入瞳处的辐照度.
从上面的两种定义中可以看出,PST为离轴角的函数,对于非圆对称光学系统,它还是方位角的函数.一般地,点源透过率随离轴角的增大衰减得很快,在点源透过率对离轴角的曲线中,纵坐标即点源透过率往往是按对数分度的.
点源透过率通过两种方法获得:一是直接测量再转化,二是通过软件分析计算.在软件分析计算中,需要知道光学系统的几何结构、镜面的杂散光属性和非镜面元件的杂散光属性.
本文主要根据软件分析计算的方法,建立空间遥感仪器成像通道的点源透过率曲线.采用软件仿真杂散光的关键步骤有两点:一是建立准确的光机系统结构,二是建立各个关键表面和重要表面的准确的光学散射特性.
若要准确地建立光机系统各个关键表面和重要表面的光学散射特性,需要引入一个重要的参数进行描述,即双向反射分布函数(bidirectional ref lectance distribution function,BRDF),它是散射辐亮度与入射辐照度的比值,其单位是球面度的倒数
式中,ES(φi,θi)为入射的辐照度,LS(φr,θr)为散射的辐亮度,φi、θi分别为入射辐射的方位角和俯仰角,φr、θr分别为出射辐射的方位角和俯仰角.
杂散光图像主要研究空间遥感仪器由于太阳光引起的杂散光,并生成杂散光的图像.通过外光源参数和点源透过率曲线,可以得到仪器各个遥感通道太阳光引起的杂散光分布.杂散光仿真流程如图1所示.
图1 基于成像过程的杂散光仿真Figure 1 Simulation of the stray light based on imaging process
利用杂散光仿真图像分析杂散光的技术特点如下:
1)在像平面分辨率上,直接按焦平面探测器的规模网格化,而且按像素在像面上的排列位置进行布置,网格化的单元可以直接模拟焦平面的像素单元.
2)在图像仿真形成方面,按统一的时间推进,可分为两种情况,一是光机扫描成像需要模拟动态扫描成像过程,由焦平面甚至单元视场镶嵌成扫描视场;二是卫星平台推扫成像需要模拟动态推扫成像过程,由焦平面(线列)条带视场镶嵌成覆盖视场.
3)在杂散光源仿真输入方面,按统一的时间推进,依照时刻确定杂散光源的方位和光照环境,模拟动态成像过程.
4)对杂散光的来源和机理进行分类处理,形成不同类别的图层(通道),杂散光图像为全部图层的叠加.
可以看出:物空间角度采样、光线追迹光线数设置以及成像过程动态建模是基于成像过程仿真的杂散光分析方法需要解决的主要问题.
1.2.1 物空间角度采样的优化
杂散光图像的分布带宽可以分为低频、中频和高频三部分.杂散光源、光机结构和成像过程都能造成不同带宽响应,如太阳、月亮、恒星等都接近点源,形成系统的高频输入.光路内的机械支撑、限制光束的光阑等,形成系统的高频响应.获得光学遥感图像之前并不知道杂散光图像带宽的分布,较直观的方法是直接按像元或亚像元级对应的角度空间(瞬时视场)划分二维角度空间网格,但会产生巨大的计算量.例如,地球静止轨道对地观测获得地球圆盘图用14µrad×14µrad瞬时视场的单元探测器采样18°×9°的覆盖视场,采样数达5.3×108个.按当前的主流计算机配置,一个采样点的计算时间约为1.5 h,由此推断出串行计算总共需要约9万年!即使使用并行计算技术,资源消耗也是惊人的.因此,有必要结合外光源和系统响应的空间频率特性确定合适的空间二维采样网格.
首先,对外光源进行傅里叶变换,计算频谱并设置阈值,获得有效的空间最高频率.一般而言,面源成像集中于低频区,因此面源成像的空间最高频率远低于高空间分辨率仪器的截至频率,其有效的空间最高频率对应的空间角度采样尺寸远大于仪器的瞬时视场,估计采样数可降低近万倍.
其次,按照对应的空间角度采样尺寸进行所有采样数的正向光线追迹,可以找到照明面.从像面逆向光线追迹,可以找到关键面.可以追迹到既是照明面又是关键面的面,一般称为重要面.通过这类方法可以获得所有采样数的重要面的公共面积,从公共面积的一阶、二阶微分获得重要面的起伏数据,以反映系统的响应频谱情况,然后设置阈值,并判断局部采样数的增减.
最后,可以根据不同的频率特性进行分类处理,形成不同的图层.通过以上方法划分物方二维角度空间网格,一般可以将二维采样数降低106倍,如果仍以地球静止轨道对地观测获得地球圆盘图为例,此时串行计算总共需要约800 h.
1.2.2 光线追迹数设置的优化
仪器背景辐射计算和杂散光仿真都需要大量的光线追迹计算,其中初始设置的光线从光源出发,经过系统后分散出更多的光线,其中只有少数到达焦平面探测器.而焦平面探测器上每一个像素收集能量后会形成二维照度分布.通常由于每个像素收集到的光线数不够,而杂散光又会在图像上表现出一定的粒度或背景噪声(可用信噪比衡量),因此会导致仿真不准.虽然追迹更多的光线和分散更多的光线可以降低粒度噪声,但该工作耗费时间较多.因此,如何确定合适的光线数是一个过程优化的问题.借鉴PST计算中不断增加光线数直到结果稳定的方法,杂散光图像仿真也采用这种反复迭代的方法,并对方法进一步改进,通过累加每一次的计算结果,减少重复的光线追迹.这样,可以有效解决光线追迹的信噪比问题.
这里涉及到初始光线数设置问题,通常直观的方法是直接按物空间最小特征尺寸在入瞳处设置光线数,这将产生很大的计算量.例如,光学口径φ600 mm的卡塞格林双反射光学系统,次(副)镜的机械支撑是近厘米量级的,需要亚毫米分布光线,光线数达3×107根.如果把次(副)镜的机械支撑放到另一个图层处理,光线数则可以降低4个数量级.因此,很有必要根据不同图层公共面积的起伏确定合适的初始光线数.按照成像关系,把重要面投影到光学系统的入瞳位置,以分类尺寸如分米级、厘米级、毫米级和亚毫米级分别设置光线数.
1.2.3 成像过程动态建模
外部光源如太阳、星星和地球的方位等都可能随时间变化,仪器视轴也会做扫描运动或指向运动.一些参与成像过程仪器的光机结构元件会随时间进行运动,同时外部光源的遮挡、重要性采样投向等也都在变化,即场景和成像都是动态变化的.因此,需要进行成像过程的动态建模.
几何模型、成像过程、外光源等动态建模问题可归结为参数化建模,且这些参数是随时间变化的,同时也要求保持像面与坐标系某个轴垂直,像面的边与另两个轴平行.为此,场景等采用反向运动技术.另外,重要性取样投射方向也需要随着时间动态改变,投射方向的角度余弦变化也需要参数化.
太空是用于开展杂散光试验验证的最佳实验室:无尘、无悬浮物和真空,除了仪器基本上没有其他散射光源;太阳是一个高功率、光谱已知、准直的光源;仪器本身就是测试对象,不需要中间环节,图像反映的是真实的杂散光.通过在轨遥感仪器进行方法验证是比较有效的方法.
风云二号多通道扫描辐射计作为核心载荷安装在地球静止轨道气象卫星上,对地球进行扫描成像.其R-C光学系统主光轴垂直于卫星自旋轴,平行于地球赤道平面,星下点指向赤道.利用卫星自旋和辐射计内的望远镜步进对地球进行两维扫描以获得地球云图,如图2所示.由于卫星自旋,辐射计实现对地球自西向东扫描,每当卫星自旋一周,辐射计望远镜就向前步进一步,实现由北向南扫描.扫描辐射计探测视场扫过地球后,望远镜自北向南步进一步,步距角140µrad,完成一条线的扫描.每幅地球圆盘图共扫2 500条线,自北向南扫过20°,可获得一幅20°×20°的地球全景图.
图2 FY-2扫描辐射计对地扫描工作原理Figure 2 Imaging principle of FY-2 meteorological satellite
作为地球静止轨道遥感卫星,风云二号气象卫星云图也存在杂散光的影响,特别是进入地球阴影区季节的午夜图像,与其他国家的气象卫星的云图一样,图像被太阳光污染严重,如2012年2月19日当地时间午夜时分的水汽通道(6.30~7.60µm)图像如图3所示.
午夜太阳光接近点源特性,为高频输入.如果系统有高频响应,图像上就有高频特征,于是可以利用这些特征与实际图像进行比对,检验方法的正确性和有效性.因此,本试验以高频特征仿真为主.
步骤1 建立扫描辐射计的光机模型[8],见图4.
图3 FY-2 E星扫描辐射计水汽通道(IR4)午夜云图Figure 3 IR4 nephogram of FY-2 E satellite
图4 FY-2扫描辐射计光机模型Figure 4 Model of optical mechanism of FY-2 VISSR
步骤2 进行图层划分.作为R-C光机系统,扫描辐射计的次镜支撑杆在光路中,太阳光直照支撑杆后出现散射,经主、次镜进入成像系统.由于次镜支撑杆是圆柱状且直径较小,散射光在反射方向很强.因为缺乏支撑杆表面的双向分布函数,这里只计算反射方向的杂散光,并只给出形状和位置信息.中等离轴角度的太阳光经过主、次镜会聚,照射到主镜室,散射后经过次镜进入系统,主镜室表面有螺纹且发黑处理,接近漫反射.小离轴角度的太阳光经过主、次镜会聚,照射到折镜光栏,散射后直接进入后光路.次镜遮光罩内壁为关键面,同时也是被照明面,中等离轴角度的太阳光经主镜会聚后,还会照亮次镜遮光罩内壁并发生散射,经过次镜或折镜进入系统.根据以上成像过程分析,分别为次镜支撑杆、主镜室、折镜光栏和次镜罩一次散射引起的杂散光建立图层.
步骤3 确定物空间角度采样.由于事先获得了云图,物空间角度采样已有先验知识.分析图3可知,图像高频特征的采样角度约为0.2°,以此作为物空间的采样角度.
步骤4 进行动态建模.坐标固连在焦平面上,辐射计东西方向的整体扫描可转化为外光源的逆向运动,每个时刻旋转一个角度采样,完成一条扫描线后,主镜、次镜和相应的结构件在南北方向旋转一个角度采样,折镜和相应的结构件也同方向旋转半个角度采样,以此完成光机运动的动态模型.重要性抽样的投射方向也需要相应的动态变化,如主镜和次镜散射后的投射方向需要在南北方向旋转一个角度采样.午夜太阳的位置可以从模拟的时间中算出,同时叠加一个东西方向的角度采样.
步骤5 进行蒙特卡罗光线追迹.
通过上述试验,得到了由次镜支撑杆、主镜室、折镜光栏和次镜罩等部件一次散射引起的杂散光图层图像,如图5所示.东西方向的次镜支撑杆影响非常明显,形状和位置非常接近于实际图像.南边的支撑杆在实际图像中有所发散,模拟的结果也有这个特征.
图5 次镜支撑杆图层图像(太阳照射到次镜支撑杆后发生散射,经过主、次镜进入系统)Figure 5 Simulation imageof knighthead for thesecond mirror
主镜室图层如图6所示,模拟的比实际的单薄,原因可能是其螺纹引起了关键面的扩散;模拟的西边亮东边暗,与实际的相反;模拟是按照设计的光路和结构件而不是实际装配的结果,存在一定的偏差,因此模拟的比实际的位置靠南;由于主镜室有一个通光孔,光线陷入孔中,成像形成了图像中的缺口,这是一个特征点,模拟图像清楚地反映了这一点.
图6 主镜室图层图像(太阳经过主、次镜会聚到主镜室后散射)Figure 6 Simulation image of the primary mirror cavity
折镜光栏图层如图7所示,模拟的比实际的单薄,原因是用平面结构来模拟实际有阶梯状的折镜光栏,与实际略有不同.
图7 折镜光栏图层图像(太阳经过主、次镜后会聚到折镜光栏并发生散射)Figure 7 Simulation image of diaphragm of the ref lector
次镜遮光罩内壁图层如图8所示,模拟与实际情况在形状和位置上存在一定的差异,这与光机模型精细化程度有关.
通过以上仿真图像与气象云图的对比,可以找到部分杂散光的来源,仿真生成的图像基本反映了实际云图的概貌,但精细化和定量化还需要开展进一步的研究.
基于成像过程的杂散光仿真图像实现了杂散光二维分布的精细表达,其结果可以与遥感图像直接比对.杂散光仿真图像具有以下优点:
图8 折镜光栏图层图像(太阳经过主、次镜后会聚到折镜光栏并发生散射)Figure 8 Simulation image of lens hood of the second mirror
1)有更为丰富的杂散光分布特征的表达力.传统的杂光系数和点源透过率都是抽象的,不能表达杂散光的分布特征,而杂散光图像达到像素级的高密度采样,与光学遥感图像具有同等的二维图像表达能力.
2)有更为高效的杂散光来源和机理评价能力.杂散光图像表达了系统在特定场景下的杂散光分布特征,“按图索骥”获得杂散光的来源和机理.通过杂散光的来源和机理可以为针对性的杂散光抑制设计提供重要的技术支持.
3)有按机理分类处理能力.杂散光图像按照机理进行图像分层,而光学遥感图像的杂散光是各种机理作用下的结果,因此可以为地面图像杂散光消除算法提供针对性的输入数据开展研究和论证.
[1]许健民,杨军,张志清,孙安来.我国气象卫星的发展与应用[J].气象,2010,36(7):94-100.
XU Jianmin,YANG Jun,ZHANG Zhiqin,SUN Anlai.Chinese meteorological satellites,achievements and applications[J].Meteorological Monthly,2010,36(7):94-100.(in Chinese)
[2]PUSCHELL J J,LOw E H A,JETER J W,et al.Japanese advanced meteorological imager[C]//Proceedings the International Society for Optical Engineering(SPIE),2005:5658.
[3]AMINOU D M,BÉZY J L,BENSI P.Meteosat third generation:preliminary imagery and sounding mission concepts and performances[C]//Proceedings the International Society for Optical Engineering(SPIE),2005,5978:K 1-K 14.
[4]BREAULT R P.Stray light technology overview of the 1980 decade(and a peek into the future)[C]//Proceedings the International Society for Optical Engineering(SPIE),1990,1331:2-11.
[5]LIKENESSB K.Stray 1ight simulation with advanced Monte Carlo techniques[C]//Proceedings the International Society for Optical Engineering(SPIE),1977,107:80-89.
[6]陈博洋,李欣耀,郭强,陈福春.FY-2E卫星杂散光评价与分析[J].红外技术,2010,32(11):636-639.
CHENBoyang,LIXinyao,GUOQiang,CHENFuchun.Estimate and analysis for stray light of FY-2E satellite[J].Infrared Technology,2010,32(11):636-639.(in Chinese)
[7]原育凯,李欣耀,贾伟.风云二号扫描辐射计可见杂散光的像元间差别[J].红外技术,2006,28(12):722-725.
YUAN Yukai,LI Xinyao,JIA Wei,Difference among pixels of FY-2 VISSR visible stray light[J].Infrared Technology,2006,28(12):722-725.(in Chinese)
[8]沈易,李欣耀,陈福春.风云二号辐射计的红外杂散光抑制研究[J].红外,2013,34(8):1-5.
SHENYi,LIXinyao,CHENFuchun.Research on suppression of stray infrared light in FY—2 radiometer[J].Infrared,2013,34(8):1-5.(in Chinese)