CN103675794A - 基于时空统一特性的航天光学遥感器成像仿真方法 - Google Patents

Abstract

基于时空统一特性的航天光学遥感器成像仿真方法，步骤为：（1）根据卫星轨道参数，计算得到在任意时刻GMT卫星位置与时间的关系以及卫星在GMT时刻星下点的经纬度；（2）根据卫星在GMT时刻星下点的经纬度，成像时刻地面目标反射率特性ρ以及太阳高度角θ，计算得到遥感器的入瞳辐亮度L(λ)；（3）结合遥感器参数和L(λ)，得到遥感器的信号S_target，将S_target依次经过线性放大、滤波和量化处理后，再叠加遥感器光学系统、遥感器探测器、遥感器电路和卫星平台的仿真MTF后，得到初始仿真图像；（4）对初始仿真图像进行压缩解压缩、辐射校正和MTF补偿后，得到用户所需的仿真图像。本发明方法可进行有效的光学遥感全链路成像仿真，大大提高光学遥感在轨成像质量。

Description

基于时空统一特性的航天光学遥感器成像仿真方法

技术领域

本发明属于航天光学遥感技术领域，涉及一种在地面仿真航天光学遥感器成像的方法。

背景技术

光学遥感图像质量是评价遥感数据好坏和遥感器成像能力的一个关键参数。随着高分辨率卫星的迫切需求，光学遥感图像质量对用户越来越重要，图像质量成为空间光学相机研制的重要评价指标之一。遥感器参数的科学合理设计成为遥感器研制初期的核心任务，是保证最终在轨获取高质量图像的关键。

当前光学遥感器信噪比的设计是依据在一定太阳高度角和一定地面反射率条件下的信噪比值来确定遥感器的硬件指标，包括遥感器口径、积分级数，使之满足总体对遥感器信噪比的要求。这种方法导致了实际成像过程中，当太阳高度角与地面反射率没有达到设计值时，遥感器在轨获取的图像信噪比将低于指标要求，从而造成图像偏暗，清晰度不够。因此迫切需要一种综合考虑时间和空间条件的光学遥感器成像仿真途径，可以根据卫星轨道实时在轨调整遥感器参数（如积分级数和增益），同时结合地面处理（相对辐射校正和MTF复原）来保证最佳的图像质量。

2007年第06期《大气与环境光学学报》上，中科院安徽光学精密机械研究所易维宁等人发表的《光学遥感图像仿真软件ORSIS介绍》一文公开介绍了光学遥感图像仿真软件ORSIS，该软件只有三个模块组成，即地面景观光学特性仿真模块、大气辐射传输仿真模块和传感器仿真模块。2008年第14期《系统仿真学报》上，中国科学院安徽光学精密机械研究所张冬英等人发表的《基于航空图像的航天光学遥感器成像的仿真》一文建立了从航空遥感图像到航天遥感图像的仿真模型，但该方法无法实现不同经纬度区域时间和空间一致变化情况下的在轨图像仿真。

上述两篇文章都只建立了从目标到遥感器的数学仿真模型，完成了信号仿真，但并未涉及整个成像链路的信号、噪声和MTF仿真。对于目前在轨成像质量问题及其优化，急需要解决两个问题：一是进行光学遥感全链路仿真技术研究，二是研究时间和空间统一情况下遥感器在轨图像质量优化的参数调研策略。目前尚未发现与此相关的文献报道。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供了一种结合统一时间和空间的航天光学遥感器成像仿真方法，可以在地面进行光学遥感器全链路仿真，从而准确的确定光学遥感器的成像参数。

本发明的技术解决方案是：基于时空统一特性的航天光学遥感器成像仿真方法，步骤如下：

（1）根据卫星轨道参数，确定卫星的地心距r和偏近点角E，

$r=\frac{a(1-e^2)}{1+e\cos f},\cos E=\frac{e-\cos f}{1+e\cos f}$

其中a为卫星轨道半长轴，e为卫星轨道偏心率，f为卫星轨道真近点角；

（2）将近地点时刻t_p作为卫星在轨道上运动时间的起算点，计算得到在任意时刻GMT卫星位置与时间的关系，

$\mathrm{GMT}-t_p=\sqrt{\frac{a^3}{\mathrm{\μ}}}(E-e\sin E)$

其中μ为地心引力常数；

（3）确定卫星在GMT时刻星下点的经纬度，其中：

经度S_lon=sin^-1[sinisin(ω+f)]

纬度S_lat=Ω+tan^-1[cositan(ω+f)]-ω_e(GMT-t_p)

i为卫星轨道倾角，ω为卫星轨道近地点幅角，Ω为卫星轨道升交点赤经，ω_e为地球自转角速度；

（4）根据卫星在GMT时刻星下点的经纬度，计算得到GMT时刻的太阳高度角θ，

θ=arcsin[sin(S_lon)×sin(S_lat)+sin(S_lon)cos(S_lat)cosGMT]；

（5）根据卫星在GMT时刻星下点的经纬度，成像时刻地面目标反射率特性ρ，以及太阳高度角θ，采用大气模型传输方法计算得到遥感器的入瞳辐亮度L(λ)，

$L\left(\mathrm{\λ}\right)=\mathrm{\τ}\frac{\mathrm{\ρ}}{\mathrm{\π}}[\cos \mathrm{\θ}{E}_0\left(\mathrm{\λ}\right)+{L_{\mathrm{solar}}}^{\↓}\left(\mathrm{\λ}\right)]+{L_{\mathrm{solar}}}^{\↑}\left(\mathrm{\λ}\right)$

其中E₀(λ)是大气层顶的太阳辐照度，τ为大气透过率，L_solar ^↓(λ)为太阳下行辐射，L_solar ^↑(λ)为太阳上行辐射；

（6）结合遥感器参数和步骤（5）所获得的入瞳辐亮度L(λ)，得到遥感器的信号

$S_{t\arg \mathrm{et}}=\frac{A_{\det \mathrm{ector}}\mathrm{\π}(1-\mathrm{\ϵ})n_{\mathrm{tdi}}{t}_{\mathrm{int}}}{4F^2\mathrm{hc}}\underset{{\mathrm{\λ}}_{\min }}{\overset{{\mathrm{\λ}}_{\max }}{\∫}}\mathrm{\η}\left(\mathrm{\λ}\right)L\left(\mathrm{\λ}\right){\mathrm{\τ}}_{\mathrm{optics}}\left(\mathrm{\λ}\right)\mathrm{\λd\λ}$

其中A_detector为探测器面积，ε为光学孔径面积遮拦比，对于无遮拦的光学系统ε=0，F为系统f数，n_tdi为TDI级数，t_int为成像系统的像元积分时间，h为普朗克常数，c为光速，η(λ)为CCD量子效率，λ_max和λ_min为光谱响应范围的上限和下限，τ_optics(λ)为光学系统的光谱透过率；

将S_target依次经过线性放大、滤波和量化处理后，叠加遥感器光学系统、遥感器探测器、遥感器电路和卫星平台的仿真MTF后，再加入探测器和电路的噪声，最终得到初始仿真图像；

（7）对初始仿真图像进行压缩解压缩、辐射校正和MTF补偿处理后，得到用户所需的仿真图像。

所述步骤（7）中的图像压缩解压缩方法为JPEG2000算法或者KAKADU算法。所述的MTF补偿方法为Wiener滤波复原或者小波复原。

本发明与现有技术相比的优点在于：本发明方法基于卫星在轨运行时时间和空间统一的特性，以准确的空间地面信息、大气信息和成像时间来仿真光学遥感器成像入瞳能量，从而开展遥感器在轨成像仿真和地面处理仿真。根据辐射传输理论，地面目标特性和大气环境特性都具有重要的时间性和空间性，本发明方法从卫星轨道特性出发获得成像时刻的时间和空间，不仅确定遥感器与太阳、地物之间的关系，还保证了光学遥感器成像时刻的地面空间信息、大气特性的正确性，从而确保了遥感载荷入瞳能量仿真的准确性。在此基础上，结合遥感器在轨成像模型，仿真不同成像条件和不同地面目标特性的图像，为在研遥感器参数优化设计提供了有效手段，为在轨遥感成像问题提供了仿真方法，为提高光学遥感卫星的成像质量打下了良好的基础。

附图说明

图1为本发明方法的流程框图。

具体实施方式

如图1所示，为本发明方法的具体实施途径。

（1）卫星轨道参数确定：

首先根据用户需求，设定卫星轨道六根数：真近点角、近地点幅角、轨道倾角、升交点赤经、偏心率和轨道半长轴，确定卫星上天后的轨道。

在计算卫星轨道参数之前，首先要介绍一下所用到的两个坐标系：地心轨道坐标系和赤道惯性坐标系。

地心轨道坐标系是：原点O在地球中心，x₀轴从地心指向航天器，也就是沿地心距矢量r，y₀轴在轨道平面内，垂直于矢量r，指向前方，z₀轴垂直于轨道平面，满足右手坐标系法则。

赤道惯性坐标系是：原点O在地球中心，x轴沿地球赤道面和黄道面的交线，指向春分点，z轴指向北极，y轴在赤道平面上垂直于x轴。

在地心轨道坐标系下，

地心距r： $r=\frac{a(1-e^2)}{1+e\cos f}$

偏近点角E： $\cos E=\frac{e-\cos f}{1+e\cos f}$

其中a为卫星轨道的半长轴，f为真近点角，e为偏心率。

（2）卫星成像时间确定：

根据卫星在轨道上运动时间的起算点，即近地点时刻t_p，可得到在GMT时卫星位置与时间关系为：

$\mathrm{GMT}-t_p=\sqrt{\frac{a^3}{\mathrm{\μ}}}(E-e\sin E)$

其中μ为地心引力常数，等于398600.44km³/s²。

（3）卫星成像空间确定：

在地心轨道坐标系中，卫星的位置坐标为：

x₀=rcosf

y₀=rsinf

z₀=0

将卫星的位置坐标由地心轨道坐标系转换成赤道惯性坐标系，则卫星的位置在赤道惯性坐标系中的坐标为：

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=\frac{a(1-e^2)}{1+e\cos f}\left[\begin{array}{c}\cos \mathrm{\Ω}\cos (\mathrm{\ω}+f)-\sin \mathrm{\Ω}\sin (\mathrm{\ω}+f)\cos i\\ \sin \mathrm{\Ω}\cos (\mathrm{\ω}+f)+\cos \mathrm{\Ω}\sin (\mathrm{\ω}+f)\cos i\\ \sin (\mathrm{\ω}+f)\sin i\end{array}\right]$

其中ω为近地点幅角，i为轨道倾角，Ω为升交点赤经。

则卫星在GMT时刻星下点经纬度为：

经度：S_lon=sin^-1[sinisin(ω+f)]

纬度：S_lat=Ω+tan^-1[cositan(ω+f)]-ω_e(GMT-t_p)

其中ω_e为地球自转角速度，为7.292×10^-5/s。

（4）成像时刻太阳高度角确定：

根据卫星成像星下点经纬度（S_lon，S_lat）和当前卫星过顶时间（GMT），得到当前成像时刻的太阳高度角θ：

θ=arcsin[sin(S_lon)×sin(S_lat)+sin(S_lon)cos(S_lat)cosGMT]

（5）遥感器入瞳辐射亮度计算：

根据卫星过顶成像的星下点经纬度，利用MODIS全球反射率产品数据即可确定成像空间中地面目标反射率特性ρ，结合太阳高度角θ，由大气模型传输方法可得到遥感器的入瞳辐亮度L(λ)。

$L\left(\mathrm{\λ}\right)=\mathrm{\τ}\frac{\mathrm{\ρ}}{\mathrm{\π}}[\cos \mathrm{\θ}{E}_0\left(\mathrm{\λ}\right)+{L_{\mathrm{solar}}}^{\↓}\left(\mathrm{\λ}\right)]+{L_{\mathrm{solar}}}^{\↑}\left(\mathrm{\λ}\right)$

其中E₀(λ)是大气层顶的太阳辐照度，τ为大气透过率，L_solar ^↓(λ)为太阳下行漫反射，L_solar ^↑(λ)为太阳散射程辐射，这些参量可以通过大气辐射传输模型MODTRAN仿真获得。

（6）在轨遥感器成像仿真：

光学遥感器仿真模型考虑光学遥感器是一个非常复杂的系统，其主要组件包括成像光学系统、CCD探测器和电子学系统三个方面，仿真模型主要描述成光学遥感器的信号、噪声和MTF。

1）信号模型：

结合遥感器参数和步骤（5）所获得的入瞳辐亮度，考虑大气的影响，以电子数表示采样成像式光学遥感器平均目标信号S_target，可由下式表示：

$S_{t\arg \mathrm{et}}=\frac{A_{\det \mathrm{ector}}\mathrm{\π}(1-\mathrm{\ϵ})n_{\mathrm{tdi}}{t}_{\mathrm{int}}}{4F^2\mathrm{hc}}\underset{{\mathrm{\λ}}_{\min }}{\overset{{\mathrm{\λ}}_{\max }}{\∫}}\mathrm{\η}\left(\mathrm{\λ}\right)L\left(\mathrm{\λ}\right){\mathrm{\τ}}_{\mathrm{optics}}\left(\mathrm{\λ}\right)\mathrm{\λd\λ}$

其中A_detector为探测器面积，ε为光学孔径面积遮拦比，对于无遮拦的光学系统ε=0，F为系统f数，n_tdi为TDI级数，t_int为成像系统的像元积分时间，h为普朗克常数，c为光速，η(λ)为CCD量子效率，λ_max和λ_min为光谱响应范围的上限和下限，τ_optics(λ)为光学系统的光谱透过率。

成像电路将对到达探测器的信号进行放大、滤波和量化等。放大对信号的作用可以近似的看作是一个线性过程，用表达式表示为：

S_amplified=kS_target+a

其中k为放大率，a为偏置量。

然后采用交流耦合放大目标与背景间的信号差，去除大范围的非均匀背景，其中交流耦合电路可以用高通滤波器实现；最后进行AD变换，将模拟信号转变为具有确定大小的数字量。

2）噪声模型：

成像噪声模型中除了光子噪声外，还有遥感器自身的噪声，主要有探测器噪声和电路噪声。在对探测器噪声进行模型化处理时将主要考虑散粒噪声、复位噪声、片内放大器噪声、片外放大器噪声、量化噪声和模式噪声。电路在完成成像功能时，也会引入各种噪声。但主要可分为两类：片外放大噪声表征电路的等效信号放大器的噪声和量化噪声表征量化的误差。具体可参考《光学技术》期刊第38卷第3期《光子成像系统信号与噪声的统计特性分析》，《红外技术》第30卷第12期《8192像元TDICCD相机信噪比的深入分析》，和《红外与激光工程》第41卷第4期《科学级CCD三维噪声模型及其评估》等三篇文章中的噪声模型，并依据冈萨雷斯的《数字图像处理》（MTALBA版）104-113页的方法进行噪声的加入。

3）MTF模型：

根据中国宇航出版社出版的《空间相设计与试验》中介绍的光学遥感器成像过程，遥感器在轨成像退化包括了遥感器的系统噪声退化、遥感器MTF退化和在轨卫星平台MTF退化。其中MTF退化包括遥感器和平台的MTF退化，如下所述：

遥感器MFT模型包括光学系统（角标optics）、探测器（角标detector）和电路（角标cir）三种MTF模型，也即

MTF_sensor=MTF_optics×MTF_detector×MTF_cir

其中光学系统MTF模型为：

MTF_optics=MTF_diff×MTF_aberration×MTF_defocus×MTF_other

其中MTF_diff为衍射MTF，MTF_aberration为像差MTF，MTF_defocus为离焦MTF，MTF_other为其它MTF。

探测器MTF模型为：

MTF_detector=MTF_sl×MTF_diffusion×MTF_{trans_ccd}×MTF'_ccd

其中MTF_sl为探测器几何传递函数，MTF_diffusion为光电子扩散传递函数，MTF_{trans_ccd}为电荷转换效率传递函数，MTF'_ccd为探测器其它传递函数。

电路MTF模型为：

MTF_cir=MTF_amp×MTF_filter

其中MTF_amp为放大器调制传递函数，MTF_filter为滤波器调制传递函数。

卫星平台MTF模型：

MTF_plat=MTF_linear×MTF_sin×MTF_random×MTF_unnormal

其中MTF_linear为卫星平台线性运动MTF，MTF_sin为余弦运动MTF，MTF_random为随机运动MTF，MTF_unnormal为非正常运动MTF。

总的MTF退化可表示：

MTF_all=MTF_optics*MTF_detector*MTF_cir*MTF_plat

其中MTF_optics，MTF_detector，MTF_cir，MTF_plat分别为光学系统、探测器、电路和卫星平台的MTF。

参考冈萨雷斯的《数字图像处理》（MTALBA版），卫星遥感器MTF退化对图像的影响可理解为对图像的一种模糊滤波，因此可采用如下方法进行实际图像的退化仿真：

假定f(x,y)为原目标图像，g(x,y)为观察到的退化图像，则

g(x,y)=f(x,y)*h(x,y)+n(x,y)

其中，h(x,y)为点扩展函数，n(x,y)为噪声函数，符号“*”表示卷积。

上式表示原图像通过与相应的点扩展函数做卷积并加上噪声而得到退化图像的过程，在实际操作中，对上式进行傅里叶后再进行执行：

G(u,v)=F(u,v)H(u,v)+N(u,v)

其中H(u,v)的值为复数，可以表示为

其中|H|为幅值，

为相位。若对幅值作归一化，使得零频率的幅值为1，则称此归一化的幅值为调制传递函数（MTF），即

MTF=|H|/k

其中，k为H在零频率的幅值。

假定MTF在作用于频率图像F时，其以频谱中心为圆心的等半径圆上的MTF值相同，则可以得到

G=F·MTF·k+N

将遥感器和平台的仿真MTF作用于原始频率图像F，然后将得到的频率上的退化图像G进行傅里叶反变换即可得到最终的光学遥感仿真图像。

（7）经过遥感器和卫星平台的仿真后，最后进行压缩解压缩、辐射校正和MTF补偿的仿真，至此完成时空统一的光学遥感全链路成像仿真。

1）压缩解压缩

采用JPEG2000和KAKADU两种算法实现不同压缩比下的图像压缩和解压缩。

2）相对辐射校正

相对辐射校正的一般公式如下：

${\mathrm{DN}}_{\mathrm{cali}}=\frac{{\mathrm{DN}}_{\mathrm{rawi}}-B_i}{{\mathrm{NG}}_i}$

其中DN_cali为相对辐射校正后第i号探元的DN值，DN_rawi为原始图像第i号探元的DN值，B_i为第i号探元的偏移值，NG_i为第i号探元的归一化增益，增益值

在1.0左右，保证定标后的DN值范围与原始数值的范围基本一致。

3）MTF补偿

MTF补偿方法可采用两种方法：Wiener滤波复原和小波复原。

Wiener滤波复原是一种对噪声起抑制和减少作用的方法，由C.W.Helstrom提出。Wiener滤波复原是寻找一个滤波器，使得复原后的图像和原图像的均方差最小。

小波分析相对于傅立叶分析而言，具有多分辨分析的特点，在空域和频域都有表征图像局部信息的能力，空间窗和频率窗都可以根据图像的具体形态动态调整，在低频部分（图像较平稳）可以采用较低的空间分辨率，而提高频率的分辨率，在高频情况下（频率变化不大）可以用较低的频率分辨率来换取精确的空间定位。

（8）通过步骤（1）～（7）就可产生在轨光学遥感器成像结果经地面处理后的仿真图像。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims (3)

1.基于时空统一特性的航天光学遥感器成像仿真方法，其特征在于步骤如下：

（1）根据卫星轨道参数，确定卫星的地心距r和偏近点角E，

$r=\frac{a(1-e^2)}{1+e\cos f},\cos E=\frac{e-\cos f}{1+e\cos f}$

其中a为卫星轨道半长轴，e为卫星轨道偏心率，f为卫星轨道真近点角；

（2）将近地点时刻t_p作为卫星在轨道上运动时间的起算点，计算得到在任意时刻GMT卫星位置与时间的关系，

$\mathrm{GMT}-t_p=\sqrt{\frac{a^3}{\mathrm{\μ}}}(E-e\sin E)$

其中μ为地心引力常数；

（3）确定卫星在GMT时刻星下点的经纬度，其中：

经度S_lon=sin^-1[sinisin(ω+f)]

纬度S_lat=Ω+tan^-1[cositan(ω+f)]-ω_e(GMT-t_p)

i为卫星轨道倾角，ω为卫星轨道近地点幅角，Ω为卫星轨道升交点赤经，ω_e为地球自转角速度；

（4）根据卫星在GMT时刻星下点的经纬度，计算得到GMT时刻的太阳高度角θ，

θ=arcsin[sin(S_lon)×sin(S_lat)+sin(S_lon)cos(S_lat)cosGMT]；

（5）根据卫星在GMT时刻星下点的经纬度，成像时刻地面目标反射率特性ρ，以及太阳高度角θ，采用大气模型传输方法计算得到遥感器的入瞳辐亮度L(λ)，

$L\left(\mathrm{\λ}\right)=\mathrm{\τ}\frac{\mathrm{\ρ}}{\mathrm{\π}}[\cos \mathrm{\θ}{E}_0\left(\mathrm{\λ}\right)+{L_{\mathrm{solar}}}^{\↓}\left(\mathrm{\λ}\right)]+{L_{\mathrm{solar}}}^{\↑}\left(\mathrm{\λ}\right)$

其中E₀(λ)是大气层顶的太阳辐照度，τ为大气透过率，L_solar ^↓(λ)为太阳下行辐射，L_solar ^↑(λ)为太阳上行辐射；

（6）结合遥感器参数和步骤（5）所获得的入瞳辐亮度L(λ)，得到遥感器的信号

$S_{t\arg \mathrm{et}}=\frac{A_{\det \mathrm{ector}}\mathrm{\π}(1-\mathrm{\ϵ})n_{\mathrm{tdi}}{t}_{\mathrm{int}}}{4F^2\mathrm{hc}}\underset{{\mathrm{\λ}}_{\min }}{\overset{{\mathrm{\λ}}_{\max }}{\∫}}\mathrm{\η}\left(\mathrm{\λ}\right)L\left(\mathrm{\λ}\right){\mathrm{\τ}}_{\mathrm{optics}}\left(\mathrm{\λ}\right)\mathrm{\λd\λ}$

其中A_detector为探测器面积，ε为光学孔径面积遮拦比，对于无遮拦的光学系统ε=0，F为系统f数，n_tdi为TDI级数，t_int为成像系统的像元积分时间，h为普朗克常数，c为光速，η(λ)为CCD量子效率，λ_max和λ_min为光谱响应范围的上限和下限，τ_optics(λ)为光学系统的光谱透过率；

将S_target依次经过线性放大、滤波和量化处理后，叠加遥感器光学系统、遥感器探测器、遥感器电路和卫星平台的仿真MTF后，再加入探测器和电路的噪声，最终得到初始仿真图像；

（7）对初始仿真图像进行压缩解压缩、辐射校正和MTF补偿处理后，得到用户所需的仿真图像。

2.根据权利要求1所述的基于时空统一特性的航天光学遥感器成像仿真方法，其特征在于：所述步骤（7）中的图像压缩解压缩方法为JPEG2000算法或者KAKADU算法。

3.根据权利要求1所述的基于时空统一特性的航天光学遥感器成像仿真方法，其特征在于：所述的MTF补偿方法为Wiener滤波复原或者小波复原。