CN104240292A - 一种交会对接逼近段成像敏感器非目标成像的仿真方法 - Google Patents

Abstract

本发明一种交会对接逼近段成像敏感器的非目标成像仿真方法，适用于针孔模型的光学系统的杂散光分析。该方法在应用迭代算法前，须导入具有分组属性的表面数据并预先计算每个待切分面元的切分信息；在迭代算法中，采用了双线程架构，并设计了包含BRDF模型输入输出结果和互遮挡控制代码的结构化面元数据，一个线程利用迭代细分算法产生这种面元数据，另一个线程以之为输入并基于辐射几何得到单个像元接收的辐照度。另外，针对杂散光计算时微面元的互遮挡问题，提出了具体可行的伪投影成像算法。本发明方法模拟的仿真图像可给出非目标成像细节，利用人眼能够直接辨认，而现有的杂散光分析软件尚未能实现这一点。

Description

一种交会对接逼近段成像敏感器非目标成像的仿真方法

技术领域

本发明属于空间光学系统的杂散光分析领域，涉及一种交会对接逼近段成像敏感器的非目标成像的仿真方法。

背景技术

针对交会对接最后逼近段光学成像敏感器的非目标成像模拟，现有的基于随机性算法的杂散光分析软件尚未能提供非目标成像的细节信息，本质上是由于随机性算法在具体应用时其成像精度受限，而目前基于确定性算法的杂散光分析软件的代码实现能力较差。现有的杂散光分析软件在模拟近距离交会对接敏感器成像方面，成像精度有限，不满足实际任务中需提供单像素级的图像模拟需求。

发明内容

本发明解决的技术问题是：克服现有技术的不足，提出了一种交会对接逼近段成像敏感器非目标成像的仿真方法，解决了传统杂散光分析软件未能提供非目标成像细节的缺陷，对全天时交会对接安全时段的选择具有参考价值。

本发明的技术方案是：一种交会对接逼近段成像敏感器的非目标成像仿真方法，步骤如下：

1)获得具有分组属性的目标航天器前端面三维表面数据；根据目标航天器前端面上各部件表面材质，为每个分组选择BRDF模型及模型系数；

2)根据基于边线中点的迭代面元细分算法，在交会对接逼近段过程中的各个位置对步骤1)获得的三维表面数据中的每个面元进行划分，每次划分将每个面元分为4ⁿ个子面元，其中n为细分次数，n的具体数值根据成像敏感器的相机参数和两对接航天器的相对距离确定；记录每个面元每次划分的细分次数n；

3)根据预设的航天器对接时间以及两航天器的轨道六要素和姿态误差，计算得到太阳以及成像敏感器的相机在目标航天器对接面坐标系下的坐标；

31)建立交会对接常用的三个坐标系{d₁},{d₂},{d_m}，分别对应目标航天器对接面坐标系、追踪航天器对接面坐标系和成像敏感器相机坐标系；

32)将坐标系{d₁}下的坐标转换到坐标系{d_m}下，具体的转换关系如下：设两航天器在各自坐标系下的姿态误差为R_ierr(i＝1,2)，即{d_ierr}＝R_ierr{d_i}，坐标系之间的转换由下式表示为：

{d_m}＝R_m{d_1err}+t_m,R_m＝t₂ ^TR_2errR_1err ^-1,t_m＝-t₂ ^T(R_2errM₁+M₂)；

其中，t₂＝[0 0 1；1 0 0；0 1 0]为常置换矩阵，M₁＝[x_r1,y_r1,z_r1]^T为坐标系{d₁}下两航天器的相对位置，M₂＝[x_m,y_m,z_m]^T为成像敏感器相机在坐标系{d₂}下的安装位置；记两航天器各自的姿态角误差为Δ_i＝[Δ_ix,Δ_iy,Δ_iz]^T(i＝1,2)，记Δθ_i ^×为Δ_i对应的反对称矩阵，则有R_ierr＝I₃-Δθ_i ^×，I₃为单位矩阵；

33)根据开普勒定律获得太阳在目标航天器本体系{d_o}下的矢量坐标为

q_s ^o＝L₀₀₁R_z(ω_o+M₀)R_x(i_o)R_z(Ω_o)R_z(-Ω_s)R_x(-i_s)R_z(-ω_s)p₀；

其中，L₀₀₁＝[0 1 0；0 0 -1；-1 0 0]，(a_o,e_o,i_o,Ω_o,w_o,M_o)为目标航天器的轨道六要素，(a_s,e_s,i_s,Ω_s,w_s,M_s)为太阳的轨道六要素；p₀＝[cosv,sinv,0]^T，ν为真近点角，R表示绕坐标轴旋转的方向余弦矩阵；由于{d_o}与{d₁}只存在平移关系，因而太阳在坐标系{d₁}下的矢量坐标也即为q_s ^o；

4)根据步骤2)获得的每个面元每次划分的细分次数n，采用双线程操作，实现整个迭代算法过程；

其中一个线程利用迭代面元细分算法获得子面元结构化面元数据，所述的结构化面元数据包括子面元的顶点坐标、法向量坐标、离轴角、面积，以及子面元到成像敏感器的相机的距离；结合该线程获得的子面元结构化面元数据法向量坐标和步骤3)的结果计算获得子面元的BRDF的输入参数；根据求得的输入参数，选用子面元对应的BRDF模型计算模型的输出值；所述BRDF输入参数的具体方法如下：

41)计算获得BRDF模型的输入参数，所述BRDF模型的输入参数包括入射方向的天顶角θ_i与反射方向的天顶角θ_r，以及入射方向的方位角φ_i与反射方向的方位角φ_r；具体方法如下：

设子面元由A,B,C三个顶点围成的三角形表示，O_i为三角形的重心，对应O_i在坐标系{d₁}下的位置矢量。在坐标系{d₁}下，定义为入射方向和反射方向的单位向量，则

${\hat{p}}_i=-R_{1\mathrm{err}}{q_s}^o,{\hat{p}}_r=-{R_m}^{-1}\frac{{\stackrel{\→}{O}}_{\mathrm{rm}}}{\left|\right|{\stackrel{\→}{O}}_{\mathrm{rm}}\left|\right|},{\mathrm{ns}}_i=\frac{\stackrel{\→}{\mathrm{AB}}\·\stackrel{\→}{\mathrm{AC}}}{\left|\right|\stackrel{\→}{\mathrm{AB}}\left|\right|\left|\right|\stackrel{\→}{\mathrm{AC}}\left|\right|};$

${\mathrm{\θ}}_i={\cos }^{-1}(-{\hat{p}}_i\·{\mathrm{ns}}_i),{\mathrm{\θ}}_r={\cos }^{-1}(-{\hat{p}}_r\·{\mathrm{ns}}_i);$

其中，为矢量的转换坐标，公式如下：

${\stackrel{\→}{O}}_{\mathrm{rm}}=R_m{R}_{\mathrm{lerr}}{\stackrel{\→}{O}}_i+t_m,R_m={t_2}^T{R}_{2\mathrm{err}}{R_{1\mathrm{err}}}^{-1},t_m=-{t_2}^T(R_{2\mathrm{err}}{M}_1+M_2);$

建立该三角形子面元的本地坐标系，三角形重心O_i为坐标系原点，三角形法线方向ns_i为z轴方向，x轴指向y轴由右手坐标系确定； $\stackrel{\→}{\mathrm{OA}}=(\stackrel{\→}{B}+\stackrel{\→}{C}-2\stackrel{\→}{A})/3,$ 则

$\hat{x}=\frac{\stackrel{-}{O}\stackrel{\→}{A}}{\left|\right|\stackrel{-}{O}\stackrel{\→}{A}\left|\right|},{\cos \mathrm{\φ}}_i=\frac{(-{\hat{p}}_i-\hat{z}\cos {\mathrm{\θ}}_i)\·\hat{x}}{\left|\right|-{\hat{p}}_i-\hat{z}\cos {\mathrm{\θ}}_i\left|\right|},{\cos \mathrm{\φ}}_r=\frac{({\hat{p}}_r-\hat{z}\cos {\mathrm{\θ}}_r)\·\hat{x}}{\left|\right|{\hat{p}}_r-\hat{z}\cos {\mathrm{\θ}}_r\left|\right|};$

42)根据求得的输入参数，选用子面元对应的BRDF模型计算模型的输出值；

另一个线程根据获得的子面元结构化面元数据，结合辐射几何、BRDF模型的输出值计算获得成像敏感器相机单个像元接收的低精度辐照度，经剔除重复后得到成像敏感器相机成像平面上单个像元接收的高精度辐照度；

5)整个迭代过程完毕后，形成完整的仿真图像。

所述步骤4)中剔除重复采用动态剔除方法，具体步骤如下：

a)反射投影图像计算

$x_p=\frac{-f^{\′}{x}_{\mathrm{mi}}}{z_{\mathrm{mi}}},y_p=\frac{-f^{\′}{y}_{\mathrm{mi}}}{z_{\mathrm{mi}}}$

根据中心投影，有 $x_u=\frac{\mathrm{sizeX}-1}{2}+\frac{x_p}{{\mathrm{dp}}_x},y_v=\frac{\mathrm{sizeY}-1}{2}+\frac{y_p}{{\mathrm{dp}}_x}$

其中，u∈[0,sizeX],v∈[0,sizeY]，表示对x按照下边界进行取整，[sizeX,sizeY]为成像敏感器相机CCD的尺寸，dpx为成像敏感器相机CCD单个像元的尺寸，(x_p,y_p)为三角形子面元dS_i的重心在坐标系{d_m}下的坐标，(x_u,y_u)为相机坐标系{d_m}对应的像素坐标系下的坐标；

以中心投影的距离作为反射投影图像的像素值，则反射投影图像的单个像元的像素值为其中，dS_mi＝[x_mi,y_mi,z_mi]^T记为三角形子面元dS_i在坐标系{d_m}下的坐标；

b)计算获得入射投影图像

b1)建立投影坐标系{d_p}，所述坐标系{d_p}为从坐标系{d₁}绕欧拉轴角(a_u,φ_u)旋转后得到的坐标系，其中坐标系{d_p}的x轴指向太阳方向，投影坐标系{d_p}原点到坐标系{d₁}原点的距离与坐标系{d₁}原点到坐标系{d_m}原点的距离相等；其中 ${\mathrm{\φ}}_u={\cos }^{-1}[-{\hat{p}}_i^T{a}_x],$ $a_u=a_x\×(-{\hat{p}}_i),$ a_x＝[1,0,0]^T；

b2)计算获得欧拉轴角(a_u,φ_u)对应的方向余弦矩阵M_proj＝cosφ_uE+(1-cosφ_u)a_ua_u ^T-sinφ_ua_u ^×，其中E是单位矩阵；

b3)将三角形子面元dS_i由坐标系{d₁}转换到坐标系{d_p}，即dS_pi＝M_projdS_i+d_M，其中，d_M＝[x_r1,0,0]^T；记dS_pi的x坐标为x_pi；

b4)按中心投影计算获得dS_pi的像素位置，计算过程中用x_r1代替x_pi；

b5)获得dS_pi像素位置的像素值I_ip(u,v)＝x_pi；

c)判断参与计算的子面元dS_i在太阳方向和成像敏感器相机的视点方向是否满足互遮挡条件是否满足互遮挡条件

根据步骤a)和步骤b)计算dS_pi和dS_mi对应的像元位置的当前像素值I_cur1和I_cur2，查询由迭代算法之前获得的dS_pi和dS_mi对应的像素位置I_ip(u,v)和I_rp(u,v)，当I_cur1>I_ip(u,v)或I_cur2>I_rp(u,v)时，认为满足互遮挡条件；

d)在每次迭代过程中，如满足互遮挡条件，则忽略当前计算出的低精度的辐照度值，否则将当前计算出的低精度辐照度值加到当前像元的像素值上。

本发明与现有技术相比的有益成果是：

1)本发明方法给出的仿真图像可显示非目标成像细节，利用人眼能够直接辨认，而现有的杂散光分析软件尚未能实现这一点。

2)本发明方法对应的确定性算法同时保持了快速计算和稳定性，与传统的杂散光分析软件相比，计算效率更高(前者单个算例耗时在二十几秒，后者在1到2小时之间。)

3)本发明方法对应的确定性算法具有良好的扩展性和实用性，可扩展应用到交会对接最后逼近段成像敏感器的地球反照成像，可切换为不同的面元切分算法。

附图说明

图1为本发明的非目标成像仿真方法的概略图；

图2为目标航天器前锥段各类表面的分布示意图；

图3为交会对接常用坐标系示意图；

图4为结构化表面微元的流程图；

图5为BRDF几何示意图；

图6为三角形面元当地坐标系示意图；

图7为离轴像点的辐照度计算示意图；

图8为伪投影成像方法流程图。

具体实施方式

如图1所示，为本发明方法所涉及的确定性算法的概略图。

本发明中，首先导入.OBJ文件格式的三维表面数据，并添加表面数据分组属性，为应用三角形面元迭代细分算法准备前提条件。为了实现快速计算和稳定性，需要保持迭代算法的独立性，从而引入了双线程架构，并设计了结构化的面元数据为双线程同步计算提供便利条件。为了进一步提高该确定性算法的准确性，针对杂散光计算时微面元互遮挡引起的重复计算问题提出了伪投影成像算法。

本发明方法的具体步骤为：

(1)在应用迭代算法前，导入图2所示的前端面表面数据作为迭代算法的输入。由于实际导入的三维数据为.OBJ文件格式(该格式存储三角形面元数据，优点是读取容易、程序扩展性好，改进后各分组数据可实现物理隔离存储，增加访问效率)，如果导入的数据没有分组属性，设置各个分组的搜索条件，并改进已有的OBJ文件格式以添加分组属性。然后，为每个分组选择合适的BRDF模型及模型系数。接着，为每个待切分的面元增加保证单个像元空间采样精度的切分次数信息。

对于图2所示的各种不规则表面，最后选用了5种表面材质，根据表面材质类型归结为3种BRDF模型(CookTorrance模型、Davies模型和5参数模型)，分别适用于表面镀膜、镁铝合金或铝板以及表面涂白漆的表面材质，各BRDF模型的相关系数选择见参考文献[1-3]。

[1]ShackR.V.,HarveyJ.E.,“Light-Scattering Characteristics of Optical Surfaces”,SPIE Vol.107,41-47,1977

[2]Bailey,W.Albert,Edward A.,KepplerK.S.,“Dynamic Bidirectional ReflectancDistribution Functions:Measurement and Representation”,J.LaserAppl.,20:22-36,2008

[3]H.Y.Wang,W.Zhang,F.G.Wang,“Visible characteristics of space-based targetsbased on bidirectional reflection distribution function”,SCIENCECHINA(TechnologicalSciences),55:982-989,2012

(2)在面元迭代细分算法中，为保持切分前后三角形的相似性选用基于边线中点的迭代细分算法(见参考文献[4])。根据该面元迭代细分算法，在交会对接逼近段过程中的各个位置对步骤1)获得的三维表面数据中的每个面元进行划分，每次划分将每个面元分为4ⁿ个子面元，其中n为细分次数，n的具体数值根据成像敏感器的相机参数和两对接航天器的相对距离确定；记录每个面元每次划分的细分次数n；

[4]Kamen Y.,Shirman L,“Triangle rendering using adaptive subdivision”,IEEEon Computer Graphics and Application,95-103,1998

(3)根据预设的航天器对接时间以及两航天器的轨道六要素和姿态误差，计算得到太阳以及成像敏感器的相机在目标航天器对接面坐标系下的坐标；

31)建立交会对接常用的三个坐标系{d₁},{d₂},{d_m}，分别对应目标航天器对接面坐标系、追踪航天器对接面坐标系和成像敏感器相机坐标系，如图3所示，三种坐标系分别定义如下：

坐标系{d₁}，原点为目标航天器对接机构对接框端面的中心，x轴与目标航天器对称纵轴重合，指向目标航天器飞行方向，y轴沿目标航天器的帆板方向，z轴由右手系决定；

坐标系{d₂}，原点为追踪航天器对接机构捕获环导向瓣根圆圆心，x轴指向追踪航天器的飞行方向，y轴沿追踪航天器的帆板方向，z轴由右手系决定；

坐标系{d_m}，原点为成像敏感器相机的光心，z轴与成像敏感器相机的主光轴重合，其正方向与入射光线反向，y轴与追踪航天器本体系的z轴平行，x轴由右手系决定；

32)将坐标系{d₁}下的坐标转换到坐标系{d_m}下，转换关系如下：设两航天器在各自坐标系下的姿态误差为R_ierr(i＝1,2)，即{d_ierr}＝R_ierr{d_i}，实际的对接坐标系之间的转换由下式表示为：

{d_m}＝R_m{d_1err}+t_m,R_m＝t₂ ^TR_2errR_1err ^-1,t_m＝-t₂ ^T(R_2errM₁+M₂)；

其中，t₂＝[0 0 1；1 0 0；0 1 0]为常置换矩阵，M₁＝[x_r1,y_r1,z_r1]^T为{d₁}坐标系下两航天器的相对位置，M₂＝[x_m,y_m,z_m]^T为成像敏感器相机在坐标系{d₂}下的安装位置；记两航天器各自的姿态角误差为Δ_i＝[Δ_ix,Δ_iy,Δ_iz]^T(i＝1,2)，记Δθ_i ^×为Δ_i对应的反对称矩阵，则有R_ierr＝I₃-Δθ_i ^×，I₃为单位矩阵；

33)根据开普勒定律获得太阳在目标航天器本体系{d_o}下的矢量坐标为

q_s ^o＝L₀₀₁R_z(ω_o+M₀)R_x(i_o)R_z(Ω_o)R_z(-Ω_s)R_x(-i_s)R_z(-ω_s)p₀；

其中，L₀₀₁＝[0 1 0；0 0 -1；-1 0 0]，目标航天器的轨道参数(a_o,e_o,i_o,Ω_o,w_o,M_o)给定，太阳的轨道参数(a_s,e_s,i_s,Ω_s,w_s,M_s)根据给定时刻由开普勒定律直接可求(参考文献[3])，p₀＝[cosv,sinv,0]^T，ν为真近点角，R表示绕坐标轴旋转的方向余弦矩阵；由于{d_o}与{d₁}只存在平移关系，因而太阳在坐标系{d₁}下的矢量坐标也即为q_s ^o；

4)根据步骤2)获得的每个面元每次划分的细分次数n，采用双线程操作，实现整个迭代算法过程；

一个线程利用迭代面元细分算法获得子面元结构化面元数据，所述的结构化面元数据包括子面元的顶点坐标、法向量坐标、离轴角、面积，以及子面元到成像敏感器的相机的距离，如图4所示；结合该线程获得的子面元结构化面元数据法向量坐标和步骤3)的结果计算获得子面元的BRDF的输入参数，如图5所示；根据求得的输入参数，选用子面元对应的BRDF模型计算模型的输出值；所述BRDF输入参数的具体方法如下：

41)计算获得BRDF模型的输入参数，所述BRDF模型的输入参数包括入射方向的天顶角θ_i与反射方向的天顶角θ_r，以及入射方向的方位角φ_i与反射方向的方位角φ_r；θ_i,θ_r在坐标系{d₁}下计算得到。φ_i,φ_r在三角形子面元的本地坐标系计算得到。具体方法如下：

设子面元由A,B,C三个顶点围成的三角形表示，如图6所示。O_i为三角形的重心，对应O_i在坐标系{d₁}下的位置矢量。在坐标系{d₁}下，定义为入射方向和反射方向的单位向量，则

${\hat{p}}_i=-R_{1\mathrm{err}}{q_s}^o,{\hat{p}}_r=-{R_m}^{-1}\frac{{\stackrel{\→}{O}}_{\mathrm{rm}}}{\left|\right|{\stackrel{\→}{O}}_{\mathrm{rm}}\left|\right|},{\mathrm{ns}}_i=\frac{\stackrel{\→}{\mathrm{AB}}\·\stackrel{\→}{\mathrm{AC}}}{\left|\right|\stackrel{\→}{\mathrm{AB}}\left|\right|\left|\right|\stackrel{\→}{\mathrm{AC}}\left|\right|};$

${\mathrm{\θ}}_i={\cos }^{-1}(-{\hat{p}}_i\·{\mathrm{ns}}_i),{\mathrm{\θ}}_r={\cos }^{-1}(-{\hat{p}}_r\·{\mathrm{ns}}_i);$

其中，为矢量的转换坐标，公式如下

${\stackrel{\→}{O}}_{\mathrm{rm}}=R_m{R}_{\mathrm{lerr}}{\stackrel{\→}{O}}_i+t_m,R_m={t_2}^T{R}_{2\mathrm{err}}{R_{1\mathrm{err}}}^{-1},t_m=-{t_2}^T(R_{2\mathrm{err}}{M}_1+M_2);$

建立该三角形子面元的本地坐标系，三角形重心O_i为坐标系原点，三角形法线方向ns_i为z轴方向，x轴指向y轴由右手坐标系确定；则

$\hat{x}=\frac{\stackrel{-}{O}\stackrel{\→}{A}}{\left|\right|\stackrel{-}{O}\stackrel{\→}{A}\left|\right|},{\cos \mathrm{\φ}}_i=\frac{(-{\hat{p}}_i-\hat{z}\cos {\mathrm{\θ}}_i)\·\hat{x}}{\left|\right|-{\hat{p}}_i-\hat{z}\cos {\mathrm{\θ}}_i\left|\right|},{\cos \mathrm{\φ}}_r=\frac{({\hat{p}}_r-\hat{z}\cos {\mathrm{\θ}}_r)\·\hat{x}}{\left|\right|{\hat{p}}_r-\hat{z}\cos {\mathrm{\θ}}_r\left|\right|};$

42)根据求得的输入参数，选用子面元对应的BRDF模型计算模型的输出值；

另一个线程根据获得的子面元结构化面元数据，结合辐射几何、BRDF模型的输出值计算获得成像敏感器相机单个像元接收的低精度辐照度，经剔除重复计算后得到成像敏感器相机成像平面上单个像元接收的高精度辐照度；

所述的单个像元接收的低精度辐照度值的计算公式如下：

${E^{\′}}_w=\mathrm{\π\τS}\left({\mathrm{\θ}}_m\right){\cos }^4{w}^{\′}\cos {\mathrm{\θ}}_i{\∫}_{{\mathrm{\λ}}_1}^{{\mathrm{\λ}}_2}F\left(\mathrm{\λ}\right)R\left(\mathrm{\λ}\right)L_{\mathrm{dSi}}({\mathrm{\θ}}_r,{\mathrm{\φ}}_r,\mathrm{\λ})\mathrm{d\λ}$

$L_{{\mathrm{\Δ}}_{\mathrm{ij}}}({\mathrm{\θ}}_r,{\mathrm{\φ}}_r,\mathrm{\λ})=s\left(\mathrm{\λ}\right)E_0\mathrm{BRDF}({\mathrm{\θ}}_r,{\mathrm{\φ}}_r,\mathrm{\λ})\cos {\mathrm{\θ}}_i$

$w^{\′}={\tan }^{-1}\left(\frac{f^{\′}\tan {\mathrm{\θ}}_m}{l_{p2}}\right),\sin {U_w}^{\′}\≈\frac{D_{p2}{\cos }^2{w}^{\′}}{{2l}_{p2}}$

上式中，E'_w为散射表面元dS_i(也即三角形切分面元Δ_ij)对应的轴外像点ds'的辐照度。其他各参量的意义如下：

τ——整个光学系统的透过率

S(θ_m)——关于离轴角θ_m的渐晕系数

θ_i——太阳矢量与前锥端表面切分后的子面元Δ_ij法向的夹角

F(λ)——为整个光学滤波片组的光谱透过率函数，λ为波长

R(λ)——为光学成像敏感器成像面阵(如CCD)的归一化光谱响应率

——表面元Δ_ij沿(θ_r,Φ_r)方向散射出的光谱辐亮度，(θ_r,Φ_r)的定义与图5相对应，分别对应散射方向的天顶角和方位角

E₀——太阳常数，取1366.9w/m²

s(λ)——为太阳的归一化光谱权重函数

BRDF(θ_r,Φ_r,λ)——前锥端表面切分后的子面元某一波长λ下的反射特性

w',U_w'——定义如图7所示，f'为像方焦距，θ_m为离轴角，l_p2为出瞳面到焦平面的距离，D_p2为出瞳直径。

上面的E'_w不能直接用来计算离轴辐照度。为不失一般性，以λ₁＝587.6nm作为BRDF计算的参考，结合光谱BRDF的转换关系并进一步简化计算，则可以得到

$E_{\mathrm{\ω}}^{\′}=u_0{\cos }^4{w}^{\′}\cos {\mathrm{\θ}}_i{\mathrm{BRDF}}_m({\mathrm{\θ}}_r,{\mathrm{\φ}}_r,{\mathrm{\λ}}_1){\∫}_{0.78}^{1.78}{\mathrm{\λ}}^{m-4}\mathrm{d\λ}$

BRDF_m＝f_r+brdf_Env

上式中，对于该光学系统而言，渐晕在半视场内S(θ_m)＝1.0，u₀＝0.0204，m由表面材质的属性决定。其他参数定义如下

f_r——某个BRDF模型的输出

brdf_Env——前端面各种表面之间发生的二次散射贡献的环境光分量，仿真时假定二次散射为各向同性，并取为0.005。

所述步骤4)中剔除重复采用动态剔除方法(图8概括为伪投影成像方法)，具体步骤如下：

a)反射投影图像计算

$x_p=\frac{-f^{\′}{x}_{\mathrm{mi}}}{z_{\mathrm{mi}}},y_p=\frac{-f^{\′}{y}_{\mathrm{mi}}}{z_{\mathrm{mi}}}$

根据中心投影，有 $x_u=\frac{\mathrm{sizeX}-1}{2}+\frac{x_p}{{\mathrm{dp}}_x},y_v=\frac{\mathrm{sizeY}-1}{2}+\frac{y_p}{{\mathrm{dp}}_x}$

其中，u∈[0,sizeX],v∈[0,sizeY]，表示对x按照下边界进行取整，[sizeX,sizeY]为成像敏感器相机CCD的尺寸，dpx为成像敏感器相机CCD单个像元的尺寸，(x_p,y_p)为三角形子面元dS_i的重心在坐标系{d_m}下的坐标，(x_u,y_u)为相机坐标系{d_m}对应的像素坐标系下的坐标；

以中心投影的距离作为反射投影图像的像素值，则反射投影图像的单个像元的像素值为

$I_{\mathrm{rp}}(u,v)=\sqrt{x_{\mathrm{mi}}^2+y_{\mathrm{mi}}^2+z_{\mathrm{mi}}^2}$

其中，dS_mi＝[x_mi,y_mi,z_mi]^T记为三角形子面元dS_i在坐标系{d_m}下的坐标；

b)计算获得入射投影图像

b1)建立投影坐标系{d_p}，所述坐标系{d_p}为从坐标系{d₁}绕欧拉轴角(a_u,φ_u)旋转后得到的坐标系，其中坐标系{d_p}的x轴指向太阳方向，投影坐标系{d_p}原点到坐标系{d₁}原点的距离与坐标系{d₁}原点到坐标系{d_m}原点的距离相等；其中 ${\mathrm{\φ}}_u={\cos }^{-1}[-{\hat{p}}_i^T{a}_x],$ $a_u=a_x\×(-{\hat{p}}_i),$ a_x＝[1,0,0]^T；

b2)计算获得欧拉轴角(a_u,φ_u)对应的方向余弦矩阵M_proj＝cosφ_uE+(1-cosφ_u)a_ua_u ^T-sinφ_ua_u ^×，其中E是单位矩阵；

b3)将三角形子面元dS_i由坐标系{d₁}转换到坐标系{d_p}，即dS_pi＝M_projdS_i+d_M，其中，d_M＝[x_r1,0,0]^T；记dS_pi的x坐标为x_pi；

b4)按中心投影计算获得dS_pi的像素位置，计算过程中用x_r1代替x_pi；

b5)获得dS_pi像素位置的像素值I_ip(u,v)＝x_pi；

c)判断参与计算的子面元dS_i在太阳方向和成像敏感器相机的视点方向是否满足互遮挡条件

根据步骤a)和步骤b)计算dS_pi和dS_mi对应的像元位置的当前像素值I_cur1和I_cur2，查询由迭代算法之前获得的dS_pi和dS_mi对应的像素位置I_ip(u,v)和I_rp(u,v)，当I_cur1>I_ip(u,v)或I_cur2>I_rp(u,v)时，认为满足互遮挡条件；

d)在每次迭代过程中，如满足互遮挡条件，则忽略当前计算出的低精度的辐照度值，否则将当前计算出的低精度辐照度值加到当前像元(切分子面元dS_i通过成像敏感器相机坐标系{d_m}投影到其成像平面上的像元位置)的像素值上。

5)整个迭代过程完毕后，形成完整的仿真图像。

本发明未具体说明部分属于本领域公知技术。

Claims (2)

1.一种交会对接逼近段成像敏感器的非目标成像仿真方法，其特征在于步骤如下：

1)获得具有分组属性的目标航天器前端面三维表面数据；根据目标航天器前端面上各部件表面材质，为每个分组选择BRDF模型及模型系数；

2)根据基于边线中点的迭代面元细分算法，在交会对接逼近段过程中的各个位置对步骤1)获得的三维表面数据中的每个面元进行划分，每次划分将每个面元分为4ⁿ个子面元，其中n为细分次数，n的具体数值根据成像敏感器的相机参数和两对接航天器的相对距离确定；记录每个面元每次划分的细分次数n；

3)根据预设的航天器对接时间以及两航天器的轨道六要素和姿态误差，计算得到太阳以及成像敏感器的相机在目标航天器对接面坐标系下的坐标；

31)建立交会对接常用的三个坐标系{d₁},{d₂},{d_m}，分别对应目标航天器对接面坐标系、追踪航天器对接面坐标系和成像敏感器相机坐标系；

32)将坐标系{d₁}下的坐标转换到坐标系{d_m}下，具体的转换关系如下：设两航天器在各自坐标系下的姿态误差为R_ierr(i＝1,2)，即{d_ierr}＝R_ierr{d_i}，坐标系之间的转换由下式表示为：

{d_m}＝R_m{d_1err}+t_m,R_m＝t₂ ^TR_2errR_1err ^-1,t_m＝-t₂ ^T(R_2errM₁+M₂)；

其中，t₂＝[0 0 1；1 0 0；0 1 0]为常置换矩阵，M₁＝[x_r1,y_r1,z_r1]^T为坐标系{d₁}下两航天器的相对位置，M₂＝[x_m,y_m,z_m]^T为成像敏感器相机在坐标系{d₂}下的安装位置；记两航天器各自的姿态角误差为Δ_i＝[Δ_ix,Δ_iy,Δ_iz]^T(i＝1,2)，记Δθ_i ^×为Δ_i对应的反对称矩阵，则有R_ierr＝I₃-Δθ_i ^×，I₃为单位矩阵；

33)根据开普勒定律获得太阳在目标航天器本体系{d_o}下的矢量坐标为

q_s ^o＝L₀₀₁R_z(ω_o+M₀)R_x(i_o)R_z(Ω_o)R_z(-Ω_s)R_x(-i_s)R_z(-ω_s)p₀；

其中，L₀₀₁＝[0 1 0；0 0 -1；-1 0 0]，(a_o,e_o,i_o,Ω_o,w_o,M_o)为目标航天器的轨道六要素，(a_s,e_s,i_s,Ω_s,w_s,M_s)为太阳的轨道六要素；p₀＝[cosv,sinv,0]^T，ν为真近点角，R表示绕坐标轴旋转的方向余弦矩阵；由于{d_o}与{d₁}只存在平移关系，因而太阳在坐标系{d₁}下的矢量坐标也即为q_s ^o；

4)根据步骤2)获得的每个面元每次划分的细分次数n，采用双线程操作，实现整个迭代算法过程；

其中一个线程利用迭代面元细分算法获得子面元结构化面元数据，所述的结构化面元数据包括子面元的顶点坐标、法向量坐标、离轴角、面积，以及子面元到成像敏感器的相机的距离；结合该线程获得的子面元结构化面元数据法向量坐标和步骤3)的结果计算获得子面元的BRDF的输入参数；根据求得的输入参数，选用子面元对应的BRDF模型计算模型的输出值；所述BRDF输入参数的具体方法如下：

41)计算获得BRDF模型的输入参数，所述BRDF模型的输入参数包括入射方向的天顶角θ_i与反射方向的天顶角θ_r，以及入射方向的方位角φ_i与反射方向的方位角φ_r；具体方法如下：

设子面元由A,B,C三个顶点围成的三角形表示，O_i为三角形的重心，对应O_i在坐标系{d₁}下的位置矢量。在坐标系{d₁}下，定义为入射方向和反射方向的单位向量，则

${\hat{p}}_i=-R_{1\mathrm{err}}{q_s}^o,{\hat{p}}_r=-{R_m}^{-1}\frac{{\stackrel{\→}{O}}_{\mathrm{rm}}}{\left|\right|{\stackrel{\→}{O}}_{\mathrm{rm}}\left|\right|},{\mathrm{ns}}_i=\frac{\stackrel{\→}{\mathrm{AB}}\·\stackrel{\→}{\mathrm{AC}}}{\left|\right|\stackrel{\→}{\mathrm{AB}}\left|\right|\left|\right|\stackrel{\→}{\mathrm{AC}}\left|\right|};$

${\mathrm{\θ}}_i={\cos }^{-1}(-{\hat{p}}_i\·{\mathrm{ns}}_i),{\mathrm{\θ}}_r={\cos }^{-1}(-{\hat{p}}_r\·{\mathrm{ns}}_i);$

其中，为矢量的转换坐标，公式如下：

${\stackrel{\→}{O}}_{\mathrm{rm}}=R_m{R}_{\mathrm{lerr}}{\stackrel{\→}{O}}_i+t_m,R_m={t_2}^T{R}_{2\mathrm{err}}{R_{1\mathrm{err}}}^{-1},t_m=-{t_2}^T(R_{2\mathrm{err}}{M}_1+M_2);$

建立该三角形子面元的本地坐标系，三角形重心O_i为坐标系原点，三角形法线方向ns_i为z轴方向，x轴指向y轴由右手坐标系确定； $\stackrel{\→}{\mathrm{OA}}=(\stackrel{\→}{B}+\stackrel{\→}{C}-2\stackrel{\→}{A})/3,$ 则

$\hat{x}=\frac{\stackrel{-}{O}\stackrel{\→}{A}}{\left|\right|\stackrel{-}{O}\stackrel{\→}{A}\left|\right|},{\cos \mathrm{\φ}}_i=\frac{(-{\hat{p}}_i-\hat{z}\cos {\mathrm{\θ}}_i)\·\hat{x}}{\left|\right|-{\hat{p}}_i-\hat{z}\cos {\mathrm{\θ}}_i\left|\right|},{\cos \mathrm{\φ}}_r=\frac{({\hat{p}}_r-\hat{z}\cos {\mathrm{\θ}}_r)\·\hat{x}}{\left|\right|{\hat{p}}_r-\hat{z}\cos {\mathrm{\θ}}_r\left|\right|};$

42)根据求得的输入参数，选用子面元对应的BRDF模型计算模型的输出值；

另一个线程根据获得的子面元结构化面元数据，结合辐射几何、BRDF模型的输出值计算获得成像敏感器相机单个像元接收的低精度辐照度，经剔除重复后得到成像敏感器相机成像平面上单个像元接收的高精度辐照度；

5)整个迭代过程完毕后，形成完整的仿真图像。

2.根据权利要求1所述的一种交会对接逼近段成像敏感器的非目标成像仿真方法，其特征在于所述步骤4)中剔除重复采用动态剔除方法，具体步骤如下：

a)反射投影图像计算

$x_p=\frac{-f^{\′}{x}_{\mathrm{mi}}}{z_{\mathrm{mi}}},y_p=\frac{-f^{\′}{y}_{\mathrm{mi}}}{z_{\mathrm{mi}}}$

根据中心投影，有 $x_u=\frac{\mathrm{sizeX}-1}{2}+\frac{x_p}{{\mathrm{dp}}_x},y_v=\frac{\mathrm{sizeY}-1}{2}+\frac{y_p}{{\mathrm{dp}}_x}$

其中，u∈[0,sizeX],v∈[0,sizeY]，表示对x按照下边界进行取整，[sizeX,sizeY]为成像敏感器相机CCD的尺寸，dpx为成像敏感器相机CCD单个像元的尺寸，(x_p,y_p)为三角形子面元dS_i的重心在坐标系{d_m}下的坐标，(x_u,y_u)为相机坐标系{d_m}对应的像素坐标系下的坐标；

以中心投影的距离作为反射投影图像的像素值，则反射投影图像的单个像元的像素值为其中，dS_mi＝[x_mi,y_mi,z_mi]^T记为三角形子面元dS_i在坐标系{d_m}下的坐标；

b)计算获得入射投影图像

b1)建立投影坐标系{d_p}，所述坐标系{d_p}为从坐标系{d₁}绕欧拉轴角(a_u,φ_u)旋转后得到的坐标系，其中坐标系{d_p}的x轴指向太阳方向，投影坐标系{d_p}原点到坐标系{d₁}原点的距离与坐标系{d₁}原点到坐标系{d_m}原点的距离相等；其中 ${\mathrm{\φ}}_u={\cos }^{-1}[-{\hat{p}}_i^T{a}_x],$ $a_u=a_x\×(-{\hat{p}}_i),$ a_x＝[1,0,0]^T；

b2)计算获得欧拉轴角(a_u,φ_u)对应的方向余弦矩阵M_proj＝cosφ_uE+(1-cosφ_u)a_ua_u ^T-sinφ_ua_u ^×，其中E是单位矩阵；

b3)将三角形子面元dS_i由坐标系{d₁}转换到坐标系{d_p}，即dS_pi＝M_projdS_i+d_M，其中，d_M＝[x_r1,0,0]^T；记dS_pi的x坐标为x_pi；

b4)按中心投影计算获得dS_pi的像素位置，计算过程中用x_r1代替x_pi；

b5)获得dS_pi像素位置的像素值I_ip(u,v)＝x_pi；

c)判断参与计算的子面元dS_i在太阳方向和成像敏感器相机的视点方向是否满足互遮挡条件是否满足互遮挡条件

根据步骤a)和步骤b)计算dS_pi和dS_mi对应的像元位置的当前像素值I_cur1和I_cur2，查询由迭代算法之前获得的dS_pi和dS_mi对应的像素位置I_ip(u,v)和I_rp(u,v)，当I_cur1>I_ip(u,v)或I_cur2>I_rp(u,v)时，认为满足互遮挡条件；

d)在每次迭代过程中，如满足互遮挡条件，则忽略当前计算出的低精度的辐照度值，否则将当前计算出的低精度辐照度值加到当前像元的像素值上。