一种大气折射修正状态下目标重建方法

摘要

本发明公开了一种大气折射修正状态下目标重建方法，属于图像处理的三维重建方法，解决现有目标重建方法未考虑大气折射对光线传播影响而导致的目标定位误差较大的问题，可用于远距离预警系统。本发明方法包括：(1)确定过目标和地心的直线方程；(2)根据相机C

说明书

技术领域

本发明属于图像处理技术领域，更具体地，涉及一种大气折射修正状 态下目标重建方法，可用于卫星导弹的预警，提高预警防御的准确性。

背景技术

高技术条件下的未来战争，突发性强，时间概念高度浓缩，对导弹的 预警显得格外重要，需要能够在第一时间得到导弹弹道、方向、速度、发 射点与弹着点及导弹数量等重要信息，以便于防御和反导。卫星导弹预警 是国家预警系统中关键的一部分，具有极其重要的地位与作用。信息处理 是卫星导弹预警系统的核心，主要包括预警卫星红外成像、弱小目标的检 测与识别、轨迹重建研究等。导弹的轨迹重建是导弹卫星预警中的关键技 术，根据导弹目标在像面成像位置，相机参数等信息重建出导弹的轨迹， 以便导弹的拦截预警。现有的轨迹重建技术都是在不考虑大气折射的影响 的情况下，按照光直线传播原理，根据目标像面位置，按照直线进行的轨 迹重建。但在实际的光线传输过程中，由于大气折射的影响，光线的传输 路径会发生偏折，光的传输路径并不是直线。所以现有的目标重建技术得 到的目标轨迹与目标实际的轨迹之间存在着误差，不能较为准确地定位目 标的位置。

发明内容

针对现有的目标重建技术的缺陷和改进需求，本发明提供一种大气折 射修正状态下目标重建技术，与现有的按照光直线传播原理进行的目标重 建技术相比，本发明考虑了大气折射影响，定位的目标位置更接近于目标 实际位置，满足同步卫星导弹预警的要求，定位精度准确、可靠性高。

本发明所提供的一种大气折射修正状态下目标重建方法，包括确定目 标和地心的直线方程，根据相机C₁确定目标实际位置P₁'，根据相机C₂确定 目标实际位置P₂'，确定目标实际位置P'步骤，本发明包括如下步骤：

S1、根据两相机位置、由两相机发出的初始光线方向确定过目标和地 心的直线方程，包括以下子步骤：

S1.1、已知相机C₁的经度Lon₁、纬度Lat₁、地心高度H₁，相机C₂的经度 Lon₂、纬度Lat₂、地心高度H₂，将两相机C₁、C₂的经纬度、高度坐标分别 转换为地心坐标系下的坐标C₁＝(X₁,Y₁,Z₁)'，C₂＝(X₂,Y₂,Z₂)'；具体转换公式 如下：

X_i＝H_i·cos(Lat_i·π/180)·cos(Lon_i·π/180)

Y_i＝H_i·cos(Lat_i·π/180)·sin(Lon_i·π/180)

Z_i＝H_i·sin(Lat_i·π/180)...i＝1,2

S1.2、根据目标在相机C₁成像的像面位置p₁＝(u₁,v₁)'，相机光轴指向 相机在地心坐标系下的坐标C₁＝(X₁,Y₁,Z₁)'及相机内部 参数，求得由相机C₁发出并最终到达目标的光线在地心坐标系下的初始方 向向量具体包括以下步骤：

S1.2.1、对相机C₁，根据C₁在地心坐标系下的坐标C₁＝(X₁,Y₁,Z₁)'，相机 光轴的指向求得相机坐标系到地心坐标系的转换矩 阵M_C1。

S1.2.2、根据目标在相机C₁成像的像面位置p₁＝(u₁,v₁)'，相机C₁内部参 数：相机所成图像的行数row、列数col，相机的像源瞬时视场角λ。计算由 相机C₁发出并最终到达目标的光线在相机坐标系下的初始方向向量

S1.2.3、由相机坐标系到地心坐标系的转换矩阵M_C1，由相机C₁发出并 最终到达目标的光线在相机坐标系下的初始方向向量求得由相机C₁发 出并最终到达目标的光线在地心坐标系下的初始方向向量即

S1.3、同步骤S1.2的方法，根据目标在相机C₂成像的像面位置 p₂＝(u₂,v₂)'，相机光轴指向相机在地心坐标系下的 坐标C₂＝(X₂,Y₂,Z₂)'及相机内部参数，求得由相机C₂发出并最终到达目标的 光线初始方向在地心坐标系下的向量

S1.4、因为大气在水平方向视为均匀分布，大气折射只发生在垂直方向 上，所以，由初始方向向量地心到相机C₁的向量所组成的平面，即 为相机C₁和目标实际位置P'所在的过地心的平面。根据方向向量相机C₁坐标C₁＝(X₁,Y₁,Z₁)'、地心O坐标O＝(0,0,0)'，求得相机C₁、目标实际位置P'以 及地心所组成的平面OC₁P'的法向量

S1.5、同步骤S1.4的方法，根据相机C₂、方向向量地心O的坐标， 求得相机C₂，目标实际位置P'以及地心所组成的平面OC₂P'的法向量${\overrightarrow{m}}_2={(m_{2x},m_{2y},m_{2z})}^{\prime }.$

S1.6、由成像原理可知，目标实际位置P'既在平面OC₁P'上，又在平面 OC₂P'上，因此目标位于两平面的交线上。根据两平面求得交线OP'的直线 方程，该方程即为目标和地心的直线方程，交线OP'的方向向量可由两 平面法向量的叉积求得，即${\overrightarrow{d}}_{{\mathrm{OP}}^{\prime }}={\overrightarrow{m}}_1\times {\overrightarrow{m}}_2={(d_x,d_y,d_z)}^{\prime }.$

S2、根据相机C₁、过目标和地心的直线方程，确定目标实际位置P₁‘；， 包括以下子步骤：

S2.1、对于相机C₁所在的平面OC₁P'，由直线OP'方程、相机C₁以及地心 O，求得目标实际位置P'与相机C₁间的地心夹角θ₁， 将地心夹角θ₁分为num₁(100≤num₁≤300)个相等的小夹角，对应的目标和相机 间的大气被划分为num₁层，每层地心夹角大小为dθ，dθ＝θ₁/num₁。因为每层 地心夹角很小，这使得每层大气的厚度很薄，所以可以取每层入射点位置 对应的大气折射率作为该层大气的折射率。

S2.2、求地心坐标系与平面坐标系之间的转换矩阵。具体包括如下步骤： 由相机C₁在地心坐标系下坐标C₁＝(X₁,Y₁,Z₁)'，交线OP'的方向向量叉积求 得相机、目标、地心所组成的平面坐标系Y_p轴的方向向量，即 X_p轴的方向向量由Y_p轴的方向向量与地心坐标 系下的Z轴的方向向量的叉积得到，即 $O{\overrightarrow{X}}_p={\overrightarrow{OY}}_p\times \overrightarrow{OZ}={(x_y,y_y,z_y)}^{\prime }\times {(0,0,1,)}^{\prime }={(x_x,y_x,z_x)}^{\prime };$根据右手螺旋定则，由X_p轴 的方向向量与Y_p轴方向向量的叉积可求得Z_p轴的方向向量， 从而求得的地心坐标系到平面坐标系的转换矩 阵M_earth-plan：

$M_{earth-plan}=\left(\begin{array}{cccc}{x}_x& y_x& z_x& 0\\ x_y& y_y& z_y& 0\\ x_z& y_z& z_z& 0\\ 0& 0& 0& 1\end{array}\right)$

平面坐标系到地心坐标系的转换矩阵M_plan-earth为：

$M_{plan-earth}=M_{earth-plan}^{\prime }=\left(\begin{array}{cccc}{x}_x& x_y& x_z& 0\\ y_x& y_y& y_z& 0\\ z_x& z_y& z_z& 0\\ 0& 0& 0& 1\end{array}\right)$

S2.3、根据地心坐标系到平面坐标系的转换矩阵M_earth-plan，将相机C₁在 地心坐标系下的坐标C₁＝(X₁,Y₁,Z₁)'转换为平面坐标系下的坐标 C_plan1＝(X_plan(1),Y_plan(1),Z_plan(1))'，转换公式为$\left(\begin{array}{c}{X}_{plan\left(1\right)}\\ Y_{plan\left(1\right)}\\ Z_{plan\left(1\right)}\\ 1\end{array}\right)=M_{earth-plan}.\left(\begin{array}{c}{X}_1\\ Y_1\\ Z_1\\ 1\end{array}\right).$因为平面 OC₁P'中任意点坐标的Y_p轴坐标始终为0，为了方便计算，在平面坐标系中， 将平面OC₁P'内点的坐标视为二维平面坐标的形式，即C_plan1＝(X_plan(1),Z_plan(1))'。

S2.4、根据地心坐标系到平面坐标系的转换矩阵M_earth-plan，将地心坐标 系下光线初始方向向量转换为平面坐标系中的向量

S2.5、在平面OC₁P'中，以相机C₁所在的位置为光线由第0层大气射向 第1层大气时，在第1层的入射点位置。根据光线的初始方向向量相 机坐标C_plan1，计算光线由第0层大气射向第1层大气时，第1层出射光线 的方向向量，以及第2层入射光线的入射点位置。具体包括以下步骤：

S2.5.1、相机C₁所在位置为光线由第0层射向第1层时的入射点位置T₁； 光线初始方向向量即为光线由第0层大气射向第1层大气时，第1 层的入射光线方向向量根据相机C₁的经纬度坐 标将相机C₁对应的地面位置的海拔高度、地面温度、地面压强、地面相对 湿度、季节、以及C₁处地心高度H₁参数，代入大气折射率Hopfiled模型， 计算C₁处的大气折射率N₁，根据大气折射指数公式n＝N×10^-6+1，计算第1 层大气折射指数n₁。考虑实际情况，第0层大气折射指数n₀和第1层大气折 射指数n₁取相同值，n₀＝n₁。光线由第0层大气射向第1层大气时不发生折 射，第1层的出射光线方向向量等于第1层的入射光线方向向量所以求 得第1层出射光线方向的单位向量

S2.5.2、由第1层的入射光线方向向量和入射点位置T₁处的法线向量 ${\overrightarrow{OT}}_1={(X_{plan1},Z_{plan1})}^{\prime },$根据公式${\alpha }_1=arccos\left(\frac{{\overrightarrow{OT}}_1·\overrightarrow{v_1}}{\left|{\overrightarrow{OT}}_1\right|·\left|\overrightarrow{v_1}\right|}\right),$计算第1层入射光线的入 射角α₁，第1层出射光线的出射角β₁等于入射角α₁。

S2.5.3、在由地心O、第1层的入射点位置T₁、第2层的入射点位置T₂组 成的三角形OT₁T₂中，由正弦定理，计算光线在第1层中的路径长度L₁，具 体公式为：

$L_1=\frac{H_1·sin\left(d\theta \right)}{sin({\beta }_1-d\theta )}$

S2.5.4、进而求得第1层出射光线方向向量该向量即为第2 层入射光线方向向量根据向量间的关系求得第2层入射点位 置T₂，$T_2=T_1+{\overrightarrow{u}}_1={(X_{plan\left(2\right)},Z_{plan\left(2\right)})}^{\prime }.$

S2.6、当光由第i(i＝1,2,3,...,num₁-1)层大气射向第i+1层大气时，计算第 i+1层出射光线方向向量以及第i+2层光线入射点位置，具体包括：

S2.6.1、当光由第i(i＝1,2,3,...,num₁-1)层大气射向第i+1层大气时，第i+1 层入射点位置T_i+1＝(X_plan(i+1),Z_plan(i+1))'、第i层大气折射指数n_i以及每层地心夹 角dθ已知，由第i+1层入射点位置T_i+1处的法线向量与入射光线的方向 向量根据公式求得第i+1层入射光线的入射角 α_i+1。利用平面坐标系到地心坐标系的转换矩阵M_plan-earth，将第i+1层入射点 位置T_i+1在平面坐标系下的坐标转换到地心坐标系下的坐标 Te_i+1＝(X_i+1,Y_i+1,Z_i+1)'，转换公式为$\left(\begin{array}{c}{X}_{i+1}\\ Y_{i+1}\\ Z_{i+1}\\ 1\end{array}\right)=M_{plan-earth}.\left(\begin{array}{c}{X}_{plan(i+1)}\\ 0\\ Z_{plan(i+1)}\\ 1\end{array}\right);$根据地心坐标到 经纬度、高度的转换公式，求出并记录第i+1层入射点位置的经纬度、高度 坐标Tl_i+1＝(Lon_i+1,Lat_i+1,H_i+1)'。地心坐标Te_i+1＝(X_i+1,Y_i+1,Z_i+1)'到经纬度、高度 Tl_i+1＝(Lon_i+1,Lat_i+1,H_i+1)'的转换公式为：

$H_{i+1}=\sqrt{X_{k+1}^2+Y_{i+1}^2+Z_{i+1}^2}$

${\mathrm{Lat}}_{i+1}=arcsin\left(\frac{Z_{i+1}}{H_{i+1}}\right)\times \frac{180}{\pi }$

${\mathrm{Lon}}_{i+1}=\left(\begin{array}{cc}\arcsin \left(\frac{Y_{i+1}}{H_{i+1}\times \cos \left({\mathrm{Lat}}_{i+1}\right)}\right)\times \frac{180}{\pi }& X_{i+1}>0,Y_{i+1}\in R\\ 180-\arcsin \left(\frac{Y_{i+1}}{H_{i+1}\times \cos \left({\mathrm{Lat}}_{i+1}\right)}\right)\times \frac{180}{\pi }& X_{i+1}\le 0,Y_{i+1}>0\\ -180-\arcsin \left(\frac{Y_{i+1}}{H_{i+1}\times \cos \left({\mathrm{Lat}}_{i+1}\right)}\right)\times \frac{180}{\pi }& X_{i+1}\le 0,Y_{i+1}\le 0\end{array}\right)$

S2.6.2、根据第i+1层入射点位置的经纬度坐标，将入射点位置T_i+1对应 地面处的海拔高度、地面温度、地面压强。地面相对湿度、入射点位置T_i+1的 地心高度H_i+1、季节参数代入到Hopfield模型中，计算第i+1层大气折射率 N_i+1，由大气折射指数公式n＝N×10^-6+1，计算第i+1层大气折射指数n_i+1。

S2.6.3、判断第i+1层大气折射指数n_i+1与第i层大气折射指数n_i之间的大 小关系，如果n_i+1小于n_i，此时可能会发生全反射，则执行步骤S2.6.4；否 则直接执行步骤S2.6.5。

S2.6.4、根据全反射临界角公式计算临界角θ。判断入 射角α_i+1与临界角θ大小关系，如果入射角α_i+1小于临界角θ，则不会发生全 反射，执行步骤S2.6.5，并且每次执行步骤S2.6.5后都要重新计算临界角并 判断与入射角的大小关系；如果入射角α_i+1大于等于临界角θ，则光线在该 层发生全反射，光线的传输方向发生变化，反射角ε_i+1等于入射角α_i+1。此时 可以将第i+1层出射角看作反射角ε_i+1的补角，求得第i+1层出射光线的出射 角β_i+1，β_i+1＝π-ε_i+1，继续执行步骤S2.6.6，此后每一层折射时不需要再判 断入射角与临界角的大小关系。

S2.6.5、根据大气球面分层的Snell定律，计算第i+1层初始光线的出射 角β_i+1，${\beta }_{i+1}=\frac{n_i}{n_{i+1}}{\mathrm{sin\alpha }}_{i+1}.$

大气球面分层的Snell定律为：当光线由大气折射指数为n₁的大气层， 射向大气折射指数为n₂的大气层时，满足下面的关系：

n₁sinφ₁＝n₂sinφ₂

式中φ₁、φ₂分别为入射角和出射角，它们分别定义为入射方向和出射方向与 界面法线的夹角，在地球大气层中，界面法线为地心与入射点的连线。

S2.6.6、由第i+1层入射点位置处的法线向量出 射光线的出射角β_i+1，利用向量内积，计算第i+1层出射光线的单位向量 具体求解方法为：

$\left(\begin{array}{cc}{\mathrm{cos\beta }}_{i+1}=\overrightarrow{u_e}·O{\overrightarrow{T}}_{i+1}=u_x*X_{plan(i+1)}+u_z*Z_{plan(i+1)}& \left(1\right)\\ u_x*u_x+u_z*u_z=1& \left(2\right)\end{array}\right)$

S2.6.7、在由地心O、第i+1层的入射点位置T_i+1、第i+2层的入射点位 置T_i+2组成的三角形OT_i+1T_i+2中，根据正弦定理，由公式计 算光线在第i+1层的路径长度L_i+1，进而求得第i+1层出射光线的方向向量 该方向向量即为第i+2层入射光线的方向向量即根 据向量间的关系计算第i+2层光线入射点位置T_i+2， $T_{i+2}=T_{i+1}+\overrightarrow{L_{i+1}{u}_e}{(X_{plan(i+2)},Z_{plan(i+2)})}^{\prime }.$

S2.7、当总层数等于num₁-1，此时对应的下一层光线入射点位置即 为通过相机C₁求得的目标位置P₁'。通过平面坐标系到地心坐标系的转换矩 阵M_plan-earth将平面坐标系下的坐标转换到地心坐标系下的坐标根 据地心坐标到经纬度、高度的转换公式，求得处的经纬度、高度坐标并 记录，此坐标即为通过相机C₁求得的目标位置P₁'。

S3、根据相机C₂、过目标和地心的直线方程，确定目标实际位置P₂‘；

对于相机C₂所在的平面OC₂P'中，求目标实际位置P₂'方法与步骤S2相 同。最终根据相机C₂求得目标实际位置为P₂'。

S4、根据两相机C₁、C₂确定的目标实际位置P₁'、P₂'，最终得到准确的 目标实际位置P'。

为了使目标定位更准确，最终目标实际位置P'取两目标位置P₁'、P₂'的 平均值，即P'＝(P₁'+P₂')/2。

本发明根据目标在相机成像的像面位置，通过三维目标重建的方法获 得目标的实际位置。相比于现有的目标重建方法，本发明考虑了光线传输 过程中大气折射的影响，并针对不同地区、不同季节、不同高度时大气的 折射率不同的问题，在计算大气折射率时，按照不同地区、不同季节、不 同高度的情况考虑；并在考虑大气折射传播时，考虑了发生全反射时的光 线传输路径。本发明重建的目标位置更接近于目标真实位置，提高了目标 定位的精确度，可更有效地应用于导弹预警，提高预警防御的准确性。

附图说明

图1本发明大气折射修正状态下目标重建方法的流程示意图；

图2大气折射修正状态下目标重建光线传输图；

图3地心坐标系与平面坐标系之间的关系图；

图4相机、目标、地心所在的平面内，光线由第0层射向第1层传输 路径图；

图5相机、目标、地心所在的平面内，光在每层大气中的传输路径图；

图6发生全反射时，光线传输路径图；

图7相机所在无人机与目标间的直线距离、两无人机间的夹角、无人 机的姿态、目标距离地面的高度、相机光轴指向一定时，不同的像面误差 与重建误差间的关系图；

图8相机所在无人机与目标间的直线距离、两无人机间的夹角、无人 机的姿态、目标距离地面的高度、目标在像面成像位置一定时，不同的相 机光轴指向误差与重建误差间的关系图；

图9相机所在无人机与目标间的直线距离、两无人机间的夹角、目标 距离地面的高度、相机光轴指向、目标在像面成像位置一定时，不同的无 人机姿态误差与重建误差间的关系图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图 及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体 实施例仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本 发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以 相互组合。

由目标发出的光线，经过大气传输，被相机接收并在像面上成像。根 据光路可逆原则，由目标发出并最终到达相机的光线传输路径与由相机发 出并最终到达目标的光线传输路径是相同的。所以在目标重建时，我们将 光线传输路径视为由相机出发并最终到达目标位置。根据目标的像面点位 置，可以得到最终入射光线方向，该方向的反方向即为从相机出发并最终 到达目标的光线初始方向，大气折射修正状态下目标重建的光线传输图， 如图2所示。本发明流程如图1所示，本发明提供一种折射修正状态下目 标重建方法，包括以下步骤：

S1、根据两相机位置、由两相机发出的初始光线方向确定过目标和地 心的直线方程；

S2、根据相机C₁、过目标和地心的直线方程，确定目标实际位置P₁‘；

S3、根据相机C₂、过目标和地心的直线方程，确定目标实际位置P₂‘；

S4、根据两相机C₁、C₂确定的目标实际位置P₁‘、P₂‘，最终得到准确的 目标实际位置P'。

其中，上述步骤S1具体包括以下子步骤：

S1.1、已知相机C₁的经度Lon₁、纬度Lat₁、地心高度H₁，相机C₂的经度 Lon₂、纬度Lat₂、地心高度H₂，将两相机C₁、C₂的经纬度、高度坐标分别 转换为地心坐标系下的坐标C₁＝(X₁,Y₁,Z₁)'，C₂＝(X₂,Y₂,Z₂)'。具体转换公式 如下：

X_i＝H_i·cos(Lat_i·π/180)·cos(Lon_i·π/180)

Y_i＝H_i·cos(Lat_i·π/180)·sin(Lon_i·π/180)

Z_i＝H_i·sin(Lat_i·π/180)...i＝1,2

S1.2、根据目标在相机C₁成像的像面位置p₁＝(u₁,v₁)'，相机光轴指向 相机在地心坐标系下的坐标C₁＝(X₁,Y₁,Z₁)'及相机内部 参数，求得由相机C₁发出并最终到达目标的光线在地心坐标系下的初始方 向向量具体包括以下步骤：

S1.2.1、对相机C₁，根据C₁在地心坐标系下的坐标C₁＝(X₁,Y₁,Z₁)'，相机 光轴的指向计算相机坐标系到地心坐标系的转换矩 阵M_C1。

下面先介绍相机坐标系的定义，再介绍相机坐标系到地心坐标系的转 换矩阵的求取方法。

(1)相机坐标系的定义

定义以相机所在位置为相机坐标系的原点O_C1，相机光轴的指向为 相机坐标系的Z_C1轴，X_C1轴垂直于相机和地心的连线OC₁与相机光轴指向 所组成的平面，Y_C1轴由右手螺旋定则确定。

(2)相机坐标系到地心坐标系的转换矩阵

相机坐标系的Z_C1轴与相机光轴指向重合，可以求得Z_C1轴的方向向量， 因为X_C1轴垂直于相机和地心的连线OC₁与相 机光轴指向所组成的平面，所以X_C1轴方向向量可由地心O到相机 C₁的向量与相机光轴的指向的叉积求得，即 根据右手螺旋定则，Y_C1轴由X_C1轴和Z_C1轴的叉 积确定，所以可求得Y_C1轴的方向向量${\overrightarrow{OY}}_{C1}={\overrightarrow{OX}}_{C1}\times {\overrightarrow{OZ}}_{C1}={({\mu }_x,{\mu }_y,{\mu }_z)}^{\prime }.$从而由相机坐标系到地心坐标系的转换矩阵为：

$M_{C1}=\left(\begin{array}{cccc}{\eta }_x& {\mu }_x& {\mathrm{vec}}_{1x}& 0\\ {\eta }_y& {\mu }_y& {\mathrm{vec}}_{1y}& 0\\ {\eta }_z& {\mu }_z& {\mathrm{vec}}_{1z}& 0\\ 0& 0& 0& 1\end{array}\right)$

S1.2.2、根据目标在相机C₁成像的像面位置p₁＝(u₁,v₁)'，相机C₁内部参 数：相机所成图像的行数row、列数col，相机的像源瞬时视场角λ；计算由 相机C₁发出并最终到达目标的光线在相机坐标系下的初始方向向量具 体计算公式如下：

γ₁＝(u₁-row/2)·λ

γ₂＝(v₁-col/2)·λ

${\overrightarrow{V}}_{C1}=\left(\begin{array}{cccc}{\mathrm{cos\gamma }}_1& 0& {\mathrm{sin\gamma }}_1& 0\\ 0& 1& 0& 0\\ -{\mathrm{sin\gamma }}_1& 0& {\mathrm{cos\gamma }}_1& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& {\mathrm{cos\gamma }}_2& -{\mathrm{sin\gamma }}_2& 0\\ 0& {\mathrm{sin\gamma }}_2& {\mathrm{cos\gamma }}_2& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{c}0\\ 0\\ 1\\ 1\end{array}\right)$

S1.2.3、由相机坐标系到地心坐标系的转换矩阵M_C1，由相机C₁发出并 最终到达目标的光线在相机坐标系下的初始方向向量计算由相机C₁发 出并最终到达目标的光线在地心坐标系下的初始方向向量即

S1.3、同步骤S1.2的方法，根据目标在相机C₂成像的像面位置 p₂＝(u₂,v₂)'，相机光轴指向相机在地心坐标系下的 坐标C₂＝(X₂,Y₂,Z₂)'及相机内部参数，计算由相机C₂发出并最终到达目标的 光线在地心坐标系下的初始方向向量

S1.4、因为大气在水平方向视为均匀分布，大气折射只发生在垂直方向 上，所以，由初始方向向量地心到相机C₁的向量所组成的平面，即 为相机C₁和目标实际位置P'所在的过地心的平面。根据方向向量相机C₁坐标C₁＝(X₁,Y₁,Z₁)'、地心O坐标O＝(0,0,0)'，计算相机C₁、目标实际位置P'所 在过地心的平面OC₁P'的法向量求平面OC₁P'的法向量具体 方法如下：

在平面OC₁P'中，由地心到相机的向量和光线初始方向向量的叉 积可以得到平面OC₁P'的法向量

S1.5、同步骤S1.4的方法，根据相机C₂、方向向量地心O的坐标， 求得相机C₂，目标实际位置P'以及地心所组成的平面OC₂P'的法向量 ${\overrightarrow{m}}_2={(m_{2x},m_{2y},m_{2z})}^{\prime }.$

S1.6、由成像原理可知，目标实际位置P'既在平面OC₁P'上，又在平面 OC₂P'上，因此目标位于两平面的交线OP'上。根据两平面求得交线OP'的直 线方程，该方程即为过目标和地心的直线方程，交线OP'的方向向量可 由两平面的法向量的叉积求得，即

步骤S2具体包括以下子步骤：

S2.1、对于相机C₁所在的平面OC₁P'，由直线OP'方程、相机C₁以及地心 O，计算目标实际位置P'与相机C₁间的地心夹角θ₁， 将地心夹角θ₁分为num₁(100≤num₁≤300)个相等的小夹角，对应的目标和相机 间的大气被划分为num₁层，每层地心夹角大小为dθ，dθ＝θ₁/num₁。因为每层 地心夹角很小，这使得每层大气的厚度很薄，所以可以取每层入射点位置 对应的大气折射率作为该层大气的折射率。

S2.2、求地心坐标系与平面坐标系之间的转换矩阵。

下面先介绍平面坐标系的定义，以及平面坐标系与地心坐标系之间的 转换矩阵：

(1)相机、目标以及地心所组成的平面坐标系定义如下：

如图3所示，设相机C在地心坐标系中坐标(x_c,y_c,z_c)'，目标P在地心坐 标系中坐标(x_p,y_p,z_p)'，地心O坐标(0,0,0)'。定义平面坐标系O-X_pY_pZ_p的原 点与地心直角坐标系的原点重合，设平面OCP与平面OXY的交线OX_p为平 面坐标系的X_p轴；Y_p轴垂直平面OCP；则由右手螺旋定则确定Z_p轴。

(2)地心坐标系与平面坐标系之间的转换矩阵

如图3所示，地心坐标下，向量OC坐标(x_c,y_c,z_c)'，向量OP坐标 (x_p,y_p,z_p)'，平面坐标系Y_p轴垂直平面OCP，所以向量OC与向量OP叉积所 得向量即为Y_p轴的方向向量因为X_p轴为平面 OCP与平面OXY的交线，所以地心坐标系中的Z轴垂直于X_p轴，由Y_p轴和Z 轴的方向向量的的叉积可以得到X_p轴的方向向量根据右手螺旋定则，X_p轴与Y_p轴的叉积即为Z_p轴， 所以可求得轴的方向向量，${\overrightarrow{OZ}}_p={\overrightarrow{OX}}_p\times {\overrightarrow{OY}}_p={(x_z,y_z,z_z)}^{\prime }.$

从而地心坐标系O-XYZ到平面坐标系O-X_pY_pZ_p的转换矩阵为：

$M_{earth-plan}=\left(\begin{array}{cccc}{x}_x& y_x& z_x& 0\\ x_y& y_y& z_y& 0\\ x_z& y_z& z_z& 0\\ 0& 0& 0& 1\end{array}\right)$

平面坐标系O-X_pY_pZ_p到地心坐标系O-XYZ的转换矩阵为：

$M_{plan-earth}=M_{earth-plan}^{\prime }=\left(\begin{array}{cccc}{x}_x& x_y& x_z& 0\\ y_x& y_y& y_z& 0\\ z_x& z_y& z_z& 0\\ 0& 0& 0& 1\end{array}\right)$

利用地心坐标系与平面坐标系的转换矩阵，可将任一点地心坐标系下 的坐标(X,Y,Z)'与平面坐标系下的坐标(X_p,Y_p,Z_p)'相互转换，其转化关系如 下：

$\left(\begin{array}{c}{X}_P\\ Y_P\\ Z_P\\ 1\end{array}\right)=M_{earth-plan}.\left(\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right)$

$\left(\begin{array}{c}X\\ Y\\ Z\\ 1\end{array}\right)=M_{plan-earth}.\left(\begin{array}{c}{X}_P\\ Y_P\\ Z_P\\ 1\end{array}\right)$

对于平面OCP中的任一点，在平面坐标系中，其Y_p轴坐标始终为0，即 平面中一点在地心坐标系下坐标(x_e,y_e,z_e)'，经转换矩阵得到平面坐标系下 坐标(x_p,0,z_p)'。且已知光线只在垂直方向上发生折射，在水平方向不发生偏 折，所以由相机发出到达目标的光线传输路径上每一点，在平面坐标系中Y_p轴坐标始终为0。

对于相机C₁所在的平面OC₁P'中，求地心坐标系与平面坐标系之间的转 换矩阵具体方法如下：由相机C₁在地心坐标系下坐标C₁＝(X₁,Y₁,Z₁)'和交线 OP'的方向向量叉积求得相机、目标、地心所组成的平面坐标系的Y_p轴 的方向向量，即X_p轴的方向向量由Y_p轴的方向 向量与地心坐标系下的Z轴的方向向量叉积得到，即 $O{\overrightarrow{X}}_p={\overrightarrow{OY}}_p\times \overrightarrow{OZ}={(x_y,y_y,z_y)}^{\prime }\times {(0,0,1,)}^{\prime }={(x_x,y_x,z_x)}^{\prime };$由X_p轴的方向向量与Y_p轴的 方向向量叉积可求得Z_p轴的方向向量，从而求 得的地心坐标系到平面坐标系的转换矩阵M_earth-plan：

$M_{earth-plan}=\left(\begin{array}{cccc}{x}_x& y_x& z_x& 0\\ x_y& y_y& z_y& 0\\ x_z& y_z& z_z& 0\\ 0& 0& 0& 1\end{array}\right)$

平面坐标系到地心坐标系的转换矩阵M_plan-earth为：

$M_{plan-earth}=M_{earth-plan}^{\prime }=\left(\begin{array}{cccc}{x}_x& x_y& x_z& 0\\ y_x& y_y& y_z& 0\\ z_x& z_y& z_z& 0\\ 0& 0& 0& 1\end{array}\right)$

S2.3、根据地心坐标系到平面坐标系的转换矩阵M_earth-plan，将相机C₁在 地心坐标系下的坐标C₁＝(X₁,Y₁,Z₁)'转换为平面坐标系下的坐标 C_plan1＝(X_plan(1),Y_plan(1),Z_plan(1))'，转换公式为$\left(\begin{array}{c}{X}_{plan\left(1\right)}\\ Y_{plan\left(1\right)}\\ Z_{plan\left(1\right)}\\ 1\end{array}\right)=M_{earth-plan}.\left(\begin{array}{c}{X}_1\\ Y_1\\ Z_1\\ 1\end{array}\right).$因为平面 OC₁P'中任意点坐标的Y_p轴坐标始终为0，为了方便计算，在平面OC₁P'坐标 系中，将平面内OC₁P'的点的坐标视为二维平面坐标的形式，即 C_plan1＝(X_plan(1),Z_plan(1))'。

S2.4、根据地心坐标系到平面坐标系的转换矩阵M_earth-plan，将光线在地 心坐标系下初始方向向量转换为平面坐标系中的向量

S2.5、在平面OC₁P'中，以相机C₁所在的位置为光线由第0层大气射向 第1层大气时，在第1层大气的入射点位置。根据光线的初始方向向量相机坐标C_plan1，计算光线由第0层大气射向第1层大气时，第1层出射光 线的方向向量，以及第2层入射光线的入射点位置。具体包括以下步骤：

S2.5.1、相机C₁所在位置为光线由第0层大气射向第1层大气时的入射 点位置T₁；光线初始方向向量即为光线由第0层大气进入第1层大气 时，第1层入射光线的方向向量根据相机C₁的经纬 度坐标将C₁对应地面位置的海拔高度、地面温度、地面压强、地面相对湿 度、季节参数、C₁处地心高度H₁，代入大气折射率Hopfiled模型，计算C₁处 的大气折射率N₁，根据大气折射指数公式n＝N×10^-6+1，计算第1层大气折 射指数n₁。考虑实际情况，第0层大气折射指数n₀和第1层大气折射指数n₁取相同值，n₀＝n₁。此时光由第0层射向第1层时不发生折射，如图4所示， 第1层出射光线的方向向量等于第1层入射光线的方向向量所以求得第 1层出射光线的单位向量

S2.5.2、由第1层入射光线的方向向量和入射点位置T₁处的法线向量 ${\overrightarrow{OT}}_1={(X_{plan1},Z_{plan1})}^{\prime },$根据公式${\alpha }_1=arccos\left(\frac{{\overrightarrow{OT}}_1·\overrightarrow{v_1}}{\left|{\overrightarrow{OT}}_1\right|·\left|\overrightarrow{v_1}\right|}\right),$计算第1层入射光线的入 射角α₁，该入射角等于第1层出射光线的出射角β₁。

S2.5.3、如图4所示，在三角形OT₁T₂中，由正弦定理，计算光线在第1 层中的路径长度L₁，具体公式为：

$L_1=\frac{H_1·sin\left(d\theta \right)}{sin({\beta }_1-d\theta )}$

S2.5.4、求得第1层出射光线的方向向量该方向向量即为第 2层入射光线的方向向量根据向量间的关系求得第2层入射 点位置T₂，$T_2=T_1+{\overrightarrow{u}}_1={(X_{plan\left(2\right)},Z_{plan\left(2\right)})}^{\prime }.$

大气折射率的Hopfield模型具体原理如下：

利用大气折射率的Hopfield模型，可以求取任意海拔高度h的大气折射 率N(h)，Hopfield模型考虑了大气含有水汽的情况，将大气折射率分为干、 湿两项，模型的具体形式为：

$\left(\begin{array}{cc}N\left(h\right)=N_d\left(h\right)+N_w\left(h\right)& \\ N_i\left(h\right)=N_{i0}{\left(\frac{H_i-h}{H_i-h_0}\right)}^4& \begin{array}{cc}i=d,w& h\le H_i\\ i=d,w& h\ge H_i\end{array}\\ N_i=0& \end{array}\right)$

式中N_d表示折射率干项，N_w表示折射率湿项；H_d表示干项等效高度， 即大气折射率中干项衰减为零的高度，H_d＝40.136+0.14872·t₀(km)；H_w表示 湿项等效高度，即大气折射率中湿项衰减为零的高度，一般取H_w＝11(km)； h为目标海拔高度；h₀为地面海拔高度。N_d0和N_w0分别表示测站地面大气折 射率的干项和湿项，N_d0＝77.6P₀/T₀，N_w0＝3.73×10⁵×e₀/T₀²，且T₀、P₀、e₀分 别表示地面温度(K),大气压(mbar)和水汽压(mbar)，T₀＝273+t₀， $e_0=f_0\times 6.11\times {10}^{(7.45t_0/(273.5+t_0\left)\right)},$f₀表示地面相对湿度。

S2.6、当光由第i(i＝1,2,3,...,num₁-1)层大气射向第i+1层大气时，计算第 i+1层出射光线的方向向量以及第i+2层入射光线的入射点位置，具体包括 如下子步骤：

S2.6.1、当光由第i(i＝1,2,3,...,num₁-1)层大气射向第i+1层大气时，如图 5所示，第i+1层的入射点位置T_i+1＝(X_plan(i+1),Z_plan(i+1))'、第i层的大气折射指数 n_i以及每层地心夹角dθ已知，由第i+1层入射点位置处的法线向量与第 i+1层入射光线的方向向量根据公式计算第i+1 层入射光线的入射角α_i+1。利用平面坐标系到地心坐标系的转换矩阵M_plan-earth， 将第i+1层的入射点位置T_i+1在平面坐标系下的坐标转换到地心坐标系下 Te_i+1＝(X_i+1,Y_i+1,Z_i+1)'，转换公式为$\left(\begin{array}{c}{X}_{i+1}\\ Y_{i+1}\\ Z_{i+1}\\ 1\end{array}\right)=M_{plan-earth}.\left(\begin{array}{c}{X}_{plan(i+1)}\\ 0\\ Z_{plan(i+1)}\\ 1\end{array}\right);$根据地心坐标到 经纬度、高度的转换公式，计算第i+1层入射点位置的经纬度、高度坐标 Tl_i+1＝(Lon_i+1,Lat_i+1,H_i+1)'，并记录第i+1层入射点位置的经纬度、高度坐标。 地心坐标Te_i+1＝(X_i+1,Y_i+1,Z_i+1)'到经纬度、高度Tl_i+1＝(Lon_i+1,Lat_i+1,H_i+1)'的转换公 式为：

$H_{i+1}=\sqrt{X_{i+1}^2+Y_{i+1}^2+Z_{i+1}^2}$

${\mathrm{Lat}}_{i+1}=arcsin\left(\frac{Z_{i+1}}{H_{i+1}}\right)\times \frac{180}{\pi }$

${\mathrm{Lon}}_{i+1}=\left(\begin{array}{cc}\arcsin \left(\frac{Y_{i+1}}{H_{i+1}\times \cos \left({\mathrm{Lat}}_{i+1}\right)}\right)\times \frac{180}{\pi }& X_{i+1}>0,Y_{i+1}\in R\\ 180-\arcsin \left(\frac{Y_{i+1}}{H_{i+1}\times \cos \left({\mathrm{Lat}}_{i+1}\right)}\right)\times \frac{180}{\pi }& X_{i+1}\le 0,Y_{i+1}>0\\ -180-\arcsin \left(\frac{Y_{i+1}}{H_{i+1}\times \cos \left({\mathrm{Lat}}_{i+1}\right)}\right)\times \frac{180}{\pi }& X_{i+1}\le 0,Y_{i+1}\le 0\end{array}\right)$

S2.6.2、根据第i+1层入射点位置的经纬度、高度坐标，将入射点T_i+1对 应地面处的海拔高度、地面温度、地面压强。地面相对湿度、入射点T_i+1处 的地心高度H_i+1、季节参数代入到Hopfield模型中，计算第i+1层大气折射 率N_i+1，由大气折射指数公式n＝N×10^-6+1，计算第i+1层大气折射指数n_i+1。

S2.6.3、判断第i+1层大气折射指数n_i+1与第i层大气折射指数n_i之间的大 小关系，如果n_i+1小于n_i，此时可能会发生全反射，则执行步骤S2.6.4；否 则直接执行步骤S2.6.5。

S2.6.4、根据全反射临界角公式求出临界角θ。判断入 射角α_i+1与临界角θ大小关系，如果入射角α_i+1小于临界角θ，则不会发生全 反射，执行步骤S2.6.5，并且每次执行步骤S2.6.5后都要重新计算临界角并 判断与入射角的大小关系；如果入射角α_i+1大于等于临界角θ，则光线在该 层发生全反射，光线的传输方向发生变化，反射角ε_i+1等于入射角α_i+1，如图 6所示，此时可以将第i+1层出射角看作反射角的补角，求得第i+1层出射 光线的出射角β_i+1，β_i+1＝π-ε_i+1，继续执行步骤S2.6.6，此后每一层折射时 不需要再判断入射角与临界角的大小关系。

S2.6.5、根据大气球面分层的Snell定律，计算第i+1层出射光线的出射 角β_i+1，${\beta }_{i+1}=\frac{n_i}{n_{i+1}}{\mathrm{sin\alpha }}_{i+1}.$

大气球面分层的Snell定律为：当光线由大气折射指数为n₁的大气层， 射向大气折射指数为n₂的大气层时，满足下面的关系：

n₁sinφ₁＝n₂sinφ₂

式中φ₁、φ₂分别为入射角和出射角，它们分别定义为入射方向和出射方向与 界面法线的夹角，在地球大气层中，界面法线为地心与入射点的连线。

S2.6.6、由第i+1层入射点位置T_i+1处的法线向量出射角β_i+1，利用向量内积，计算第i+1层出射光线的单位向量具体求解方法为：

$\left(\begin{array}{cc}{\mathrm{cos\beta }}_{i+1}=\overrightarrow{u_e}·O{\overrightarrow{T}}_{i+1}=u_x*X_{plan(i+1)}+u_z*Z_{plan(i+1)}& \left(1\right)\\ u_x*u_x+u_z*u_z=1& \left(2\right)\end{array}\right)$

S2.6.7、如图5所示，在三角形OT_i+1T_i+2中，根据正弦定理，由公式 求得光线在第i+1层中的路径长度L_i+1，进而求得第i+1层 出射光线的方向向量该方向向量即为第i+2层入射光线的方向向量 即${\overrightarrow{v}}_{i+2}=\overrightarrow{L_{i+1}{u}_e}.$根据向量间的关系$T_{i+2}-T_{i+1}=\overrightarrow{L_{i+1}{u}_e},$计算第i+2层光线入 射点T_i+2，$T_{i+2}=T_{i+1}+\overrightarrow{L_{i+1}{u}_e}={(X_{plan(i+2)},Z_{plan(i+2)})}^{\prime }.$

S2.7、当总层数等于num₁-1时，此时对应的下一层光线入射点位置即为通过相机C₁求得的目标位置P₁'。通过平面坐标系到地心坐标系的转换 矩阵M_plan-earth将在平面坐标系下的坐标转换到地心坐标系下的坐标根据地心坐标到经纬度、高度的转换公式，求得处的经纬度、高度坐标 并记录，此坐标即为通过相机C₁求得的目标位置P₁'。

步骤S3根据相机C₂成像信息确定目标位置P₂'，求目标位置P₂'的方法与 步骤S2中根据相机C₁求目标位置P₁'的方法相同。最终根据相机C₂求得的目 标位置P₂'。

步骤S4具体包括以下子步骤：

通过相机C₁、C₂分别确定了目标的位置P₁'、P₂'，由于计算时会存在误 差，为了使目标定位更精确，最终目标实际位置P'取两平面求得的目标位 置P₁'、P₂'的平均值，即P'＝(P₁'+P₂')/2。

实验结果分析：

图7～9表示分别对目标在像面位置、相机光轴指向、相机所在的无人 机平台的姿态添加不同的误差时，利用本发明的方法进行目标位置重建， 得到重建的目标位置与目标实际位置间的误差，称该误差为重建误差。以 判断本专利方法在重建时有外界干扰情况下，对目标定位的精确度。

图7为对目标的像面成像位置分别添加3σ＝0,0.5,1,1.5个像素的像面 位置误差时，得到的重建误差结果；图8为对相机的光轴指向添加3σ从0.001° 到0.01°间隔为0.001°的指向误差，得到的重建误差结果；图9为对相机 所在的无人机平台的姿态添加3σ从0.005°到0.02°间隔0.005°的姿态误 差，得到的重建误差结果。

从上面的实验结果可以看出，本发明不仅重建的目标位置精确度高， 而且抗外界环境干扰能力强，在由外界误差干扰的情况下，目标重建误差 在1km以内，这对提高导弹的预警防御具有重要作用。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已， 并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等 同替换和改进等，均应包含在本发明的保护范围之内。