单向精密测距三角高程测量方法

摘要

本发明提供了一种单向精密测距三角高程测量方法，主要解决测距三角高程测量中大气垂直折光系数的确定，以及测距三角高程网平差中的确权问题。主要包括：基准点及待测点的布设及其高程确定，用于折光系数相关关系确定的样本数据获取，各观测方向折光系数的计算，各方向折光系数相关关系的确定及其方差估计，单向测距三角高程的计算，以及单向测距三角高程网的确权及平差。本发明适于实施重复及高危环境下的高程测量工作，可解决不具备双向观测条件项目的高程获取问题，或者通过确定折光系数差提高三角高程测量精度的项目，或者为节约成本而采用单向代替双向三角高程测量的项目。

说明书

技术领域

本发明属于测距三角高程测量技术领域，具体涉及一种单向精密测距 三角高程测量方法。

背景技术

虽然目前GNSS技术的持续发展为空间数据的获取提供了极大便利，在 一些特定条件下，其平面相对定位精度甚至达到了亚毫米，但由于受大地 高向正常高转换过程中高程异常获取精度及其自身测量精度等因素的影响， 限制了其在更高精度高程测量领域的推广，因此高程仍然大量采用传统的 测量手段，目前较高精度高程测量主要有两种方式，一种是水准测量，另 外一种是精密测距三角高程，其中水准测量精度最高，在一些高精度高程 测量方面几乎无可替代，所以，目前一、二等高程获取主要依靠精密水准 测量实现，但在一些地区，比如高山峡谷地区，几何水准测量高程传递是 非常困难的，部分地区甚至是不可能的，但采用测距三角高程测量就比较 容易实现，三角高程具有测定高差速度快,测量简单,不受地形条件限制等 优点,特别是在高差较大的地区优势更加明显，特别是近年来，测量设备测 角、测距精度不断提高，计算方法也不断得到完善，使得测距三角高程测 量精度有了显著提高，目前采用测距三角高程测量方法获取三等、四等高 程应用已经成熟，国内、外一些学者也开展过采用测距精密三角高程代替 二等甚至一等水准的研究乃至生产，如辽宁工程大学的沈忱、中国地质大 学的韩昀、天津水电院的周水渠等人均进行了相关研究及实践，对于进一 步发挥三角高程的优势，持续挖掘高效的三角高程测量技术及测量数据处 理方法，具有重要的研究和应用价值。

但三角高程也有其自身的缺陷,其影响因素主要有两个方面，一方面受 到垂直角测量精度的限制，另一方面主要受到大气折射及垂线偏差因素的 影响。垂直角测量和大气垂直折光差是精密测距三角高程的主要误差源(应 用中，测距边长一般较短，故垂线偏差影响未予单独考虑)，随着观测仪 器精度的不断提升，特别是具有自动目标识别、自动照准、自动测角与测 距测量机器人的不断发展，垂直角观测精度已经大幅度提升，为此，大气 折射的影响成为三角高程测量重要的误差源，为了削弱折光对三角高程测 量精度的影响，一般采用对向观测的手段对折光影响进行削弱，如果地形 对称，同时段观测，确保测量的环境条件较接近，这种削弱效果是比较理 想的，但如果地形条件不对称，或者不能在同一个时间段进行往、返观测， 则在往、返测的过程中环境发生了较大变化，则难以较好抵消大气折光对 三角高程的影响，部分条件下甚至会加大误差。三角高程折光关系见附图1。 图中展示为自A点向B点观测示意图，i为仪器高，s为目标高，a为垂直 角，平距为D，高差值为h_AB。由于垂直大气折光的影响，视线发生向地面 或者反方向的折射，视线不是直线而是弧线。图中虚线即为受到大气折光 影响的视线方向。

因此，如果项目不具备对向观测条件，特别对于一些人员不便到达及 不宜暴露或不宜频繁暴露的项目，如高危滑坡体监测项目，就限制了水准 及对向测距三角高程测量的应用。当采用单向测距三角高程测量时，由于 受到大气折射等因素的影响，单向测距三角高程精度难以保证。

发明内容

本发明是提供一种单向精密测距三角高程测量方法，目的是克服背景 技术中所述单向测距三角高程测量现有技术所存在的不足。主要解决测距 三角高程测量中大气垂直折光系数(隐含垂线偏差等要素)的确定，以及 三角高程网平差中的确权问题。

本发明的技术方案是：单向精密测距三角高程测量方法，包括如下步 骤：

步骤1：布设基准点及待测点

本步骤包括如下具体步骤：

1).首先根据工程需求布设待测点，即布设观测目标点；在待测点确定 的基础上开展基准点的布设，包括用于对待测点直接观测的基准点的布设 以及用于测定实时大气垂直折光系数的基准方向点的布设；

2).确保基准点在基准点的一个复测周期内不发生显著变形，这里采用 t检验法对各基准点进行位移显著性检验以确定是否发生显著变形；

3).所布设的基准方向，应与待测方向视线所穿越的地形结构及朝向情 况接近，确保两者近地气层结构相似，其中基准方向指由基准点至基准方 向点所确定的方向，待测方向指由基准点至待测点所确定的方向；

4).如果滑坡体周边没有满足步骤2)和步骤3)要求的基准点，则在 异地独立建造辅助基准点与基准方向，要求异地所建的辅助基准点与基准 方向应与基准点至待测点方向视线所穿越的地形结构及朝向一致，确保两 地地形及朝向一致，光照条件一致；

步骤2：获取基准点及待测点高程

基准点、待测点布设完成后，联测基准点及待测点已知高程，采用水 准或者精密测距三角高程测量的方式联测，基准点高程要准确测定，测量 精度高于待测点一个等级；待测点初始高程精度需高于设计的精度要求， 对于部分获取待测点精确高程存在困难的项目，则采用单向测距三角高程 测量结果，每次反演大气垂直折光系数所采用的基准高程应一致。

步骤3：获取样本数据用于确定基准方向与待测方向大气垂直折光系数相关 关系

取多个环境下的测距三角高程测量样本数据来建立基准方向与待测方 向的相关关系；通过连续观测基准点及待测点方向测距三角高程数据，利 用步骤2所获取的高程数据，同步计算出持续的折光系数；通过对折光系 数的统计、降噪处理，建立基准方向与待测方向的相关关系；观测时间避 开阴天，以保证复杂的温度梯度变化，确保在这个观测条件下获取日常观 测需要的参数，采用等步长观测，在大气变换剧烈的情况下加密观测。

对于步骤2中提出的待测点采用单向测距三角高程的情况，需保证在 数据采集的起、止时间段内待测点不会发生显著变形，对于发生变形的点 位，可通过在起、止时间分别观测高程，按照时间内插各个时刻待测点高 程；

对于拥有多期测距三角高程测量数据的项目，直接对这些数据进行精 度及可靠性检验后，作为样本数据进行计算分析；

无论是何种方式获取的测距三角高程测量样本数据，均需要保证每组 所采集的基准方向与待测方向测距三角高程测量数据在时间上一一对应；

步骤4：确定各观测方向大气垂直折光系数

依据步骤2所提供的基准点及待测点高程值，以及步骤3得到的测距 三角高程测量样本数据，按照如下公式求解各组基准方向及待测方向大气 垂直折光系数k值：

k值计算公式如下：

k＝1-[(h-S_斜×Sinα-I+L)2R/(S_斜×cosα)²](1)

上式中：h为已知高差值，S_斜为斜边，α为垂直角，I为仪器高,L为 标高,R为参考椭球曲率半径；

对所计算的大气垂直折光系数采用小波进行降噪处理，或者在回归方 程参数估计过程中采用整体最小二乘方法进行参数估计，大气垂直折光系 数方差m_k²估算公式如下：

$\begin{array}{c}{m_k}^2={\left(\frac{2R}{{(S·\cos \alpha )}^2}\right)}^2{m_h}^2+{\left[\frac{2R}{{(S·\cos \alpha )}^2}\left(2\right(h-I+L)/S-\sin \alpha )\right]}^2{m_s}^2\\ +{\left[\frac{2R}{S^2·{\cos }^3\alpha }\left(2\right(h-I+L)sin\alpha -S({\cos }^2\alpha +2{\sin }^2\alpha \left)\right)\right]}^2{m_{\alpha }}^2+{\left(\frac{2R}{{(S·\cos \alpha )}^2}\right)}^2{m_I}^2+{\left(\frac{2R}{{(S·cos\alpha )}^2}\right)}^2{m_L}^2\end{array}---\left(2\right)$

公式中m_h²、m_s²、m_a²、m_I²、m_L²分别为对应的高差h、斜距s、垂直 角α、仪器高I、标高L的方差；

用待测点连续观测进行计算时，短期内变形点如果发生变形，则在观 测起、止时间分别测定待测点高程，通过时间内插求出指定时刻的变形点 高程；

步骤5：确定基准方向与待测方向大气垂直折光系数的相关关系及折光系数 方差估计

在步骤4确定的大气垂直折光系数k及其方差基础上，建立起基准方 向与待测方向大气垂直折光系数的相关关系模型；两者相关性建立有两种 方法：

方法一：利用回归方程建立基准方向与各个待测方向相关关系；

方法二：结合当地近地气层温度梯度分布特点，以方位、高度角要素 确定任意方向待测点大气垂直折光系数计算模型；

基于步骤4所计算的大气垂直折光系数k，利用回归方程建立基准方向 与各个待测方向的相关关系；首先通过绘制散点图，确定回归方程选配线 型，然后分析基准方向与各个待测方向大气垂直折光系数的相关性，根据 项目精度要求进行相关性满足要求的待测方向筛选，筛选的原则是精度要 求越高则对相关性要求越高；通过相关性检验后，建立基准方向与各个待 测方向折光系数的回归方程，并进行回归方程参数的计算及折光系数方差 的估算；

如果观测点较多，不便于逐点确定，则按照方法二，以方位、高度角 要素确定任意方向待测点大气垂直折光系数计算模型；

步骤6单向测距三角高程的确定

步骤5确定出了基准方向与待测方向的相关关系，每次测量过程中按 照如下公式确定基准方向大气垂直折光系数：

KH＝1-[(h-S_斜×Sinα-I+L)2R/(S_斜×cosα)²](6)

每次监测过程中，通过观测出基准方向斜距S_斜、垂直角α，仪器高I 和标高L，计算出大气垂直折光系数，

将计算出的系数KH_i代入下式：

KB_i＝a+bKH_i(7)

每次测量过程中，通过测出KH_i，就可以求出KB_i，将KB_i代入下式，确 定出各个方向的修正后三角高差值：

h_i＝S_斜×Sinα_i+I_i-L_i+(1-KB_i)×(S_斜×Cosα_i)²/2R(8)

步骤7单向测距三角高程网的确权及平差

通过两个或者多个测站对待测点进行测量，能够形成附合或者闭合路 线，则依据高程闭合差进一步修正大气垂直折光系数。计算过程如下：

设测量两个基准点至待测点B1的高差分别为h₁和h₂，两个基准点高程 为H₁和H₂，则存在如下关系：

ω＝H₁+h₁-h₂-H₂(9)

式中的ω就是闭合差，理论上ω＝0，然后将ω按照距离确权分配到高 差h₁和h₂上，再将这个数值再代回式(1)进行k值迭代计算；

然后建立高差观测值间接平差模型，确定各个高差的最优估计值，模 型如下：

$V=B\hat{x}-f---\left(13\right)$

V为高差观测值改正数；B为平差模型的系数矩阵；为高程未知参 数；-f为平差模型的常数项；由此模型确定高程参数

其中高程平差按照测距三角高程方差进行确权；

上述步骤1中所述的复测周期指由于地质运动以及工程施工等因素影 响，可能会影响到测量基准点的稳定性，需要定期进行测量控制网测量， 以便及时发现基准点的变形，并对发生变形的基准点赋予新的测量成果； 这里所说的一个复测周期就是指重复测量经历的时间间隔。

上述步骤5中，相关性确定方法具体是以最常用的线性回归方程来实 现，具体如下：

首先，求回归线的前提是变量待测点大气垂直折光系数KB与基准方向 大气垂直折光系数KH必须存在线性相关，线性相关的指标是相关系数ρ：

$\rho =\frac{\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}({\mathrm{KH}}_i-\overline{KH})({\mathrm{KB}}_i-\overline{KB})}{\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{({\mathrm{KH}}_i-\overline{KH})}^2\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{({\mathrm{KB}}_i-\overline{KB})}^2}}---\left(3\right)$

其次，令：

$\begin{array}{c}KB=\left(\begin{array}{c}K{B}_1\\ {\mathrm{KB}}_2\\ ...\\ K{B}_n\end{array}\right),KH=\left(\begin{array}{cc}1& {\mathrm{KH}}_1\\ 1& {\mathrm{KH}}_2\\ ...& ...\\ 1& {\mathrm{KH}}_n\end{array}\right),\beta =\left(\begin{array}{c}a\\ b\end{array}\right),e=\left(\begin{array}{c}{e}_1\\ e_2\\ ...\\ e_n\end{array}\right)\\ KB=KH\beta +e\end{array}---\left(4\right)$

当有效KB大于2个，需要进行参数最优估计；

根据最小二乘原理，直接给出参数估值最小二乘求解公式：

β＝(KH^TPKH)^-1KH^TPKB(5)

解法方程组得：

$\left(\begin{array}{c}\hat{a}=\frac{\Sigma {{\mathrm{KH}}_i}^2\Sigma {\mathrm{KB}}_i-\Sigma {\mathrm{KH}}_i\Sigma {\mathrm{KB}}_i{\mathrm{KH}}_i}{n\Sigma {{\mathrm{KH}}_i}^2-{(\Sigma {\mathrm{KH}}_i)}^2}\\ \hat{b}=\frac{n\Sigma {\mathrm{KB}}_i{\mathrm{KH}}_i-\Sigma {\mathrm{KB}}_i\Sigma {\mathrm{KH}}_i}{n\Sigma {{\mathrm{KH}}_i}^2-{(\Sigma {\mathrm{KH}}_i)}^2}\end{array}\right)---\left(6\right)$

用模型的复相关系数对回归模型拟合程度综合度量：

$R^2=\frac{\Sigma {(\hat{K}{B}_i-\overline{{\mathrm{KB}}_i})}^2}{\Sigma {({\mathrm{KB}}_i-\overline{{\mathrm{KB}}_{i_i}})}^2}---\left(7\right)$

其中：为待测方向大气垂直折光系数平均值，为待测方向大气 垂直折光系数计算值；R²越大，即越接近于1，则模型拟合效果越好，方程 的偏离度越小；

最后，评定其方差：

设残差值为V，则：

V＝KHβ-KB(8)

待测方向折光系数KB的方差估计值σ²为：

${\hat{\sigma }}^2=\frac{V^{t}V}{n-2}---\left(9\right)$

上述步骤7中，依据测距三角高程方差进行确权的方法如下：

单向测距三角高程方差估计公式为：

公式中m_s²、m_a²、m_I²、m_L²分别为对应的高差h、斜距s、垂 直角α、折光系数k、仪器高I及标高L的方差值，其中m_s²、m_a²、m_I²、m_L²方差采用标称或估算精度计算，直接采用待测方向折光系数KB的方差估 计值进行计算，$m_k^2={\sigma }_k^2;$

设定单位权中误差为σ₀，确定各个测段三角高程的权值为：

$P_i={\sigma }_0^2/m_{hi}^2,i=1,2,\mathrm{...},n---\left(11\right)$

本发明的有益效果：本发明是设计一套方法，以提高单向测距三角高 程测量精度。该方法主要解决测距三角高程测量中大气垂直折射系数(隐 含垂线偏差要素)的确定，以及测距三角高程网平差中的确权等问题。该 方法是依据所布置的待测点布置若干观测基准点，并通过观测获取基准点 及待测点高程。获取各组基准方向与待测方向三角高程同步观测数据，并 依据基准点与待测点高程计算各组基准方向与待测方向大气垂直折光系数， 并对所计算的各方向大气垂直折光系数进行降噪处理。根据降噪处理后的 各组大气垂直折光系数，建立基准方向与待测方向大气垂直折光系数相关 关系模型，并可估计出大气垂直折光系数方差。在应用中，当测定出基准 方向的折光系数，便可依据所建立的模型求出各个待测方向的折光系数， 将求出的待测方向折光系数代入三角高程计算公式中，便可求出基准点至 待测点高差及高差方差值。最后依据高差方差进行三角高程网确权，并进 行单程测距三角高程网或高程混合网平差。

本发明的优点如下：

1、适于重复及高危环境下的高程测量工作，比如高程监测项目，为其 提供了一套高效率、高精度的测距三角高程解决方案；

2、测量效率更高，由于不采用双向观测，可显著的降低野外数据采集 时间；

3、测量精度更高，如测量结果及选用模型可靠，可显著削弱由于往返 测环境条件不对称，所造成的折光系数抵消不彻底对结果的影响；

4、测量作业过程更安全，由于其不需要进行往返测，待测点可以放置 棱镜后，无人值守，为危险区域高精度高程数据的连续获取提供了可能；

5、测量成本更低，由于效率的显著提高，为重复高程测量提供了更低 成本的解决方案。

以下将结合附图对本发明做进一步详细说明。

附图说明

图1是三角高程折光关系示意图。

图2是模拟实验项目示意图。

图3是本发明方法的作业流程。

图4是降噪后的基准方向及待测方向大气垂直折光系数变化过程线。

具体实施方式

如图3所示为本发明的作业流程，本发明的具体实施步骤如下：

1基准方向及待测点的布设

基准方向及待测点的布设应根据项目实际地形特征而建造，主要原则 如下：

1).根据工程需求布设待测点(观测目标点)，在待测点确定的基础上 开展基准点的布设，基准点包括用于对待测点直接观测的基准点以及用于 测定实时大气垂直折光系数的基准方向点(注：文中提出的基准方向指由 基准点至基准方向点所确定的方向，待测方向指由基准点至待测点所确定 的方向)；

2).基准点应有较高的稳定性，能够确保在基准点一个复测周期内不 发生显著变形，这里采用t检验法对各基准点进行位移显著性检验以确定 是否发生显著变形；

3).所布设的基准点至基准方向，应与基准点至待测点方向视线所穿 越的地形结构及朝向情况接近，确保两者近地气层结构相似，以保证基准 方向与待测方向具有相似的温度梯度结构。

4).如果滑坡体周边没有满足以上第2条和第3条要求的稳定基准方 向点，可在异地独立建造辅助基准点与基准方向，要求异地所建的辅助基 准点与基准方向应与基准点至待测点方向视线所穿越的地形结构及朝向情 况接近，并不宜远离测区，确保两地地形及朝向接近，日照基本同步。下 边模拟一个项目以介绍普遍的布设原则。

图2是模拟实验项目示意图。

图中，白线范围内为滑坡体，按照监测计划，在变形体上建造B1、B2、 B3、B4，在这个基础上，布设观测基准点H1和H2，其中H1为观测基准点， H2为辅助基准点，基准方向H1→H2地形结构与变形点方向相似，使之间折 光关系存在相关性。

2基准点及待测点高程的获取

基准点及待测点布设完成后，联测基准点及待测点已知高程，采用水 准测量或者精密三角高程测量方法，准确测定基准点高程，测量精度应高 于待测点一个等级。对于获取待测点精确高程存在困难的项目，如仅需要 相对变形量(如滑坡变形监测等)，可采用单程测量结果，但每次反演折 光系数所采用的基准高程应一致，通过估算，只要各次反演折光系数所采 用的基准高程是一致的，其对所解算的相对折光系数变化影响量值很小。

3获取折光系数相关关系确定的样本数据

基准方向与待测方向的相关关系的建立，应获取多个环境下的观测样 本，通过连续观测基准点及待测点方向测距三角高程数据，利用步骤2所 获取的高程数据，可同步计算出持续的折光变化数据，连续观测(自上午 太阳升起之前至太阳落山后)基准点及待测点，获取持续的折光数据，对 这些数据统计并降噪处理，用于建立基准方向与待测方向的相关关系。观 测时间应避开阴天，保证复杂的温度梯度变化，确保在这个观测条件下获 取日常观测需要的参数，在大气变换剧烈的情况下加密观测。

太阳升起和太阳落山前后，折光变化迅速，推荐20分钟测量一次，每 次观测不少于4个测回(这个可以根据精度灵活确定)，其余时间可1个 小时观测一次。

对于步骤2中提出的待测点采用单向测距三角高程的情况，需保证在 数据采集的起、止时间段内待测点不会发生显著变形,对于发生变形的点位， 可通过在起、止时间分别观测高程，采用时间内插出各个时刻待测点高程。

在实际生产中，对于拥有多期测距三角高程测量数据的项目，可直接 对这些数据进行精度及可靠性检验后，作为样本数据进行计算分析。

无论是何种方式获取的测距三角高程测量样本数据，均需要保证每组 所采集的基准方向与待测方向测距三角高程测量数据在时间上一一对应。

观测仪器的选择根据观测精度要求来确定。

4大气垂直折光系数的确定

针对工程中常用的短距离测距三角高程，垂直角折光角可转换为垂直 折光系数。

两点间测距三角高程测量单向高程基本公式为： h＝S_斜×Sinα+I-L+(1-k)×(S_斜×Cosα)²/2R(1)

式(1)中：S_斜为斜边，α为垂直角，仪器高I、标高L系观测量，曲 率半径R为已知。

依据步骤2所提供的基准点及待测点高程值，以及步骤3得到的测距 三角高程测量样本数据，按照如下公式求解各组基准方向(含辅助基准方 向)及待测方向大气垂直折光系数k值：

k＝1-[(h-S_斜×Sinα-I+L)2R/(S_斜×cosα)²](2)

其中，h为已知高差值，S_斜为斜边，α为垂直角，I为仪器高,L为标 高,R为参考椭球曲率半径。

这样就可以确定出各个时间点选定方向的折光系数，由于垂直角，边 长等量观测依然存在误差，导致折光系数不可避免的含有噪声误差，必要 的情况下需要进行降噪处理，在实验阶段的案例中我们对计算的折光系数 采用小波进行了降噪处理。或者在回归方程参数估计过程中，采用整体最 小二乘方法进行参数估计，这就需要估算出折光系数的方差，按照误差传 播率，推出折光系数k方差估算公式如下：

$\begin{array}{c}{m_k}^2={\left(\frac{2R}{{(S·\cos \alpha )}^2}\right)}^2{m_h}^2+{\left[\frac{2R}{{(S·\cos \alpha )}^2}\left(2\right(h-I+L)/S-\sin \alpha )\right]}^2{m_s}^2\\ +{\left[\frac{2R}{S^2·{\cos }^3\alpha }\left(2\right(H-I+L)sin\alpha -S({\cos }^2\alpha +2{\sin }^2\alpha \left)\right)\right]}^2{m_{\alpha }}^2+{\left(\frac{2R}{{(S·\cos \alpha )}^2}\right)}^2{m_I}^2+{\left(\frac{2R}{{(S·cos\alpha )}^2}\right)}^2{m_L}^2\end{array}---\left(3\right)$

公式中m_s²、m_a²、m_I²、m_L²分别为对应的高差h、斜距s、垂 直角α、折光系数k、仪器高I及标高L的方差值。

用待测点连续观测进行计算时，短期内变形点一般不会发生显著变形， 如发生变形，可在观测起、止时间分别测定待测点高程，通过时间内插求 出指定时刻的变形点高程。

5基准方向与待测方向大气垂直折光系数相关关系的确定

基准方向与待测方向大气垂直折光系数相关关系的确定及其方差估计

在步骤4确定的大气垂直折光系数k及其方差基础上，建立起基准方 向与待测方向大气垂直折光系数的相关关系模型。两者相关性建立有两种 方法：

1利用回归方程建立基准方向与各个待测方向相关关系；

2结合当地近地气层温度梯度分布特点，以方位、高度角要素确定任 意方向待测点大气垂直折光系数计算模型。

利用回归方程建立基准方向与各个待测方向相关关系，利用步骤4所 计算的大气垂直折光系数k，通过绘制散点图，选配回归方程线型，然后分 析基准方向与各个待测方向大气垂直折光系数的相关性，线性相关的指标 是相关系数ρ，这里用如下式子计算：

如果两站之间地形结构对称，又是在同一个时刻进行往返测，那么就 可以有效削弱折光影响，否则往返测对折光的削弱有限。采用回归方程建 立相关关系方法如下：

首先通过绘制散点图，根据散点图选配线型，并建立回归方程，方程 建立应尽量简单，否则会给模型可靠性评价带来更大的难度。

一般称$\left(\begin{array}{c}Y=X\beta +ϵ\\ E\left(ϵ\right)=0,COV(ϵ,ϵ)={\sigma }^2{I}_n\end{array}\right)$

为高斯—马尔柯夫线性模型(k元线性回归模型)，并简记为(Y,Xβ,σ²I_n)

下边我们以最常用的线性回归方程为例介绍其实现过程。

设定基准方向大气垂直折光系数为KH，变形点或者待测点折光系数为 KB，线性关系可设定为如下的形式：

KB_i＝a+bKH_i+e_i(4)

其中：a、b为系数，e_i为随机误差。KB与KH服从一元线性回归模型。

求回归直线的前提是变量KB与KH必须存在线性相关，否则所配直线 就无实际意义，线性相关的指标是相关系数ρ，这里用如下式子计算：

$\rho =\frac{\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}({\mathrm{KH}}_i-\overline{KH})({\mathrm{KB}}_i-\overline{KB})}{\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{({\mathrm{KH}}_i-\overline{KH})}^2\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{({\mathrm{KB}}_i-\overline{KB})}^2}}---\left(5\right)$

|ρ|越接近于1，线性回归分析的效果越好；特别，当|ρ|＝1时， (KB₁,KH₁),(KB₂,KH₂),…,(KB_n,KH_n)将全部落在直线KB＝a+bKH上。

一般认为：

如果ρ∈[0.75，1]，那么正相关很强；如果ρ∈[-1，-0.75]，那么负相关很 强；如果ρ∈(-0.75，-0.30]或r∈[0.30，0.75)，那么相关性一般；如果 ρ∈[-0.25，0.25]，那么相关性较弱。

通过相关性检验后，令：

$\begin{array}{c}KB=\left(\begin{array}{c}K{B}_1\\ {\mathrm{KB}}_2\\ ...\\ K{B}_n\end{array}\right),KH=\left(\begin{array}{cc}1& {\mathrm{KH}}_1\\ 1& {\mathrm{KH}}_2\\ ...& ...\\ 1& {\mathrm{KH}}_n\end{array}\right),\beta =\left(\begin{array}{c}a\\ b\end{array}\right),e=\left(\begin{array}{c}{e}_1\\ e_2\\ ...\\ e_n\end{array}\right)\\ KB=KH\beta +e\end{array}---\left(6\right)$

当有效KB大于2个，需要进行参数最优估计。

根据最小二乘原理，直接给出参数估值最小二乘求解公式：

β＝(KH^TPKH)^-1KH^TPKB(7)

解法方程组得：

$\left(\begin{array}{c}\hat{a}=\frac{\Sigma {{\mathrm{KH}}_i}^2\Sigma {\mathrm{KB}}_i-\Sigma {\mathrm{KH}}_i\Sigma {\mathrm{KB}}_i{\mathrm{KH}}_i}{n\Sigma {{\mathrm{KH}}_i}^2-{(\Sigma {\mathrm{KH}}_i)}^2}\\ \hat{b}=\frac{n\Sigma {\mathrm{KB}}_i{\mathrm{KH}}_i-\Sigma {\mathrm{KB}}_i\Sigma {\mathrm{KH}}_i}{n\Sigma {{\mathrm{KH}}_i}^2-{(\Sigma {\mathrm{KH}}_i)}^2}\end{array}\right)---\left(8\right)$

下面用模型的复相关系数对回归模型拟合程度综合度量：

$R^2=\frac{\Sigma {(\hat{K}{B}_i-\overline{{\mathrm{KB}}_i})}^2}{\Sigma {({\mathrm{KB}}_i-\overline{{\mathrm{KB}}_{i_i}})}^2}---\left(9\right)$

其中：为待测方向大气垂直折光系数平均值，为待测方向大气 垂直折光系数计算值。R²越大(越接近于1)，则模型拟合效果越好，方程 的偏离度越小。

然后将上式参数代入式(3)就可确定出基准方向与待测点方向折光关 系式。然后评定其中误差：

设残差值为V，则：

V＝KHβ-KB(10)

待测方向折光系数KB的方差估计值σ²为：

${\hat{\sigma }}^2=\frac{V^{t}V}{n-2}---\left(11\right)$

方差项可以较直观反映出模型计算的精度情况。

另外统计出方差分析的F统计量

$F=\frac{SSR/f_R}{SSE/f_E}=\frac{SSR}{SSE/(n-2)}---\left(12\right)$

方差分析的显著性概率

p＝P(F(1,n-2)＞F)(13)

以上数据计算及分析可在Matlab等软件实现。

如果观测点较多，不便于逐点确定，可以按照第2种方法，以方位、 高度角要素确定任意方向待测点大气垂直折光系数计算模型。具体选择那 种方法要根据精度要求和环境特征来确定。

另外，结合本步骤内容再谈谈步骤1中的点位布置，基准点布置应根 据初选的拟合模型，评估由于基准点原始参考高程测量误差对拟合精度的 影响程度，使用中应区别不同特征项目进行分别研究。

6单向测距三角高程的确定

当进行待测点观测时，只要观测出基准方向斜距S_斜、垂直角α，仪器 高I和标高L，便可通过公式(14)计算出测量时刻的基准方向大气垂直折 光系数：

KH＝1-[(h-S_斜×Sinα-I+L)2R/(S_斜×cosα)²](14)

由于在步骤5中已经确定出了基准方向与待测方向折光系数k的相关 关系。将式(14)计算出的系数KH_i代入式(15)：

KB_i＝a+bKH_i(15)

便可求出待测方向折光系数KB_i，将KB_i代入式(16)，就可解算出各个 方向修正后精密测距三角高差值：

h_i＝S_斜×Sinα_i+I_i-L_i+(1-KB_i)×(S_斜×Cosα_i)²/2R(16)

7单向测距三角高程网的定权及平差

如果通过两个或者多个测站对待测点进行测量，可形成附合或者闭合 路线，可依据高程闭合差进一步修正折光系数。计算过程如下：

设通过两个基准点测量的至待测点B1高差分别为h₁和h₂，两个基准点 高程为H₁和H₂，则存在如下关系：

ω＝H₁+h₁-h₂-H₂(17)

式中的ω就是闭合差，理论上ω＝0，然后将ω按照距离确权分配到高 差h₁和h₂上，再将这个数值再代回式(2)进行K值迭代计算。关于平差中 确权的问题，在这里有必要探讨一下：

根据误差传播定律，设h、S、a、k、I及L的中误差分别为m_h、m_s、m_a、 m_k、m_I、m_L，据式(1)按照误差传播定律：

$m_{{}_h}^2={\left(\frac{\partial h}{\partial s}\right)}^2{m}_{{}_s}^2+{\left(\frac{\partial h}{\partial a}\right)}^2{m}_{{}_a}^2+{\left(\frac{\partial h}{\partial K}\right)}^2{m}_k^2+{\left(\frac{\partial h}{\partial I}\right)}^2{m}_{{}_I}^2+{\left(\frac{\partial h}{\partial L}\right)}^2{m}_{{}_L}^2---\left(18\right)$

下边忽略推导过程，直接给出单向高差精度公式：