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

Abstract

本发明公开了一种大气折射修正状态下目标重建方法，属于图像处理的三维重建方法，解决现有目标重建方法未考虑大气折射对光线传播影响而导致的目标定位误差较大的问题，可用于远距离预警系统。本发明方法包括：(1)确定过目标和地心的直线方程；(2)根据相机C₁成像信息确定目标位置P₁ˊ；(3)根据相机C₂成像信息确定目标位置P₂ˊ；(4)得到准确的目标实际位置Pˊ。本发明考虑了大气折射的影响，并针对不同地区、不同季节、不同高度时大气的折射率不同的问题，计算大气折射率时，按照不同地区、不同季节、不同高度的情况考虑，本发明重建得到的目标位置更接近目标实际位置，具有更高的目标定位精度,为高精度预警信息的获取提供了保障,可提升预警系统的整体性能。

Description

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

技术领域

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

背景技术

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

发明内容

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

本发明所提供的一种大气折射修正状态下目标重建方法，包括确定目标和地心的直线方程，根据相机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'的法向量 ${\stackrel{\→}{m}}_2={(m_{2x},m_{2y},m_{2z})}^{\′}.$

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

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{\stackrel{\→}{X}}_p={\stackrel{\→}{OY}}_p\×\stackrel{\→}{OZ}={(x_y,y_y,z_y)}^{\′}\×{(0,0,1,)}^{\′}={(x_x,y_x,z_x)}^{\′};$ 根据右手螺旋定则，由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}^{\′}=\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₁处的法线向量 ${\stackrel{\→}{OT}}_1={(X_{plan1},Z_{plan1})}^{\′},$ 根据公式 ${\mathrm{\α}}_1=arccos\left(\frac{{\stackrel{\→}{OT}}_1\·\stackrel{\→}{v_1}}{\left|{\stackrel{\→}{OT}}_1\right|\·\left|\stackrel{\→}{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\mathrm{\θ}\right)}{sin({\mathrm{\β}}_1-d\mathrm{\θ})}$

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

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)\×\frac{180}{\mathrm{\π}}$

${\mathrm{Lon}}_{i+1}=\left\{\begin{array}{cc}\arcsin \left(\frac{Y_{i+1}}{H_{i+1}\×\cos \left({\mathrm{Lat}}_{i+1}\right)}\right)\×\frac{180}{\mathrm{\π}}& X_{i+1}>0,Y_{i+1}\∈R\\ 180-\arcsin \left(\frac{Y_{i+1}}{H_{i+1}\×\cos \left({\mathrm{Lat}}_{i+1}\right)}\right)\×\frac{180}{\mathrm{\π}}& X_{i+1}\≤0,Y_{i+1}>0\\ -180-\arcsin \left(\frac{Y_{i+1}}{H_{i+1}\×\cos \left({\mathrm{Lat}}_{i+1}\right)}\right)\×\frac{180}{\mathrm{\π}}& X_{i+1}\≤0,Y_{i+1}\≤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， ${\mathrm{\β}}_{i+1}=\frac{n_i}{n_{i+1}}{\mathrm{sin\α}}_{i+1}.$

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

n₁sinφ₁＝n₂sinφ₂

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

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

$\left\{\begin{array}{cc}{\mathrm{cos\β}}_{i+1}=\stackrel{\→}{u_e}\·O{\stackrel{\→}{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}+\stackrel{\→}{L_{i+1}{u}_e}{(X_{plan(i+2)},Z_{plan(i+2)})}^{\′}.$

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轴的方向向量 ${\stackrel{\→}{OY}}_{C1}={\stackrel{\→}{OX}}_{C1}\×{\stackrel{\→}{OZ}}_{C1}={({\mathrm{\μ}}_x,{\mathrm{\μ}}_y,{\mathrm{\μ}}_z)}^{\′}.$ 从而由相机坐标系到地心坐标系的转换矩阵为：

$M_{C1}=\left[\begin{array}{cccc}{\mathrm{\η}}_x& {\mathrm{\μ}}_x& {\mathrm{vec}}_{1x}& 0\\ {\mathrm{\η}}_y& {\mathrm{\μ}}_y& {\mathrm{vec}}_{1y}& 0\\ {\mathrm{\η}}_z& {\mathrm{\μ}}_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)·λ

${\stackrel{\→}{V}}_{C1}=\left[\begin{array}{cccc}{\mathrm{cos\γ}}_1& 0& {\mathrm{sin\γ}}_1& 0\\ 0& 1& 0& 0\\ -{\mathrm{sin\γ}}_1& 0& {\mathrm{cos\γ}}_1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& {\mathrm{cos\γ}}_2& -{\mathrm{sin\γ}}_2& 0\\ 0& {\mathrm{sin\γ}}_2& {\mathrm{cos\γ}}_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'的法向量 ${\stackrel{\→}{m}}_2={(m_{2x},m_{2y},m_{2z})}^{\′}.$

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轴，所以可求得轴的方向向量， ${\stackrel{\→}{OZ}}_p={\stackrel{\→}{OX}}_p\×{\stackrel{\→}{OY}}_p={(x_z,y_z,z_z)}^{\′}.$

从而地心坐标系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}^{\′}=\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{\stackrel{\→}{X}}_p={\stackrel{\→}{OY}}_p\×\stackrel{\→}{OZ}={(x_y,y_y,z_y)}^{\′}\×{(0,0,1,)}^{\′}={(x_x,y_x,z_x)}^{\′};$ 由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}^{\′}=\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₁处的法线向量 ${\stackrel{\→}{OT}}_1={(X_{plan1},Z_{plan1})}^{\′},$ 根据公式 ${\mathrm{\α}}_1=arccos\left(\frac{{\stackrel{\→}{OT}}_1\·\stackrel{\→}{v_1}}{\left|{\stackrel{\→}{OT}}_1\right|\·\left|\stackrel{\→}{v_1}\right|}\right),$ 计算第1层入射光线的入射角α₁，该入射角等于第1层出射光线的出射角β₁。

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

$L_1=\frac{H_1\·sin\left(d\mathrm{\θ}\right)}{sin({\mathrm{\β}}_1-d\mathrm{\θ})}$

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

大气折射率的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\≤H_i\\ i=d,w& h\≥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\×6.11\×{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)\×\frac{180}{\mathrm{\π}}$

${\mathrm{Lon}}_{i+1}=\left\{\begin{array}{cc}\arcsin \left(\frac{Y_{i+1}}{H_{i+1}\×\cos \left({\mathrm{Lat}}_{i+1}\right)}\right)\×\frac{180}{\mathrm{\π}}& X_{i+1}>0,Y_{i+1}\∈R\\ 180-\arcsin \left(\frac{Y_{i+1}}{H_{i+1}\×\cos \left({\mathrm{Lat}}_{i+1}\right)}\right)\×\frac{180}{\mathrm{\π}}& X_{i+1}\≤0,Y_{i+1}>0\\ -180-\arcsin \left(\frac{Y_{i+1}}{H_{i+1}\×\cos \left({\mathrm{Lat}}_{i+1}\right)}\right)\×\frac{180}{\mathrm{\π}}& X_{i+1}\≤0,Y_{i+1}\≤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， ${\mathrm{\β}}_{i+1}=\frac{n_i}{n_{i+1}}{\mathrm{sin\α}}_{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\β}}_{i+1}=\stackrel{\→}{u_e}\·O{\stackrel{\→}{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层入射光线的方向向量即 ${\stackrel{\→}{v}}_{i+2}=\stackrel{\→}{L_{i+1}{u}_e}.$ 根据向量间的关系 $T_{i+2}-T_{i+1}=\stackrel{\→}{L_{i+1}{u}_e},$ 计算第i+2层光线入射点T_i+2， $T_{i+2}=T_{i+1}+\stackrel{\→}{L_{i+1}{u}_e}={(X_{plan(i+2)},Z_{plan(i+2)})}^{\′}.$

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以内，这对提高导弹的预警防御具有重要作用。

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

Claims (10)

1.一种大气折射修正状态下的目标重建方法，其特征在于，所述方法包括步骤：

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

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

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

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

2.如权利要求1所述方法，其特征在于，所述步骤S1包括以下子步骤：

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

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

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

S1.4、根据初始方向向量相机C₁坐标C₁＝(X₁,Y₁,Z₁)'、地心O坐标O＝(0,0,0)'，计算相机C₁、目标实际位置P'以及地心所组成的平面OC₁P'的法向量：

S1.5、同步骤S1.4的方法，根据方向向量相机C₂坐标C₂＝(X₂,Y₂,Z₂)'、地心O的坐标O＝(0,0,0)'，计算相机C₂、目标实际位置P'以及地心所组成的平面OC₂P'的法向量：

S1.6、根据OC₁P'和OC₂P'两平面法向量，计算交线OP'的直线方程，该方程即为过目标和地心的直线方程。

3.如权利要求2所述方法，其特征在于，所述步骤S1.2具体包括以下步骤：

S1.2.1、根据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₁发出并最终到达目标的光线在地心坐标系下的初始方向向量 。

4.如权利要求1或2所述方法，其特征在于，所述步骤S2包括以下子步骤：

S2.1、由直线OP'方程、相机C₁以及地心O，计算目标实际位置P'与相机C₁间的地心夹角θ₁，并将地心夹角分为num₁(100≤num₁≤300)个相等的小夹角；

S2.2、计算地心坐标系与平面坐标系之间的转换矩阵M_earth-plan；

S2.3、根据地心坐标系到平面坐标系的转换矩阵M_earth-plan，将相机C₁在地心坐标系下的坐标C₁＝(X₁,Y₁,Z₁)'转换为平面坐标系下的坐标，并记为二维平面坐标的形式，即C_plan1＝(X_plan(1),Z_plan(1))'；

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

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

S2.6、当光由第i(i＝1,2,3,...,num₁-1)层大气射向第i+1层大气时，计算第i+1层大气的出射光线方向向量以及光线在第i+2层大气的入射点位置，其中num₁为步骤S2.1中将地心夹角等分的小夹角个数；

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

5.如权利要求2所述方法，其特征在于，所述步骤S1.6中的交线OP'的方向向量由OC₁P'和OC₂P'两平面的法向量的叉积求得，即

6.如权利要求4所述方法，其特征在于，所述步骤S2.5具体包括如下步骤：

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

S2.5.2、由第1层的入射方向向量和入射点位置T₁处的法线向量根据公式求得第1层入射光线的入射角α₁，该入射角等于第1层出射光线的出射角β₁；

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

S2.5.4、求得第1层大气的出射光线方向向量该方向向量即为第2层大气的入射光线方向向量根据向量间的关系计算第2层的入射点位置T₂，

7.如权利要求2所述方法，其特征在于，所述步骤S2.6中具体包括如下子步骤：

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)'，转换公式为根据地心坐标到经纬度、高度的转换公式，求得第i+1层入射点的经纬度、高度坐标Tl_i+1＝(Lon_i+1,Lat_i+1,H_i+1)'，并记录第i+1层的入射点位置处的经纬度、高度坐标；

S2.6.2、根据第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，继续执行步骤S2.6.6，此后每一层折射时不需要再判断入射角与临界角的大小关系；

S2.6.5、根据大气球面分层的Snell定律，计算第i+1层出射光线的出射角β_i+1，

S2.6.6、由第i+1层入射点T_i+1处的法线向量出射角β_i+1，利用向量内积，求得第i+1层出射光线的单位向量

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的坐标，

8.如权利要求7所述方法，其特征在于，所述步骤S2.6.6中的具体求解方法为：

9.如权利要求2所述方法，其特征在于，所述步骤S4具体为：通过相机C₁、C₂分别确定了目标的位置，最终目标实际位置P'取两平面求得的目标位置P₁'、P₂'的平均值，即P'＝(P₁'+P₂')/2。

10.如权利要求7所述方法，其特征在于，所述步骤S2.6.5中的大气球面分层的Snell定律为：当光线由折射指数为n₁的大气层，射向折射指数为n₂的大气层时，满足下面的关系：

n₁sinφ₁＝n₂sinφ₂

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