CN110081905B - 一种基于单站光电经纬仪的光波大气折射误差计算方法 - Google Patents

Abstract

本发明提供了一种基于单站光电经纬仪的光波大气折射误差计算方法，其特征在于，还需要雷达对飞行目标的距离进行测量，具体包括以下步骤：S1)建立站心坐标系，将光电经纬仪位置、雷达位置转换到站心坐标系；……S4)从光电经纬仪的位置起将

按n(r)的分层数量分为L段，计算光电经纬仪与飞行目标之间的视在距离

S5)从光电经纬仪的位置起将r_M按n(r)的分层数量分为L段计算光电经纬仪与飞行目标之间的地心角φ；S6)根据下式(V)计算飞行目标相对于光电经纬仪的真实仰角E′。本发明的光波大气折射误差计算方法能够在单站光电经纬仪观测数据条件下，同时通过优化视在距离、地心角定积分计算方法，进行精度更高的光波大气折射误差计算。

Description

一种基于单站光电经纬仪的光波大气折射误差计算方法

技术领域

本发明属于光电测量技术领域，具体涉及一种基于单站光电经纬仪的光波大气折射误差计算方法。

背景技术

光电经纬仪是一种带伺服系统的大型光学测量设备，具有对空中运动目标进行自动跟踪的功能。主要由光电经纬仪组成的光学测量系统，适用于对空中轨迹段较长的运动目标进行目标观测、轨迹测量等工作。光波通过大气传播时，因大气层折射率的不同，其传播速度已不是匀速直线运动，路径也发生了弯曲，这就引起光波信号的折射误差。实际观测数据的计算结果表明，当目标到探测站距离为100km左右，俯仰角为5°～10°时，折射误差可达几十米，高低角数据的折射误差在1角分以上。因此，在实际观测数据处理中，光电经纬仪观测数据必须进行大气折射修正。

若光电测量设备安装了激光测距仪，可采用射线描迹法进行折射误差修正(参考文献：成求青，等.导弹测控系统总体设计原理与方法[M].北京清华大学出版社，2014.；张瑜，甘利萍.大气折射误差修正研究现状与展望[J].河南师范大学学报(自然科学版)，2016，44(4)：45-52.)，该方法的基础是射线理论的费马原理，在大气球面分层条件下导出光波射线满足的Snell定理，再根据光电设备和目标的几何关系推导出光波射线描迹基本公式，从而推出大气折射误差公式。同时，在实用天文学和大地测量学中主要研究大气折射映射函数，将天文大气折射、大气延迟表示为映射函数与相应参数的乘积，映射函数的精度确定大气折射误差计算精度(参考文献：严豪健.大气折射母函数方法、大气折射解析解和映射函数[J].世界科技研究与发展，2016(01)：851-862)。而光电经纬仪提供设备指向目标的方位和俯仰角，不直接提供目标距离信息，在GJB 2234A-2014中采用两台设备的观测值，用方向余弦的方法进行交会，得到目标的初始位置，按照无测距数据时修正公式进行光波折射误差修正(参考文献：中华人民共和国国家军用标准.光电经纬仪事后数据处理方法[S].GJB 2234A-2014)。该方法需要两台光电经纬仪同时提供测量数据，而光电经纬仪为大型贵重设备，在试验的同时段、同区域存在仅有单台设备的情况时，该方法不适用。同时，在GJB 2234A-2014(4.10光波折射误差修正)中视在距离、地心角积分计算公式采用简化近似公式(GJB 2234A-2014附录C)，该方法计算的折射误差精度不高。

发明内容

本发明的发明目的是提供一种基于单站光电经纬仪的光波大气折射误差计算方法，能够在单站光电经纬仪观测数据条件下，同时通过优化视在距离、地心角定积分计算方法，进行精度更高的光波大气折射误差计算。

本发明的具体技术方案是一种基于单站光电经纬仪的光波大气折射误差计算方法，用于对从地面发射的飞行目标进行跟踪测量，其特征在于，还需要雷达对飞行目标进行测距，具体包括以下步骤：

S1)建立站心坐标系，将光电经纬仪位置、雷达位置转换到站心坐标系；

S2)建立光电经纬仪上空的光波折射指数剖面分布n(r)，r为地心距；

S3)利用光电经纬仪测量到的飞行目标的方位角A和俯仰角E以及雷达测距值R，进行交会计算，得到飞行目标的坐标位置，计算飞行目标近似地心距

和光电经纬仪与飞行目标之间的斜距R'₁；

S4)从光电经纬仪的位置起将

按n(r)的分层数量分为L’段，按下式(I)计算光电经纬仪与飞行目标之间的视在距离

其中，

按下式(II)计算，

式中，r1_i为各分段节点的地心距，i＝0,1,…,L’，r₀为光电经纬仪所在位置的地心距，n₀为光电经纬仪所在位置的地面光波折射指数，

为第i段对应的n(r)的分层的光波折射指数，x_r1为积分变量，

S41)进行迭代计算，若

则飞行目标地心距

进入S44),

S42)若

则令

若

则令

其中l为迭代次数，

第一次迭代时，ε⁽¹⁾为

所在n(r)的分层的层距，

S43)将新得到的

代入公式(I)中，重新计算

返回S41)，

S44)结束迭代；

S5)从光电经纬仪的位置起将r_M按n(r)的分层数量分为L”段，按下式(III)计算光电经纬仪与飞行目标之间的地心角φ，

其中，Δφ_i按下式(IV)计算，

式中，r2_i为各分段节点的地心距，i＝0,1,…,L”，

为第i段对应的n(r)的分层的光波折射指数，x_r2为积分变量；

S6)根据下式(V)计算飞行目标相对于光电经纬仪的真实仰角E′，

最后得到折射修正量ΔE＝E-E'。

更进一步地，所述的步骤S1)中，采用飞行目标的起飞点为原点建立站心坐标系。

更进一步地，所述的步骤要S2)中，光电经纬仪上空的光波折射指数剖面分布n(r)可以采用折射率分布模型或气象数据，根据气象数据按下式(VI)计算光波折射指数：

n(r)＝1+N(r)·10^-6……(VI)

其中，N(r)按下式(VII)计算，

式中，P_r为压强hPa，t_r为气温，单位℃。

更进一步地，所述的步骤3)中进行进行交会计算的方法是，采用最小二乘法的迭代算法，

将观测值方程组线性化得到以下公式(VIII)，

ΔL＝CΔD+ζ……(VIII)

式中，ΔL＝[A-A⁰,E-E⁰,R-R⁰]^T，A⁰、E⁰为为飞行目标相对于光电经纬仪的方位角和俯仰角的迭代初值，R⁰为飞行目标相对于雷达的距离的迭代初值，A⁰、E⁰、R⁰由飞行目标的坐标的迭代初值x⁰、y⁰、z⁰计算得到，ΔD＝[Δx,Δy,Δz]^T为飞行目标轨迹函数在三个坐标轴上的微分，C为A、E、R对飞行目标变量(x、y、z)偏导数的雅可比矩阵，ζ是光电经纬仪和雷达的观测数据的随机误差向量，即

设各观测数据的随机误差不相关，则其协方差阵记为

由高斯－马尔科夫估计法可得ΔD＝(C^TP^-1C)^-1C^TP^-1ΔL，则飞行目标的位置坐标

按下式(IX)计算，

式中，

[x⁰]＝(x⁰,y⁰,z⁰)^T，

进行迭代计算，如果结果不收敛，将公式(IX)计算得到的飞行目标的位置坐标赋与[x⁰]，再次按公式(IX)进行计算。

更进一步地，所述的A⁰、E⁰、R⁰由飞行目标的坐标的迭代初值x⁰、y⁰、z⁰计算得到的公式为下式(X)、(XI)、(XII)，

R⁰＝[(x⁰-x₂)²+(y⁰-y₂)²+(z⁰-z₂)²]^1/2……(X)

E⁰＝arctan[(y⁰-y₁)/L₁]……(XII)

式中，L₁＝[(x⁰-x₁)²+(z⁰-z₁)²]^1/2，(x₁,y₁,z₁)为光电经纬仪的坐标，(x₂,y₂,z₂)为雷达的坐标。

更进一步地，所述的

其中，

本发明的有益效果是1)本发明的光波大气折射误差计算方法，采用单站光电经纬仪测量数据结合雷达测距数据解算光波大气折射误差，能够适用于仅有单台光电经纬仪设备或在试验的同时段、同区域存在仅有单台设备数据的情况；2)对视在距离、地心角的定积分计算进行了优化，同现有标准的简化近似公式相比，提高了折射误差计算精度；3)利用最小二乘法给出求目标坐标的迭代算法，具有计算精度高、收敛速度快的优点，相对于现有直接算法，能够提离飞行目标的测算精度。

附图说明

图1为本发明所述的基于单站光电经纬仪的光波大气折射误差计算方法的流程简图；

图2为光电经纬仪与雷达交会示意图，其中O1为光电经纬仪，O2为雷达，T为交会目标位置，A为光电经纬仪方位角，E为光电经纬仪俯仰角，R为雷达测距值；

图3为本发明的光波大气折射误差计算方法计算的折射误差。

具体实施方式

下面结构说明书附图对本发明的具体技术方案作进一步地描述。

如附图1所示，为了更形象说明本发明的实现方式，采用单台光电光电经纬仪对飞行目标进行跟踪测量，雷达提供对目标的距离测量数据，飞行目标上搭载了高精度的GPS设备(该设备定位数据作为轨迹标准源)，光电经纬仪与雷达交会示意如图2所示。本发明的一种基于单站光电经纬仪的光波大气折射误差计算方法，具体包括以下步骤：

S1)建立站心坐标系，将光电经纬仪位置、雷达位置转换到站心坐标系。站心坐标系的坐标原点可以是飞行目标的起飞点；

S2)建立光电经纬仪上空的光波折射指数剖面分布n(r)，r为地心距。光电经纬仪上空的光波折射指数剖面分布n(r)可以采用折射率分布模型或气象数据，根据气象数据按下式(VI)计算光波折射指数：

n(r)＝1+N(r)·10^-6……(VI)

其中，N(r)按下式(VII)计算，

式中，P_r为压强hPa，t_r为气温，单位℃，r＝r₀+H，单位m；r₀为测站的地心距，单位m；H为离测站的高度，单位m；

S3)利用光电经纬仪测量到的飞行目标的方位角A和俯仰角E以及雷达测距值R，进行交会计算，得到飞行目标的坐标位置，计算飞行目标近似地心距

和光电经纬仪与飞行目标之间的斜距R'₁，

r_b＝R₀+h_b，r_b发射坐标系原点的地心距；R₀地球平均半径；h_b发射坐标系原点的大地高，

S4)从光电经纬仪的位置起将

按n(r)的分层数量分为L’段，按下式(I)计算光电经纬仪与飞行目标之间的视在距离

其中，

按下式(II)计算，

式中，r1_i为各分段节点的地心距，i＝0,1,…,L’，r₀为光电经纬仪所在位置的地心距，n₀为光电经纬仪所在位置的地面光波折射指数，

为第i段对应的n(r)的分层的光波折射指数，x_r1为积分变量，

S41)进行迭代计算，若

则飞行目标地心距

进入S44),

S42)若

则令

若

则令

其中l为迭代次数，

第一次迭代时，ε⁽¹⁾为

所在n(r)的分层的层距，

S43)将新得到的

代入公式(I)中，重新计算

返回S41)，

S44)结束迭代；

S5)从光电经纬仪的位置起将r_M按n(r)的分层数量分为L”段，按下式(III)计算光电经纬仪与飞行目标之间的地心角φ，

其中，Δφ_i按下式(IV)计算，

式中，r2_i为各分段节点的地心距，i＝0,1,…,L”，

为第i段对应的n(r)的分层的光波折射指数，x_r2为积分变量；

S6)根据下式(V)计算飞行目标相对于光电经纬仪的真实仰角E′，

最后得到折射修正量ΔE＝E-E'。

所述的步骤3)中进行进行交会计算的方法是，采用最小二乘法的迭代算法，

将观测值方程组线性化得到以下公式(VIII)，

ΔL＝CΔD+ζ……(VIII)

式中，ΔL＝[A-A⁰,E-E⁰,R-R⁰]^T，A⁰、E⁰为为飞行目标相对于光电经纬仪的方位角和俯仰角的迭代初值，R⁰为飞行目标相对于雷达的距离的迭代初值，A⁰、E⁰、R⁰由飞行目标的坐标的迭代初值x⁰、y⁰、z⁰计算得到，ΔD＝[Δx,Δy,Δz]^T为飞行目标轨迹函数在三个坐标轴上的微分，C为A、E、R对飞行目标变量(x、y、z)偏导数的雅可比矩阵，ζ是光电经纬仪和雷达的观测数据的随机误差向量，即

设各观测数据的随机误差不相关，则其协方差阵记为

由高斯－马尔科夫估计法可得ΔD＝(C^TP^-1C)^-1C^TP^-1ΔL，则飞行目标的位置坐标

按下式(IX)计算，

式中，

[x⁰]＝(x⁰,y⁰,z⁰)^T，

进行迭代计算，如果结果不收敛，将公式(IX)计算得到的飞行目标的位置坐标赋与[x⁰]，再次按公式(IX)进行计算。

所述的A⁰、E⁰、R⁰由飞行目标的坐标的迭代初值x⁰、y⁰、z⁰计算得到的公式为下式(X)、(XI)、(XII)，

R⁰＝[(x⁰-x₂)²+(y⁰-y₂)²+(z⁰-z₂)²]^1/2……(X)

E⁰＝arctan[(y⁰-y₁)/L₁]……(XII)

式中，L₁＝[(x⁰-x₁)²+(z⁰-z₁)²]^1/2，(x₁,y₁,z₁)为光电经纬仪的坐标，(x₂,y₂,z₂)为雷达的坐标。

所述的

其中，

将飞行器遥测数据作为轨迹基准数据，计算目标相对于光电经纬仪的真实俯仰角，进而评估折射误差修正效果。图3为本发明所述单站光电经纬仪测量计算的折射误差；横轴是观测时间，纵轴(左)为大气折射误差，纵轴(右)为光电经纬仪的真实俯仰角；该图可直观了解不同俯仰角对应的大气折射误差效果。表1所示为抽取的俯仰角E与大气折射误差ΔE统计表。

表1抽取的俯仰角E与大气折射误差ΔE统计表

E/°	2.7417	4.6361	6.2538	8.8577	9.5487	11.4270
							ΔE/”	40.339	43.165	41.944	39.674	39.1302	37.981
E/°	16.7090	20.7903	24.7656	26.4063	30.8635	33.2855
							ΔE/”	33.259	29.106	25.833	24.529	20.094	21.393

从图3可以看到，光电经纬仪在仰角比较低时，大气折射误差ΔE明显较大；而在仰角比较高时，大气折射误差ΔE明显减小；上述情况在表1中也能直观反映出来。0～10°俯仰角时，大气折射误差峰值高于43″，范围在39″～44″；俯仰角大于30°，大气折射误差在21秒附近。从上面的计算结果看，该方法对单站光电经纬仪的大气折射误差能够进行定性分析。

本实施案例提供一种光电经纬仪测量中对光波大气折射误差的计算方法，针对仅有单台光电经纬仪设备或在试验的同时段、同区域存在仅有单台设备数据的情况，结合雷达测距数据能够定性分析单站光电经纬仪的大气折射误差。针对大气折射修正“视在距离”、“地心角”理论公式进行优化，优化后算法在时间效率上优于常用数值积分方法；优化后的“视在距离”、“地心角”定积分计算公式同GJB 2234A-2014中简化近似公式相比，提高了折射误差计算精度。

虽然本发明已经以较佳实施例公开如上，但实施例并不是用来限定本发明的。在不脱离本发明之精神和范围内，所做的任何等效变化或润饰，同样属于本发明之保护范围。因此本发明的保护范围应当以本申请的权利要求所界定的内容为标准。

Claims (6)

1.一种基于单站光电经纬仪的光波大气折射误差计算方法，用于对从地面发射的飞行目标进行跟踪测量，其特征在于，还需要雷达对飞行目标进行测距，具体包括以下步骤：

S1)建立站心坐标系，将光电经纬仪位置、雷达位置转换到站心坐标系；

S2)建立光电经纬仪上空的光波折射指数剖面分布n(r)，r为地心距；

S3)利用光电经纬仪测量到的飞行目标的方位角A和俯仰角E以及雷达测距值R，进行交会计算，得到飞行目标的坐标位置，计算飞行目标近似地心距

和光电经纬仪与飞行目标之间的斜距R₁'；

S4)从光电经纬仪的位置起将

按n(r)的分层数量分为L’段，按下式(I)计算光电经纬仪与飞行目标之间的视在距离

其中，

按下式(II)计算，

式中，r1_i为各分段节点的地心距，i＝0,1,…,L’，r₀为光电经纬仪所在位置的地心距，n₀为光电经纬仪所在位置的地面光波折射指数，

为第i段对应的n(r)的分层的光波折射指数，x_r1为积分变量，

S41)进行迭代计算，若

则飞行目标地心距

进入S44),

S42)若

则令

若

则令

其中l为迭代次数，

第一次迭代时，ε⁽¹⁾为

所在n(r)的分层的层距，

S43)将新得到的

代入公式(I)中，重新计算

返回S41)，

S44)结束迭代；

S5)从光电经纬仪的位置起将r_M按n(r)的分层数量分为L”段，按下式(III)计算光电经纬仪与飞行目标之间的地心角φ，

其中，Δφ_i按下式(IV)计算，

式中，r2_i为各分段节点的地心距，i＝0,1,…,L”，

为第i段对应的n(r)的分层的光波折射指数，x_r2为积分变量；

S6)根据下式(V)计算飞行目标相对于光电经纬仪的真实仰角E′，

最后得到折射修正量ΔE＝E-E'。

2.根据权利要求1所述的一种基于单站光电经纬仪的光波大气折射误差计算方法，其特征在于，所述的步骤S1)中，采用飞行目标的起飞点为原点建立站心坐标系。

3.根据权利要求1所述的一种基于单站光电经纬仪的光波大气折射误差计算方法，其特征在于，所述的步骤S2)中，光电经纬仪上空的光波折射指数剖面分布n(r)可以采用折射率分布模型或气象数据，根据气象数据按下式(VI)计算光波折射指数：

n(r)＝1+N(r)·10^-6……(VI)

其中，N(r)按下式(VII)计算，

式中，P_r为压强hPa，t_r为气温，单位℃。

4.根据权利要求1所述的一种基于单站光电经纬仪的光波大气折射误差计算方法，其特征在于，所述的步骤3)中进行交会计算的方法是，采用最小二乘法的迭代算法，

将观测值方程组线性化得到以下公式(VIII)，

ΔL＝CΔD+ζ……(VIII)

式中，ΔL＝[A-A⁰,E-E⁰,R-R⁰]^T，A⁰、E⁰为飞行目标相对于光电经纬仪的方位角和俯仰角的迭代初值，R⁰为飞行目标相对于雷达的距离的迭代初值，A⁰、E⁰、R⁰由飞行目标的坐标的迭代初值x⁰、y⁰、z⁰计算得到，ΔD＝[Δx,Δy,Δz]^T为飞行目标轨迹函数在三个坐标轴上的微分，C为A、E、R对飞行目标变量(x、y、z)偏导数的雅可比矩阵，ζ是光电经纬仪和雷达的观测数据的随机误差向量，即

设各观测数据的随机误差不相关，则其协方差阵记为

由高斯－马尔科夫估计法可得ΔD＝(C^TP^-1C)^-1C^TP^-1ΔL，则飞行目标的位置坐标

按下式(IX)计算，

式中，

[x⁰]＝(x⁰,y⁰,z⁰)^T，

进行迭代计算，如果结果不收敛，将公式(IX)计算得到的飞行目标的位置坐标赋与[x⁰]，再次按公式(IX)进行计算。

5.根据权利要求4所述的一种基于单站光电经纬仪的光波大气折射误差计算方法，其特征在于，所述的A⁰、E⁰、R⁰由飞行目标的坐标的迭代初值x⁰、y⁰、z⁰计算得到的公式为下式(X)、(XI)、(XII)，

R⁰＝[(x⁰-x₂)²+(y⁰-y₂)²+(z⁰-z₂)²]^1/2……(X)

E⁰＝arctan[(y⁰-y₁)/L₁]……(XII)

式中，L₁＝[(x⁰-x₁)²+(z⁰-z₁)²]^1/2，(x₁,y₁,z₁)为光电经纬仪的坐标，(x₂,y₂,z₂)为雷达的坐标。

6.根据权利要求4所述的一种基于单站光电经纬仪的光波大气折射误差计算方法，其特征在于，所述的

其中，

Images