Correction of atmospheric refraction error based on correlation analysis of star measurement residual data

靶场采用光电经纬仪测角数据交汇测量的方法得到空中飞行目标的外弹道参数^[1-4]。因此，光电经纬仪的测角精度^[5-7]是其非常重要的指标。通常采用夜间拍星^[8-11]的方式计算光电经纬仪的静态测角误差。拍星前，需要对光电经纬仪进行调平，并测量垂直轴倾斜角、照准差、零位差等参数，用于修正光电经纬仪对星体的测角数据。根据星体相对于光电经纬仪的理论方位、俯仰角和经纬仪测量值计算得到测量残差，光电经纬仪的静态测角总误差、系统误差和随机误差都是基于测量残差的统计量^[8]。

目前，对拍星流程、星体理论值计算、经纬仪测角精度计算的研究较多，相关的理论也比较成熟。但是，对俯仰角的大气折光差^[12]修正均采用经验公式^[6,8]，由于经验公式中的参数统计区域与实际拍星地域和大气状态的不一致，导致实际拍星解算得到的俯仰角测角总误差超差。且由于大气折光差修正误差较大，影响了测量残差中其他误差因素的表征，不利于外场通过拍星方法分析导致设备测角精度超差的误差因素。

为了进行准确的大气折光差修正，通常需要使用探空气球、气象飞机采集不同高度层大气参数^[13-14]。这种方法组织复杂，难以在靶场外场条件下常态化开展。为此，文中提出了一种基于多台光电经纬仪同步拍星数据相关性^[15]分析的大气折光差修正方法。根据经纬仪拍星方位角和俯仰角测量残差模型推导得到了大气折光差修正误差模型。根据该误差模型，利用分布在同一区域不同点位的多台光电经纬仪拍星俯仰角残差数据，采用最小二乘法拟合^[16]得到大气折光差的修正系数并修正俯仰角测量残差数据。实测数据表明，采用该方法对俯仰角测量残差进行大气折光差修正后，光电经纬仪俯仰角测角总误差显著降低，且由垂直轴倾斜误差修正错误引起的方位角和俯仰角残差特性得以显现。文中提出的方法无需使用探空气球等获取大气参数，即可对同一地域分布的多台光电经纬仪拍星俯仰角残差数据进行修正，修正后计算得到的设备俯仰角测角总误差更为合理可靠，且修正后的数据可用于进一步分离其它误差因素，具有较强的工程应用价值。

光电经纬仪拍摄的第i颗星第j画幅的方位角A_ij和俯仰角E_ij的测量值为^[8]：

式中：$ A_{ij}'' $为第i颗星第j画幅的编码器方位角；$ E_{ij}'' $为第i颗星第j画幅的编码器俯仰角；$ \Delta {a_{ij}} $为第i颗星第j画幅的方位角脱靶量；$ \Delta {e_{ij}} $为第i颗星第j画幅的俯仰角脱靶量；c为照准差；b为水平轴倾斜角；I为垂直轴倾斜角；g为方位角定向差；h为零位差；a_h为垂直轴倾斜方位角；$ \Delta {P_{ij}} $为大气折光差修正值。

第i颗星第j画幅的方位角A_ij和俯仰角E_ij的测量残差为：

式中：$ A_{ij}^0 $和$ E_{ij}^0 $分别为星体方位、俯仰角理论值^[8]。综合公式（1）和公式（2）即为方位角和俯仰角拍星测量残差模型。

光电经纬仪方位测角总误差$ {\sigma _A} $、俯仰测角总误差$ {\sigma _E} $以及系统误差和随机误差都是基于测量残差$ \Delta {A_{ij}} $和$ \Delta {E_{ij}} $的统计量。理想情况下，测量残差应该是均值为0的随机序列。但是，从公式（1）可以看出，轴系误差计算错误、编码器误差计算错误、大气折光差修正错误等都可能导致星体角度测量值出现误差，从而给测量残差序列带来非随机性。

假设垂直轴倾斜角真值为I₀，然而实际代入公式（1）用于误差修正的垂直轴倾斜角测量值为I，则垂直轴倾斜角修正错误给方位角测量残差带来的修正误差为：

式中：$ \Delta I = I - {I_0} $。

垂直轴倾斜角修正错误给俯仰角测量残差带来的修正误差为：

公式（3）和公式（4）即为垂直轴倾斜角修正误差模型。

根据垂直轴倾斜角修正误差模型，垂直轴倾斜角修正错误引起的方位与俯仰角残差是成对出现的，且随星体方位角的变化分别呈现正弦和余弦变化规律。因此，当测量残差的主要误差因素是垂直轴倾斜角修正误差时，方位与俯仰方向测量残差的互相关曲线应符合$ \sin $函数与$ -\cos $函数的互相关曲线特征。该特征为在0处相关性为0，在相位$ \pm {\text{9}}{{\text{0}}^{\text{o}}} $处相关性为正负极大值，且互相关曲线关于零点旋转对称。

根据上述特征，定义方位和俯仰角测量残差归一化互相曲线$ {M_{corr}} $，该曲线的峰值为$ {M_{\max }} $，该值越大，说明方位和俯仰角测量残差相关性越强。

大气折光差修正值通常采用经验公式计算，该公式为^[1]：

式中：P₀为光电经纬仪处大气压强；T₀为大气温度。

公式（5）表明，大气折光差修正值是一个与星体俯仰角余切成正比例关系的数。但是该公式中的比例系数是经验值，当光电经纬仪所处地域大气环境特征与该经验值统计区域存在较大差异时，该公式给出的俯仰角修正值误差较大。

假设大气折光差修正值与星体俯仰角余切的比例系数存在误差量$ \Delta k $，则该误差引起的大气折光差修正误差为：

式中：$ \Delta k $为大气折光差修正系数。上式表明大气折光差修正误差与星体俯仰角余切线性相关。

上述大气折光差修正误差模型是基于公式（1）推导得到的。但是，公式（1）描述的光电经纬仪测角误差修正模型只考虑了经纬仪自身轴系误差、编码器误差和大气折射率的影响。在靶场实际工程应用中，往往容易出现经纬仪基座晃动、机械结构松动、变形、光机结构调整不到位等误差因素。因此，如果只有1台光电经纬仪拍星的俯仰角残差数据与俯仰角余切呈现正相关性，很难判断这是由大气折光差修正错误引起的还是由其他尚未发现的，在公式（1）中没有考虑到的误差因素引起的。

如果同一区域同一时段有多台光电经纬仪拍星，且该区域内大气环境较为稳定、一致，那么大气折光差修正错误对各台设备的影响规律是一致的，各设备的俯仰角测量残差均会表现出与星体俯仰角余切的相关性，且各设备的大气折光差修正系数$ \Delta k $应较为一致。因此，对每台设备的俯仰测角残差数据，将$ \cot (E_{ij}'') $作为自变量，$ \Delta {E_{ij}} $作为因变量，采用最小二乘法进行线性拟合，得到每台设备的大气折光差修正系数$ \Delta k $。取其平均值$ \Delta \overline k $作为综合修正系数，按照下式对每台设备俯仰角测量残差进行修正：

式中：$ \Delta {E'_{ij}} $为修正后的俯仰角测量残差。此处即为大气折光差修正误差模型。

2022年6月，在西北荒漠戈壁地区组织6台光电经纬仪拍星。上述设备分布在半径$10\; {\rm{km}}$的圆域内。

拍星前，所有设备拍摄方位标，调整大盘水平并记录调整数据。拍星时间为21~23时。这段时间内拍星区域无云。首先采用正镜、倒镜分别拍摄北极星，然后根据事前编制的拍星表，统一调度所有设备按四个象限依次摄、记录同一星体图像数据。拍星结束后判读拍星图像，并根据方位标判读结果、大盘水平调整结果、北极星判读结果计算三差、垂直轴倾斜角以及水平轴倾斜误差等参数。最后根据公式（1）和公式（2）计算测量残差并统计得到测角总误差等参数。

理想情况下，如果光电经纬仪的各项系统误差和大气折光差修正正确，那么方位和俯仰角测量残差应该是均值为0的随机序列。但是，从公式（6）可以看出，如果大气折光差修正错误，俯仰角测量残差将同星体俯仰角余切值线性相关。

以2号光电经纬仪的残差数据为例，其俯仰角测角总误差为$ {\sigma _{Em = 2}} = 17.0'' $（m为设备编号），该值相较于设备设计指标明显超差。以$ \cot (E_{ij}'') $为横轴，以俯仰角残差$ \Delta {E_{ij}} $为纵轴作图，如图1所示。

Figure 1.  Pitch angle residual and atmospheric refraction error fitting curve

由此可以看出，俯仰角测量残差中存在明显的随$ \cot (E_{ij}'') $线性变化的分量。采用最小二乘法进行直线拟合，得到拟合直线如图1红色直线所示，拟合得到的直线斜率$ \Delta {k_{m = {\text{2}}}}{\text{ = 15}}{\text{.9}} $即为大气折光差修正系数。该组数据拟合残差的均方根误差（Root mean squared error, RMSE）为4.8。RMSE越小，说明仰角测量残差与星体俯仰角余切值相关性越强，拟合效果越好。

因为在西北荒漠地区夜间拍星，拍星期间各个测量点位测得的环境温度、气压等参数较为接近，且拍星期间整片区域无云静风，大气环境较为稳定，因此该片区域大气环境一致性较好。对6台设备俯仰角测量残差按照上述步骤进行线性拟合得到6个大气折光差修正系数，取其平均值$ \Delta \overline k $作为综合修正系数，按照公式(7)对每台设备俯仰角测量残差进行修正。2号光电经纬仪修正后的俯仰角测量残差$ \Delta {E'_{ij}} $随$ \cot (E_{ij}'') $的变化关系见图2。从图2可以看出，采用综合修正系数修正大气折光差后，设备2的俯仰角拍星残差$ \Delta {E'_{ij}} $与$ \cot (E_{ij}'') $的相关性显著降低，俯仰角测角总误差从17.0″降低到6.1″，符合设备指标要求。

每台设备的大气折光差修正系数$ \Delta {k_m} $，拟合残差的RMSE，综合修正系数$ \Delta \overline k $以及使用$ \Delta \overline k $修正前后的俯仰测角总误差$ {\sigma _{Em}} $和$ {\sigma '_{Em}} $如表1所示。

Figure 2.  Pitch angle residual after atmospheric refraction error correction

综合图2和表1可以看出，6台设备俯仰角残差均与$ \cot (E_{ij}'') $呈现较强的线性相关性，且拍星计算得到的6台设备俯仰角静态测角总误差相较于设备设计指标均严重超差，这显然是不合理的。从图2可以看出，按照大气折光差修正误差模型修正后的俯仰角残差与$ \cot (E_{ij}'') $的相关性消失。表1数据表明，经过修正后，6台设备俯仰角静态测角总误差显著降低，满足设备指标要求。

根据垂直轴倾斜角修正误差模型，垂直轴倾斜角修正错误引起的方位与俯仰角残差是成对出现的，且随星体方位角的变化分别呈现正弦和余弦变化规律。

图3中，蓝色点为2号光电经纬仪方位角测量残差$ \Delta {A_{ij}} $随星体方位角A_ij的变化数据点，红色实线为数据点的正弦拟合曲线。图4中，蓝色点和红色实线分别为2号光电经纬仪俯仰角测量残差$ \Delta {E_{ij}} $随星体方位角A_ij的变化数据点及其余弦拟合曲线。由此可以看出，$ \Delta {A_{ij}} $随A_ij呈现正弦变化规律，然而$ \Delta {E_{ij}} $随A_ij的余弦变化规律不明显。这与垂直轴倾斜角修正误差模型的结论不相符。

Figure 3.  $ \Delta {A_{ij}} $ versus A_ij

Figure 4.  $ \Delta {E_{ij}} $ versus A_ij

结合大气折光差修正误差模型及综合修正系数$ \Delta \overline k $不难看出，由大气折光差修正错误引起的俯仰角残差变化幅度达到$\pm {\text{12}}{\text{.9″}}$，与该设备的俯仰角残差$ \Delta {E_{ij}} $变化幅度相当。该误差的存在导致垂直轴倾斜角修正错误引起的俯仰角残差随方位角呈余弦变化的规律被淹没。

按照2.2节方法，经过大气折光差修正后的2号光电经纬仪俯仰角测量残差$ \Delta {E'_{ij}} $及其余弦拟合曲线见图5。对比图4和图5可以看出，经过大气折光差修正后，$ \Delta {E'_{ij}} $随方位角呈现明显的余弦变化规律。

Figure 5.  $ \Delta {E'_{ij}} $ versus A_ij

图6为$ \Delta {A_{ij}} $与$ \Delta {E_{ij}} $的归一化互相关曲线$ {M_{corr}}(\Delta {A_{ij}},\Delta {E_{ij}}) $和$ \Delta {A_{ij}} $与$ \Delta {E'_{ij}} $的归一化互相关曲线$ {M_{corr}}(\Delta {A_{ij}},\Delta {E'_{ij}}) $。从图6可以看出，大气折光差修正前，$ \Delta {A_{ij}} $与$ \Delta {E_{ij}} $的归一化互相关曲线峰值为0.66。而进行大气折光差修正后，$ \Delta {A_{ij}} $与$ \Delta {E'_{ij}} $的归一化互相关曲线峰值相关系数提升至0.77，即大气折光差修正后，$ \Delta {A_{ij}} $与$ \Delta {E'_{ij}} $的相关性增强，垂直轴倾斜角修正误差特性凸显。

Figure 6.  Normalized cross-correlation curves of $ \Delta {A_{ij}} $, $ \Delta {E_{ij}} $ and $ \Delta {A_{ij}} $, $ \Delta {E'_{ij}} $

6台光电经纬仪大气折光差修正前的方位角残差与俯仰角残差的归一化互相关曲线峰值$ {M_{\max }} $以及修正后的归一化互相关曲线峰值$ {M'_{\max }} $如表2所示。

从表2可以看出，大气折光差修正前，各设备的方位角测量残差和俯仰角测量残差的归一化互相关曲线峰值普遍低于0.7。即受大气折光差修正误差的影响，方位角测量残差和俯仰角测量残差的相关性降低，垂直轴倾斜角修正误差特征不显著，不利于发现该误差因素对测角总误差的影响。大气折光差修正后，各设备的方位角测量残差和俯仰角测量残差的归一化互相关曲线峰值均大于0.7，表明方位和俯仰角残差的相关性增强，该相关性主要是由垂直轴倾斜角修正误差引起的。通过对$ \Delta {A_{ij}} $与$ \Delta {E'_{ij}} $的进一步拟合、分析，可以得到垂直轴倾斜角修正误差幅度和相位等信息，可用于指导设备进行误差分离，调整大盘水平，并再次组织拍星检测系统测角精度。

文中提出了一种基于多台光电经纬仪同步拍星数据相关性分析的大气折光差修正方法。针对靶场条件下分布在一定区域内不同点位的多台光电经纬仪同步开展拍星的情况，若各设备点测得的地面大气参数一致性较好，则可利用该区域内多台光电经纬仪拍星俯仰角残差数据，根据文中提出的大气折光差修正模型，采用最小二乘法拟合得到大气折光差的修正系数并修正俯仰角测量残差数据。

实测数据和分析结果表明，设备俯仰角残差与$ \cot (E_{ij}'') $呈现较强的线性相关性，由大气折光差修正错误引起的俯仰角残差变化幅度达到$\pm {\text{12}}{\text{.9″}}$，与光电经纬仪的拍星俯仰角残差幅度相当。该误差的存在导致计算得到的俯仰角静态测角总误差严重超差，且垂直轴倾斜角修正错误引起的俯仰角残差随方位角呈余弦变化的规律被淹没，方位角和俯仰角残差归一化互相关曲线峰值均小于0.7。采用文中方法对俯仰角测量残差进行大气折光差修正后，光电经纬仪俯仰角静态测角总误差显著降低，且由垂直轴倾斜误差修正错误引起的方位角和俯仰角残差的相关性得以显现，修正后方位角和俯仰角残差归一化互相关曲线峰值大于0.7。

文中提出的方法无需使用探空气球、气象飞机等获取大气参数，即可对同一地域分布的多台光电经纬仪拍星俯仰角残差数据进行修正，剔除俯仰角测量残差中大气折光差修正错误的影响。基于修正后的残差数据计算得到的俯仰方向静态测角总误差更为合理，且修正后的数据可用于进一步分离测量残差中的垂直轴倾斜角修正误差等误差因素，具有较强的工程应用价值。后续的研究重点是根据残差数据特性，进一步分离轴系误差等误差因素，为靶场外场开展光电经纬仪静态精度鉴定与设备状态调整提供理论指导。