The construction method of space-based digital imaging link mathematical model

随着数字化、智能化设备建设的快速发展，构建可用于模拟空间成像的天基数字成像系统尤为重要。系统地研究数字成像系统可为后续数字孪生技术的应用提供先期的理论模型。

国外对天基数字成像链路的研究起步较早，目前已有的工程化仿真软件包括美国Satellite tool kit (STK)^[1]、Digital imaging and remote sensing image generation (DIRSIG)^[2]、Time-domain analysis simulation for advanced tracking (TASAT)^[3]，法国Synthetic environment workbench (SE-Workbench)^[4]、德国Software environment for the simulation of optical remote sensing systems (SENSOR)^[5]等，且均已应用于军事项目实践。中国在该领域的研究起步较晚，例如解放军信息工程大学蓝朝桢^[6]采用Geometry behavior model language(GBML)进行复杂目标光学特性建模，结合Open graphics library (OpenGL)实现空间目标的数字成像过程；天津大学杨晋生等^[7]、中国科学院大学许兴星^[8]均采用STK与OpenGL局部光照算法结合的方法进行星载可见光相机仿真成像；华中科技大学杨长才等^[9]通过STK、3D Studio Max及全局光照光线跟踪算法实现目标表面微面元传递到电荷耦合器件(CCD)镜头处能量值的计算；北京理工大学王昊等^[10]采用Optix并行图形处理器(GPU)光线跟踪引擎进行大场景红外成像仿真；北京航天飞行控制中心以天宫二号伴星对组合体的观测任务为需求^[11]，将Visual C++和OpenGL相结合计算辐射量完成像素着色。文献[12-14]在空间目标成像中采用全局光照算法进行目标辐亮度传输，但并未涉及目标与相机的轨道及姿态特性研究。由此可见，国内大多数研究均是基于国外现有仿真软件的二次开发或特定简化模型的应用性分析，未建立起数字成像链路的系统性分析体系。

文中以相机与目标的轨道参数为基础，在满足目标可见性观测条件下，光线从光源出发，经空间目标反射通过空间传输到相机入瞳面，最后由光电转换及成像调制生成最终探测器输出图像。文中针对整个成像环节构建空间目标成像链路模型，进行天基相机数字化系统设计与集成。通过对可视时间段内目标成像结果的分析，验证了数字成像系统的正确性。该成像链路数理模型的构建对空间目标数字成像技术的发展具有重要的参考价值。

文中构建的数理模型主要包括相机与目标光学观测可见性模型、目标几何与光学特性模型、成像辐射传输模型、光电成像转换模型及成像调制模型。主要技术框图如图1所示，几何传输路径中通过观测可见性模型确定目标相对相机的位姿及可视时间段；能量传输路径通过目标几何与光学特性模型、成像辐射传输模型、光电能量转换及成像调制模型确定目标在相机像元面的成像灰度值。

Figure 1.  Technical diagram of the space-based digital imaging system for planar targets

卫星位姿解算的前提是确定各坐标系之间的数学关联。文中以地心赤道惯性坐标系${O_i}XYZ$($J2000.0$)作为基准坐标系构建目标、相机及太阳的位置关系。通过公式(1)将轨道六根数转换为目标$J2000.0$坐标系下的位置矢量。$J2000.0$和质心轨道坐标系$O{X_o}{Y_o}{Z_o}$的位置矢量转换关系如公式(2)所示。$O{X_o}-{Y_o}{Z_o}$坐标系和卫星本体坐标系$O{X_b}{Y_b}{Z_b}$的位置矢量转换关系如公式(3)所示。${O_i}XYZ$、$O{X_o}{Y_o}{Z_o}$和$O{X_b} {Y_b}{Z_b}$坐标系三者关系如图2所示。

式中：$a$为轨道长半轴；$e$为轨道偏心率；$\varOmega $为升交点赤经；$i$为轨道倾角；$u = \omega + \xi $；$\omega $为近地点幅角，$\xi $为真近点角；${R_x}({\textit{Λ}} )$、${R_y}({\textit{Λ}} )$和${R_z}({\textit{Λ}} )$分别为绕$ X、Y、Z $轴旋转${\textit{Λ}} $角(逆时针为正)的基元三维旋转矩阵^[6]；${r_{eo}}$为卫星到地心的距离。

式中：$ {T_{eo}} = \left[ {\begin{array}{*{20}{c}} 0&1&0 \\ 0&0&{ - 1} \\ { - 1}&0&0 \end{array}} \right]{R_z}(u){R_x}(i){R_z}(\varOmega ) $，为惯性坐标系$J2000.0$到质心轨道坐标系$O{X_o}{Y_o}{Z_o}$的坐标旋转矩阵；${L_{eo}} = {\left[ {\begin{array}{*{20}{c}} 0&0&{{r_{eo}}} \end{array}} \right]^{\rm{T}}}$。

式中：$\theta $为俯仰角；$\varphi $为滚动角；$\psi $为偏航角。

Figure 2.  Illustration of coordinate position relationships

太阳相对于地球呈椭圆轨道运动，计算太阳平均轨道要素如表1所示。

首先将太阳轨道根数通过公式(1)转换到$J2000.0$坐标${x_1}{y_1}{z_1}$，再经公式(4)计算太阳平赤经${\alpha _s}$和平赤纬${\delta _s}$。

根据黑体辐射理论，太阳在波长[${\lambda _1},{\lambda _2}$]内的辐射出射度${M_e}$为：

式中：${c_1} = 3.742 \times {10^{ - 16}}\;{\text{ W}} \cdot {{\text{m}}^{\text{2}}}$为第一黑体辐射常数，${c_2} = 1.438\;8 \times {10^{ - 2}}\;{\text{ m}} \cdot {\text{K}}$为第二黑体辐射常数。通常太阳辐射可认为是温度为${T_0} = 5\;900{\text{ K}}$的黑体辐射。

假设太阳在空间各个方向上均匀辐射能量，根据距离平方反比定律，波长450~850 nm范围内太阳在空间目标距离处的辐照度${E_0}$为：

式中：$I$为太阳的发光强度；$\varPhi $为太阳的辐射总通量；${R_s} = 6.959\;9 \times {10^8}\;{\text{m}}$为太阳半径；$r$为太阳与空间目标的距离，由于目标与太阳的距离相对目标与地球的距离大的多，因此可将其视为地日距离$r = 1.521 \times {10^{11}}\;{\text{m}}$。

空间目标的可视条件包括地球遮挡与地光可见$ {G_{E{\text{arth}}}} $、地影可见$ {G_{Eclipse}} $、日光可见${G_{sun}}$、相机设备可见$ {G_P} $等，其几何通视关系如图3所示。

Figure 3.  Diagram illustrating line-of-sight position relationships

探测器的实际可观测区域为${G_{view}} = {G_{Earth}} \cap {G_{Eclipse}} \cap {G_{sun}} \cap {G_p}$，各区域计算如表2所示。

对于天基成像场景，光与目标的相互作用仅发生在物体表面，同时考虑数据兼容性问题，文中研究主要通过点表、面表的多边形表示法构建三维物体的几何信息和拓扑信息，即obj文件表示法，并且该格式表示法有利于编程中对几何形体的解析。

成像系统空间分辨率指标对目标模型的空间颗粒度提出要求，同时空间目标表面结构、包覆材料等特性复杂，一般较难给出其表面准确的反射率表示。因此目标三维几何结构属性采用离散化网格三角面元几何模型表示。

文中研究针对距离函数表示的圆滑曲面几何边界，首先进行概率筛选布点，通过Delaunay构建三角网格，基于桁架结构力平衡原理的均匀化光顺算法^[15]将边界外的点通过梯度函数拉回边界。针对分段表示的几何形体，通过“从外向内”推进波前技术^[16]，采用最小角优先生成原则，解决距离函数进行剖分时棱边脚点不收敛的问题，不同几何表面网格划分结果如图4所示。

Figure 4.  Partitioning of geometric surface mesh based on different distance functions. (a) Square; (b) Sphere; (c) Circular ring; (d) Cube

目标材料光学特性是影响成像速度和成像质量的重要因素，目前主要采用基于双向反射分布函数(BRDF)来表征其光学散射特性。文中目标材质采用的微表面五参量BRDF模型如下：

式中：${\alpha _r}$为微平面法线与平均法线夹角；${\gamma _r}$为入射光矢量与微平面法线夹角；$G({\theta _i},{\theta _r},{\varphi _r})$为遮蔽函数；${k_b},{k_d},{k_r}, {a_r},{b_r}$为待定参数，不同材质拟合参数不同，参量详细说明见文献[17]。

常见的空间目标表面主要由太阳能电池板和聚酰亚胺包覆层组成，其中包覆层材质表现为褶皱的微表面几何特性，并非符合漫反射或者镜面反射。表3给出文中使用的硅电池板材质和聚酰亚胺包覆层针对五参数模型的拟合参数值^[7]。

由上述1.1节可得目标、相机及光源在$J2000.0$坐标系下的位置，通过矩阵${T_{ec}}$将其转换到相机坐标系。相机像素与目标表面点的位置由相机焦距、像元尺寸及目标位姿决定，通过矩阵${T_{ep}}$将目标在相机坐标系下的空间点映射到像素坐标。

式中：${R_x},{R_y},{R_z}$为光源或目标在$J2000.0$坐标系下的位置；${\varTheta _x} = \arccos ({{\boldsymbol{r}}_x}/\sqrt {{{\boldsymbol{r}}_x}^2 + {{\boldsymbol{r}}_y}^2} ) \cdot {\rm sign}({{\boldsymbol{r}}_y})$；${\varTheta _z} = \arccos ({{\boldsymbol{r}}_z}/ \sqrt {{{\boldsymbol{r}}_x}^2 + {{\boldsymbol{r}}_y}^2 + {{\boldsymbol{r}}_z}^2} )$，$({{\boldsymbol{r}}_x},{{\boldsymbol{r}}_y},{{\boldsymbol{r}}_z})$为目标相对于相机的方向矢量。

式中：$({u_0},{v_0})$为图像中心点对应的像素坐标；dx,dy像素大小；$f$相机焦距；$Z'$物点距离。

由渲染公式(10)可知，基于BRDF的渲染方程属于困难积分，无法精确计算其原函数，因此为获得更加逼真的目标图像，文中采用基于蒙特卡洛路径追踪算法的全局光照技术。在不考虑自发辐射的前提下，可将渲染方程改写成公式(11)。

式中：$L(p,{\omega _o})$为p点沿${\omega _o}$方向出射的辐射亮度；${\omega _i}$为入射光方向；${\omega _o}$为出射方向；f表示为BRDF的值；${\theta _i}$为入射方向和表面法线的夹角；${L_e}(p,{\omega _o})$为p点指向${\omega _o}$方向的自发辐亮度，在天基成像场景中可忽略；${L_i}(p,{\omega _i})$为p点沿${\omega _i}$方向入射的辐射亮度。

式中：$N$为单个像元的采样频率；$PDF({\omega _i})$为入射光线概率密度函数。

对于简单随机路径追踪，从相机出发的光线经多次弹射后由于随机采样未击中光源面，则该路径对入瞳面的辐亮度无贡献，计算效率低。因此对光源方向进行重要性采样以提高图像渲染速度，其辐射传输路径如图5所示。

Figure 5.  Schematic diagram of light source importance sampling for radiative transfer path

针对天基成像场景，因目标和光源距离较远，可将太阳及地反光源视为方向光源。在光源方向检查可见性，将渲染方程改为直接光照与间接光的集合，如公式(12)所示：

式中：$M$为成像场景中的光源数；$V(p,{q_m})$为点$p$和光源${q_m}$的可见函数，由光线碰撞检测计算得到，两点可见时为1，否则为0。

路径追踪过程中光线和大多数三角形不相交，使用层次包围盒(BVH)根据光线和场景的三维空间关系减少冗余求交。其时间复杂度为${ O}({\log _2}n)$，图(6)所示为不同递归深度包围盒数量差异。

Figure 6.  Illustration of the number of bounding boxes for different layer depths. (a) Layer depth 1; (b) Layer depth 2; (c) Layer depth 4; (d) Layer depth 6

由于太阳能电池板和聚酰亚胺包覆层反射能量主要集中在镜面反射方向±15°范围内^[7]，为了提高收敛速度，光线击中点的反射方向需在该抽样区间进行采样。采样方向的概率密度函数为$f(\omega ) = 1/(2\pi (1 - \cos {\theta _1}))$，将其转换到球面坐标得$(\theta ,\phi )$的联合密度概率函数为$f(\theta ,\phi ) = \sin \theta /(2\pi (1 - \cos {\theta _1}))$，${\theta _1} = 15^\circ $。根据边缘及条件概率密度函数$f(\phi ),f(\theta \left| \phi \right.)$可得$\theta ,\phi $在球面坐标下的采样。最后通过逆分布函数获得抽样方向的三维坐标如公式(14)所示。

式中：$w = 1/(1 - \cos {\theta _1})$，${\varepsilon _1},{\varepsilon _2}$为$\left[ {0,1} \right]$均匀分布的随机数。

由1.3节计算可得面目标在相机入瞳处的辐亮度$L$，根据公式(15)计算到达像面的光子数${N_{o - }}$。

式中：$D/f'$为相机的相对孔径；$T$为曝光时间；${\tau _0}$为光学系统透过率；$S$为单个像素像元面积；$h$为普朗克常量；$v$为光波频率。

经光电转换和量化像元的灰度值表示为${G_{gray}} = ({N_{o - }}\eta /{N_{full}}) \times {2^{GRD}}$，${N_{full}}$为饱和电子数，$GRD$为量化位数，$\eta $为量子效率。文中选取的成像仿真参数如表4所示。

成像系统的空间调制传递特性在最终成像结果上主要表现为图像模糊，主要的空间调制效应包括相机的像差、离焦，成像平台的线性运动及高频、随机抖动，探测器的采样频率等。根据线性滤波理论，光电成像链路系统可视为线性时空不变系统^[18-19]，分别建立系统各环节的传递函数模型并在频域进行级联相乘可得到系统整体传函，即各效应对成像结果的影响可表示为调制传递函数(MTF)在频域上对二维图像的乘积。成像过程中的噪声主要包括光子噪声、暗电流噪声、探测背景噪声、读出噪声、量化噪声等，各噪声对成像的影响可看作泊松和高斯分布函数在频域上的叠加。加性噪声及线性MTF对成像链路的影响如图7所示。

由傅里叶变换原理及表4参数可知，频域采样间隔$\Delta d$及截止频率${f_c}$分别为如下式所示：

式中：${N_s}$为时域信号采样点数；${F_s}$为时域信号的采样频率。

Figure 7.  Modulation transfer function and mathematical representation of noise in imaging process

文中天基数字成像系统由Visual Studio 2022开发，主要用于实现对天基面目标的模拟成像。首先输入成像平台和目标的轨道六参数，通过光学观测可见性模型获取目标的可成像时间段；其次导入目标三维几何和BRDF材质拟合参数，经成像辐射传输模型计算相机入瞳处辐亮度；基于输入的镜头及探测器参数，根据光电能量转换及成像调制模型，最终输出空间面目标的数字成像结果。

文中通过1.1节可见性模型构建方法，由表5给出的轨道六根数计算相机及目标在$J2000.0$坐标系下的位置矢量，与STK中Two Body模式仿真结果对比如图8所示。24 h及15 d内位置误差在±0.005 m、±0.02 m的波动范围。

根据太阳位置计算模型，获取2022年1~6月1日零时太阳位置，与天文年历表结果对比如表6所示，可得太阳位置的误差在秒量级。

Figure 8.  (a), (b) Represent the position error values between the calculated results of the camera and target for 24 hours and 15 days, respectively, using the mathematical model, and the simulated results from STK

对目标进行几何通视、地影、探测距离（300 km内）条件下15 d可见性分析并与STK仿真结果进行对比，由图9可知，可视时间段无误差。由以上分析可知，相机、目标可见性模型构建过程中坐标转换及位置计算方法正确。

Figure 9.  Results of target visibility time periods within 15 days. (a) Visibility model calculation result; (b) STK simulation result

根据表5、6所示的轨道及成像参数设置，在可见时间段内对空间目标进行数字成像。目标姿态为对地定向，光源设置为太阳光和地反光。太阳辐照度计算如1.1.2节所示，设地球为反射率为0.3的反射体，计算得地球的反射辐照度为190.39 W/m²。

在可视时间段内每间隔3 s的成像结果如图10所示。由图可知，可见时间段内目标的相对位置、姿态及像元辐亮度表现符合真实空间场景成像过程。为观察成像阴影投射的准确性，调整目标的成像距离，不同姿态、不同光照下的成像结果如图11所示，由图中红色框图可以看出成像阴影显示正确，黄色框图体现出全局光照算法的多次弹射结果。

Figure 10.  Imaging results at different distances within the visible time period

Figure 11.  Imaging results under different poses and lighting directions

通过辐射传输模型获取目标辐亮度图像，根据1.4节中空间频率传递特性及噪声模型，对目标辐亮度图像进行处理生成最终的传感器输出图像。由图7可知导致图像质量退化的因素较多，以下给出成像平台高频振动及光子噪声对图像质量的影响分析及仿真结果。

成像平台高频振动MTF可看作零阶贝塞尔函数，其表达式如(18)所示，三维示意图及截止频率内的频谱图如图12所示。高频振动中不同振幅对MTF的影响如图13所示，由图可知振动频率越高，MTF下降越快。图14(b)为高频振动作用后的成像结果。

Figure 12.  (a) Three-dimensional schematic of the high-frequency vibration transfer function MTF for the imaging platform; (b) Spectrum plot within the cutoff frequency

Figure 13.  The influence of different amplitudes of high-frequency vibrations on MTF

式中：${f_v}$为空间频率；$A$为振幅。

光子噪声属于白噪声，服从泊松分布，其方差${n_s} = \sqrt {{N_s}} $，图14(c)、(d)分别为添加光子噪声及光子噪声和高频振动共同作用后的成像结果。

由添加高频振动及光子噪声后的成像结果可以看出，对目标辐亮度图像通过频域MTF乘积和时域噪声累加可以实现图像像质的退化。整个成像过程是多环节耦合的结果，尽可能全面的考虑各环节的MTF和噪声可使仿真结果更逼真。

Figure 14.  (a) Target radiance image; (b) Image with added high-frequency platform vibrations; (c) Image with added photon noise; (d) Image with the combined effects of high-frequency vibrations and photon noise

文中从相机及目标的轨道参数出发，经相机及目标可见性模型、目标几何和光学特性模型、成像辐射传输模型、光电能量转换和成像调制模型，最终输出空间面目标数字成像结果。仿真结果表明：成像平台及目标15 d内的位置误差在±0.02 m波动范围，太阳位置误差为秒量级，通过对比STK可视时间段，验证了相机及目标可见性模型的准确性；模拟成像结果中目标姿态及像元辐亮度表现符合真实空间成像过程；目标成像阴影投射准确，符合光线多次弹射效果；通过成像调制模型可模拟在轨成像中不同程度的像质退化效果。文中构建的成像链路数理模型可实现在不依赖于其他轨道及成像软件如STK、OpenGL前提下开展特定成像场景的数字化成像技术研究，为实体相机参数设计、探测器选型及天基成像数字孪生系统核心模块构建提供参考和依据。

后续可开展的工作包括：1)进行空间背景模型构建，添加恒星库，计算相机视场内恒星的在探测像面的位置及灰度值；2)研究无参图像质量评价算法，对仿真图像成像结果与实测图像进行客观质量评价。