CN106652016B - 一种基于辐射几何一体化设计的光学遥感图像仿真方法 - Google Patents

Abstract

本发明一种基于辐射几何一体化设计的光学遥感图像仿真方法。基于辐射几何一体化设计的光学遥感图像仿真方法，根据仿真姿态数据、轨道数据、光学系统焦距及内方位元素、探测器尺寸、积分时间和级数建立探测器每个探元的成像几何模型，每个探元由成像几何模型在仿真场景高程数据上计算其物方坐标系三维坐标矩阵，在仿真场景反射率数据中获取三维矩阵中每个物方三维点的反射率，由计算太阳福照度和大气透过率计算得出仿真场景入瞳福亮度数据，入瞳福亮度数据经过光学仿真系统得到到达探测器的能量，探测器能量经过探测器仿真系统得到最终的仿真图像。本发明同时考虑几何成像和辐射成像关系，对空间遥感光学成像仿真程度高。

Description

一种基于辐射几何一体化设计的光学遥感图像仿真方法

技术领域

本发明属于图像仿真领域，涉及一种基于辐射几何一体化设计的光学遥感图像仿真方法。

背景技术

遥感图像仿真能够对遥感器在空间环境条件运行中遇到的各种状态进行模拟。通过仿真分析，设计人员可以在遥感器研制初期可靠地掌握有关信息，为提高遥感器性能，获得用户满意的图像而服务。

目前由于计算机仿真技术的快速发展，可以比较逼真地模仿各种目标与背景的各种辐射特性，模拟各种照明条件，对不同目标进行各种数学仿真。美国柯达公司设计了遥感器仿真软件，从1979年开始应用共经过了超过20000幅(次)成像验证。软件考虑了成像过程中的15个环节，主要关注辐射传输链路的模拟，成像几何关系语焉不详。北京理工大学开发了图像仿真系统，由实验室内可见光—近红外运动目标生成装置、动态特性测试系统、光量子噪声生成装置等组成，但不是针对真实地面三维信息进行遥感成像仿真。中科院长春光机所开发了图像仿真软件系统，将选定的输入图像目标经过处理，考察在假定的条件下，图像目标的质量下降(退化)程度，用于图像质量的预估。该软件也仅考虑辐射传输的情况。

随着卫星平台机动能力的增强，多角度全方位高时间分辨率的要求日益增加，成像几何关系对图像仿真的影响也越来越高，所以急需将遥感成像几何辐射联合仿真，才能更真实的反映成像条件，提高图像仿真精度。

发明内容

本发明所要解决的问题是：顾及成像几何对图像仿真的影响，提供一种基于辐射几何一体化设计的光学遥感图像仿真方法，提高光学遥感图像仿真精度。

本发明的技术方案是：一种基于辐射几何一体化设计的光学遥感图像仿真方法，包括以下步骤：

1)确定成像起始时间t₀，成像时间t，时标间隔Δt，采样点数n＝t/Δt，n取整为自然数，第p个采样点处对应的时间为t₀+p×Δt；以时间为顺序，在每个采样点处仿真得到包括成像时侧摆角、俯仰角和偏航角的姿态参数和成像时以摄影中心三维坐标表示的轨道参数(X_i,Y_i,Z_i)，其中i＝1,2…n，为采样点数；将每个姿态参数拟合为以时间为自变量，姿态参数为因变量的姿态仿真三次多项式；将每个轨道参数拟合为以时间为自变量，轨道参数为因变量的轨道仿真三次多项式；

2)确定光学系统焦距f，积分级数integral_level，积分时间integral_time,每一行探元数量为num_CCD，仿真图像行数num_line＝t/(integral_level×integral_time)，num_line取整数，第l行成像时刻为t_l＝t₀+l×(integral_level×integral_time)；l行对应的姿态数据和轨道数据(X_l,Y_l,Z_l)由t_l带入步骤1)中的姿态仿真三次多项式和轨道仿真三次多项式获得；

3)设定探元尺寸p，轨道高度H，地面采样间隔GSD＝pH/f，成像范围为num_line×num_CCD×GSD；

4)利用起始行轨道数据(X₁,Y₁,Z₁)，得到仿真场景范围为X方向X₁-num_CCD×GSD/2到X₁+num_CCD×GSD/2，Y方向为Y₁到Y₁+num_line×GSD；对已有的地面数字高程模型在仿真场景范围内进行分辨率为GSD的重采样，得到仿真场景高程数据；

5)对第l行建立任一探元的成像几何模型：

其中s＝1,2,……num_CCD，为探元编号，(ψ_x(s),ψ_y(s))为探元s在相机坐标系下的指向角，λ＝H/f为摄影比例因子，

为第l行相方坐标系到物方坐标系的旋转矩阵，(X_g,Y_g,Z_g)为探元s物方坐标系下的坐标，(x₀,x₁,x₂,x₃,y₀,y₁,y₂,y₃)为内定标参数模型的系数；

第l行任一探元s，将仿真场景高程数据平均值为Z₀代入成像几何模型计算得到该探元初始地面点物方平面坐标(X_g',Y_g')，在仿真场景高程数据中读取(X_g',Y_g')对应的高程Z_g'，再由成像几何模型计算获得新的平面坐标(X_g",Y_g")，在仿真场景高程数据中读取(X_g",Y_g")对应的高程Z_g"，迭代计算直到本次高程计算结果与上一次高程计算结果两次计算得到的高程之差小于0.1米，计算结束，将本次计算结果作为得到探元s所对应的物方坐标(X_g,Y_g,Z_g)；

遍历所有探元得到每个探元对应的大小为num_line×num_CCD×3的物方坐标三维矩阵RG；

6)以GSD大小为采样间隔实地测得仿真场景中各点三维点的反射率，得到仿真场景反射率数据，按照RG中每个点的三维坐标重采样得到其对应的反射率值，根据地理位置获取太阳辐照度及大气透过率，太阳辐照度乘以大气透过率再乘以反射率，得到各点入瞳辐亮度Ri，输出仿真场景入瞳辐亮度数据；

7)给定探测器最大辐亮度Rmax和最小辐亮度Rmin，杂光系数SLC，渐晕系数VC，面遮拦因子BC，光学系统平均透过率PR，大气MTF矩阵aMTF，光学系统MTF矩阵oMTF，仿真场景入瞳辐亮度数据中各点值Ri到探测器的能量QP＝(Ri×PR×SLC×VC×BC)*aMTF*oMTF，其中*为卷积符号，若QP≥Rmax，则QP＝Rmax，若QP≤Rmin，则QP＝Rmin；

8)给定光电转换效率LER，平均量子效率QE，探测器MTF矩阵dMTF；探测器获得的电子数量NE＝(QP×LER×QE)*dMTF，设定饱和电子数fNE，若NE≥fNE，则NE＝fNE，若NE＜fNE，则NE＝(QP×LER×QE)*dMTF；设定光谱响应函数fs(NE)，量化噪声矩阵QNoise，探测器获得的DN值DN＝fs(NE)*QNoise；

9)输出结果图像。

本发明与现有技术相比的优点在于：首先建立探测器各探元成像几何模型，由几何模型获取成像位置三维坐标，再由三维坐标确定该探元对应的仿真场景反射率，考虑太阳辐照度和大气透过率，在经过光学仿真系统和探测器仿真系统得到仿真图像。同时考虑几何成像和辐射成像关系，提高空间遥感光学成像仿真精度。

附图说明

图1为本发明的操作流程图；

具体实施方式

下面结合附图对本发明做进一步详细描述，结合图1本发明流程图，概括起来，本方法的实施可以分为三个步骤：

1.成像几何建模

1.1实施例中确定成像起始时间t₀＝0，成像时间t＝144s，时标间隔Δt＝1s，采样点数n＝t/Δt＝144，则生成的姿态文件包含144组姿态参数生成的轨道文件包含144组轨道参数(X_i,Y_i,Z_i)，其中i＝1,2…144；将144个姿态参数拟合为以时间为自变量，姿态参数为因变量的姿态仿真三次函数；将144个轨道参数拟合为以时间为自变量，轨道参数为因变量的轨道仿真三次函数；

1.2探元尺寸p＝10μm，轨道高度H＝500km，焦距f＝5m，则地面采样间隔GSD＝pH/f＝1m。图像积分级数为integral_level＝48，integral_time＝0.3ms，每一行探元数量为10000个，则仿真图像行数num_line＝t/(integral_level×integral_time)＝10000。

第l行成像时刻为t_l＝l×(48×0.3ms)；l行对应的姿态数据和轨道数据(X_l,Y_l,Z_l)由t_l带入步骤1.1中的姿态仿真三次函数和轨道仿真三次函数获得；

成像范围为num_line×num_CCD×GSD＝10000m×10000m。

1.3利用起始行轨道数据(X₁,Y₁,Z₁)中的沿轨方向值Y₁和Y₁+10000×GSD作为仿真场景范围Y方向的上下限，X₁-5000×GSD和X₁+5000×GSD作为仿真场景范围X方向的上下限，得到仿真场景范围。对已有的地面数字高程模型在仿真场景范围内进行分辨率为GSD的重采样，得到仿真场景高程数据。

1.4由式建立第l行每一个探元的成像几何模型。其中s为探元号，为探元s在相机坐标系下的指向角，λ＝H/f为摄影比例因子，为第l行相方坐标系到物方坐标系的旋转矩阵，由第l行姿态数据计算得到，(X_g,Y_g,Z_g)为探元s物方坐标系下的坐标，(X_l,Y_l,Z_l)第l行轨道数据，(x₀,x₁,x₂,x₃,y₀,y₁,y₂,y₃)为内定标参数模型的系数。

第l行任一探元s，将仿真场景高程数据平均值为Z₀代入成像几何模型计算得到该探元初始地面点物方平面坐标(X_g',Y_g')，在仿真场景高程数据中读取(X_g',Y_g')对应的高程Z_g'，再由成像几何模型计算获得新的平面坐标(X_g",Y_g")，在仿真场景高程数据中读取(X_g",Y_g")对应的高程Z_g"，迭代计算直到本次高程计算结果与上一次高程计算结果两次计算得到的高程之差小于0.1米，计算结束，将本次计算结果作为得到探元s所对应的物方坐标(X_g,Y_g,Z_g)；

遍历所有探元得到每个探元对应的物方坐标三维矩阵RG(10000×10000×3)。

2.成像辐射建模

2.1以GSD大小为采样间隔实地测得仿真场景中各点三维点的反射率，得到仿真场景反射率数据，按照RG中每个点的三维坐标双三次重采样得到其对应的反射率值，根据地理位置获取太阳辐照度及大气透过率，太阳辐照度乘以大气透过率再乘以反射率，得到各点入瞳辐亮度Ri，输出仿真场景入瞳辐亮度数据。

2.2给定探测器最大辐亮度Rmax＝75和最小辐亮度Rmin＝8，杂光系数SLC＝0.04，渐晕系数VC＝0.2，面遮拦因子BC＝0，光学系统平均透过率PR＝0.35，大气MTF矩阵aMTF，光学系统MTF矩阵oMTF，仿真场景入瞳辐亮度数据中各点值Ri到探测器的能量QP＝(Ri×PR×SLC×VC×BC)*aMTF*oMTF，其中*为卷积符号，若QP≥Rmax，则QP＝Rmax，若QP≤Rmin，则QP＝Rmin。

2.3给定光电转换效率LER＝11×10^-6，平均量子效率QE＝0.25，探测器MTF矩阵dMTF。探测器获得的电子数量NE＝(QP×LER×QE)*dMTF，设定饱和电子数fNE＝86，若NE≥fNE，则NE＝fNE，若NE＜fNE，则NE＝(QP×LER×QE)*dMTF，设定光谱响应函数fs(NE)，量化噪声矩阵QNoise，探测器获得的DN值DN＝fs(NE)*QNoise。

2.4按照10000×10000每个像素的排列顺序输出DN值文件，即为仿真图像。

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

Claims (1)

1.一种基于辐射几何一体化设计的光学遥感图像仿真方法，其特征在于包括以下步骤：

1)确定成像起始时间t₀，成像时间t，时标间隔Δt，采样点数n＝t/Δt，n取整为自然数，第p个采样点处对应的时间为t₀+p×Δt；以时间为顺序，在每个采样点处仿真得到包括成像时侧摆角、俯仰角和偏航角的姿态参数和成像时以摄影中心三维坐标表示的轨道参数(X_i,Y_i,Z_i)，其中i＝1,2…n，为采样点数；将每个姿态参数拟合为以时间为自变量，姿态参数为因变量的姿态仿真三次多项式；将每个轨道参数拟合为以时间为自变量，轨道参数为因变量的轨道仿真三次多项式；

2)确定光学系统焦距f，积分级数integral_level，积分时间integral_time,每一行探元数量为num_CCD，仿真图像行数num_line＝t/(integral_level×integral_time)，num_line取整数，第l行成像时刻为t_l＝t₀+l×(integral_level×integral_time)；l行对应的姿态数据和轨道数据(X_l,Y_l,Z_l)由t_l带入步骤1)中的姿态仿真三次多项式和轨道仿真三次多项式获得；

3)设定探元尺寸p，轨道高度H，地面采样间隔GSD＝pH/f，成像范围为num_line×num_CCD×GSD；

4)利用起始行轨道数据(X₁,Y₁,Z₁)，得到仿真场景范围为X方向X₁-num_CCD×GSD/2到X₁+num_CCD×GSD/2，Y方向为Y₁到Y₁+num_line×GSD；对已有的地面数字高程模型在仿真场景范围内进行分辨率为GSD的重采样，得到仿真场景高程数据；

5)对第l行建立任一探元的成像几何模型：

其中s＝1,2,……num_CCD，为探元编号，(ψ_x(s),ψ_y(s))为探元s在相机坐标系下的指向角，λ＝H/f为摄影比例因子，

为第l行相方坐标系到物方坐标系的旋转矩阵，(X_g,Y_g,Z_g)为探元s物方坐标系下的坐标，(x₀,x₁,x₂,x₃,y₀,y₁,y₂,y₃)为内定标参数模型的系数；

第l行任一探元s，将仿真场景高程数据平均值为Z₀代入成像几何模型计算得到该探元初始地面点物方平面坐标(X_g',Y_g')，在仿真场景高程数据中读取(X_g',Y_g')对应的高程Z_g'，再由成像几何模型计算获得新的平面坐标(X_g",Y_g")，在仿真场景高程数据中读取(X_g",Y_g")对应的高程Z_g"，迭代计算直到本次高程计算结果与上一次高程计算结果两次计算得到的高程之差小于0.1米，计算结束，将本次计算结果作为得到探元s所对应的物方坐标(X_g,Y_g,Z_g)；

遍历所有探元得到每个探元对应的大小为num_line×num_CCD×3的物方坐标三维矩阵RG；

6)以GSD大小为采样间隔实地测得仿真场景中各点三维点的反射率，得到仿真场景反射率数据，按照RG中每个点的三维坐标重采样得到其对应的反射率值，根据地理位置获取太阳辐照度及大气透过率，太阳辐照度乘以大气透过率再乘以反射率，得到各点入瞳辐亮度Ri，输出仿真场景入瞳辐亮度数据；

7)给定探测器最大辐亮度Rmax和最小辐亮度Rmin，杂光系数SLC，渐晕系数VC，面遮拦因子BC，光学系统平均透过率PR，大气MTF矩阵aMTF，光学系统MTF矩阵oMTF，仿真场景入瞳辐亮度数据中各点值Ri到探测器的能量QP＝(Ri×PR×SLC×VC×BC)*aMTF*oMTF，其中*为卷积符号，若QP≥Rmax，则QP＝Rmax，若QP≤Rmin，则QP＝Rmin；

8)给定光电转换效率LER，平均量子效率QE，探测器MTF矩阵dMTF；探测器获得的电子数量NE＝(QP×LER×QE)*dMTF，设定饱和电子数fNE，若NE≥fNE，则NE＝fNE，若NE＜fNE，则NE＝(QP×LER×QE)*dMTF；设定光谱响应函数fs(NE)，量化噪声矩阵QNoise，探测器获得的DN值DN＝fs(NE)*QNoise；

9)输出结果图像。