机场道面高程图（也称为板角高程图）测量一般在机场改扩建设计前或者机场建设竣工后进行。板角高程测量精度至毫米^[1]，一般尾数为估读，且板角点相对临近水准点高程中误差不大于5 mm，极限误差不大于1 cm，因此可认为，板角高程测量的精度要求优于1 cm。传统的测量方法将平面数据和高程数据分开进行采集，然后一一对应合并（一个机场大概3万个数据），过程繁琐，工作效率低下，且很容易出现错误。地面三维激光扫描仪（terrestrial laser scan⁃ ner，TLS）发展已经30多年，特别是近10年来，精度高、测速快的扫描仪已经逐渐走进各个应用领域^[2]，比如工程测量^[3]、文物保护^[4]、工业测量^[5]、交通道路测量^[6-8]等，极大提高了工作效率。

但根据机场道面高程测量中的精度要求，采用常规扫描仪流程，主要存在以下3个主要问题：一是靶标控制点的高程精度要足够高，且测量过程中需要尽量减少“量高误差”；二是由于道面测量对平面精度和高程精度要求不一致，特别是扫描仪没有整平装置时，平面坐标误差会影响高程匹配的精度，且配套软件没有对权值进行设定；三是由于机场周边地势平坦，且无连续制高点，距离仪器越远，道面测量点的反射角就越小，大于70 m几乎就没有回波数据，这样就使得相邻测站的距离不能太大。因此采用常规的三维扫描仪测量方法，不仅效率低，而且也不能满足厘米级以内的高程精度。

本文首先分析了使用加权整体最小二乘^[9-11]进行坐标匹配转换的算法，该算法公式经过推导和向量化后，采用LM（Levenberg⁃Marquardt）算法求出最优转换参数，并对目标坐标和测站坐标的误差来源进行分析，确定了权值；其次，为了减少靶标量高误差和提高地面回波率，设计了一种控制点布设方案和测量方法；最后，通过实例数据的分析和比较，得出结论。

1 坐标匹配算法和权值的确定 1.1 加权整体最小二乘法坐标匹配算法

TLS外业测量时采用的是测站坐标系统，而目标坐标为测图时采用的坐标系统，两者之间需要转换。由于测站坐标和目标坐标均含有随机误差，直接采用经典最小二乘法求解转换参数时，求出的结果会有偏差^[12]。因此需要使用变量中存在误差的EIV（errors⁃in variables）模型，求出最优估计参数解。从测站坐标到目标坐标系统EIV模型公式为：

$ \left[\begin{array}{l} X \\ Y \\ Z \end{array}\right]+\left[\begin{array}{c} V_{X} \\ V_{Y} \\ V_{Z} \end{array}\right]=\lambda \boldsymbol M\left[\begin{array}{c} x+V_{x} \\ y+V_{y} \\ z+V_{z} \end{array}\right]+\left[\begin{array}{c} \Delta x \\ \Delta y \\ \Delta z \end{array}\right] $

(1)

式中, $\left[\begin{array}{lll}X \mathrm{YZ}\end{array}\right]^{\mathrm{T}}$为目标坐标向量; [$\left.\begin{array}{lll}V_{x} & V_{y} & V_{z}\end{array}\right]^{\mathrm{T}}$为目标坐标误差向量; $[x y z]^{\mathrm{T}}$表示测站坐标向量; $\left[\begin{array}{lll}v_{x} & v_{y} & v_{z}\end{array}\right]^{\mathrm{T}}$为测站坐标误差向量, $[\Delta x \Delta y \Delta z]^{\mathrm{T}}$为平移参数向量; $\lambda$代表尺度缩放参数, TLS的坐标转换中, 一般认为$\lambda$等于$1 ; \boldsymbol M$为3个旋转角组成的旋转矩阵。

$ \boldsymbol M=\boldsymbol R\left(\varepsilon_{z}\right) \boldsymbol R\left(\varepsilon_{y}\right) \boldsymbol R\left(\varepsilon_{x}\right) $

(2)

根据式（1）、式（2）可以得出：

$ \begin{aligned} Z+&V_{Z}=\sin \varepsilon_{y}\left(x+V_{x}\right)-\cos \varepsilon_{y} \sin \varepsilon_{x}\left(y+V_{y}\right)+\\ &\cos \varepsilon_{x} \cos \varepsilon_{y}\left(z+V_{z}\right)+\Delta z \end{aligned} $

(3)

由式(3) 可以看出, 在仪器整平的情况下, $\varepsilon_{x} 、\varepsilon_{y}$为小角度, $\sin \varepsilon_{x} 、\sin \varepsilon_{y}$也为极小值, $\cos \varepsilon_{x} 、\cos \varepsilon_{y}$接近1, 这时, 目标坐标高程主要与测站坐标的$z$方向坐标分量、误差、平移参数相关。但是, TLS一般没有整平装置, 这时, 目标坐标高程与测站3个坐标分量以及误差、平移参数均相关。

由式（1）、式（2）可以看出旋转参数和系数矩阵呈非线性关系，无法直接使用整体最小二乘算法求解，一般将式（1）转换成线性化的模型求解，但该方法仅仅适合于两套坐标系间的旋转角较小时的转换，若旋转角较大时，模型误差会很大，甚至完全得不到正确的结果^{[13, 14]}。

$ \left[\begin{array}{c} x \\ y \\ z \end{array}\right]+\left[\begin{array}{c} V_{x} \\ V_{y} \\ V_{z} \end{array}\right]=\left[\begin{array}{l} x+V_{x} \\ y+V_{y} \\ z+V_{z} \end{array}\right] $

(4)

将式（1）和式（4）联立，坐标转换模型转化为非线性GM模型：

$ \boldsymbol l+\boldsymbol V=f(\beta) $

(5)

式中, $\boldsymbol l$为式(1) 和式(4) 中目标坐标和测站坐标的向量化; $\boldsymbol V$为误差的向量化。假设$A$为测站坐标系统的$n$个坐标, $\beta_{c}$为6个坐标转换参数, vec表示逐行向量化, 则变量参数$\boldsymbol \beta$可表示为:

$ \boldsymbol \beta=\operatorname{vec}\left[\begin{array}{l} \beta_{c} \\ \operatorname{vec}(A) \end{array}\right] $

(6)

其中，f(β)为式（1）和式（4）联立后的方程的右半部分。

经过向量化变化后，式（5）转换为非线性方程的最优化问题，加权整体最小二乘的最优估计的函数可表示为：

$ \boldsymbol \beta=\arg \min ( F(\boldsymbol \beta)) $

(7)

$ F(\boldsymbol \beta)=\boldsymbol V^{\mathrm{T}} \boldsymbol P \boldsymbol V=(f(\boldsymbol \beta)-l)^{\mathrm{T}} \boldsymbol{P}(f(\boldsymbol \beta)-l) $

(8)

式中，P表示权值。

对于非线性方程的最优化问题，常见的方法有梯度法和高斯牛顿法。梯度法主要从负梯度方向进行迭代搜索，并最终找到极小值，该方法开始时下降较快，但靠近极小值时收敛速度很慢。高斯牛顿法中由雅克比（Jacobian）矩阵近似计算海森（Hessian）矩阵，避免了直接计算海森矩阵的复杂性，收敛速度快，但是该算法对初始值较敏感，容易陷入局部极小值。LM算法引入阻尼因子^[15]，每次迭代时调节阻尼因子μ，使算法介于梯度法和高斯牛顿法之间，能有效克服常见缺点。对于本文的加权非线性方程的求解中，在梯度向量G和海森矩阵H更新时需要加入权值项。

每一次迭代时，h表示β值的更新值：

$ \left.\boldsymbol h=\beta_{i+1}-\beta_{i}=-(\boldsymbol H+\mu I)^{-1} \boldsymbol G\right) $

(9)

当梯度的无穷范数小于预定的一个极小值ε1，或者通过更新值h和变量参数的2⁃范数进行比较，当小于变量参数的一个极小比例ε2时，程序找到非线性最小二乘的最优解β，由式（6）可知，解的前6个值就是旋转角和平移量参数。本文中ε1、ε2均设置为1×10^-10。

1.2 权值的确定

机场道面很平整，横坡和纵坡最大也不超过2%^[16]，即平面偏差5 cm，对高程影响也就1 mm。为了保证高精度的道面高程测量成果，同时又提高效率。靶标控制点的平面坐标使用GPS⁃RTK采集，高程则使用二等水准进行测量。道面测量采用RIEGL VZ ⁃ 1000三维激光扫描仪（角度分辨率1. 8″，小于100 m距离的测距精度5 mm），是一款高密度、高精度、高效率的扫描设备。

目标坐标系统下靶标点的测量误差分为平面和高程两个部分。由于平面坐标由GPS⁃RTK采集，有固定误差，对中杆手握对中以及其他因素引起的偶然误差，要准确计算每一个点的平面误差比较困难。由于机场环境空旷，高度角大于15°的卫星一般都能接收到信号，本文计算中统一设置为2. 1 cm中误差（两个坐标分量都是1. 5 cm的中误差）。由于电子水准仪的普遍应用，水准测量的效率和精度都得到了很大提高，有的仪器甚至已经达到了0. 2 mm/km的高程中误差，本案例的控制点采用二等水准测量，平差后高程中误差为0. 6 mm，加上放置靶标的误差后，统一设置为2 mm。

测站坐标系统是扫描仪外业测量时采用的坐标系统，靶标位置在该系统的误差包括两部分：一是仪器的测量误差，另一个是靶标识别的误差。其中在测站坐标系统下，如图 1所示，仪器的观测量为：水平角α、垂直角θ、以及斜距S。则仪器的测量中误差在3个方向上分别为^[17]：

$ \sigma_{x 1}=\pm \sqrt{(\cos \theta \cos \alpha)^{2} \sigma_{s}{ }^{2}+(S \sin \theta \cos \alpha)^{2} \sigma_{a}{ }^{2}+(S \sin \theta \cos \alpha)^{2} \sigma_{\theta}{ }^{2}} $

(10)

$ \sigma_{y 1}=\pm \sqrt{(\cos \theta \sin \alpha)^{2} \sigma_{s}{ }^{2}+(S \cos \theta \cos \alpha)^{2} \sigma_{\alpha}{ }^{2}+(S \sin \theta \sin \alpha)^{2} \sigma_{\theta}{ }^{2}} $

(11)

$ \sigma_{z 1}=\pm \sqrt{(\sin \theta)^{2} \sigma_{s}{ }^{2}+(S \cos \theta)^{2} \sigma_{\theta}{ }^{2}} $

(12)

靶标识别误差主要和靶标上的点云密度有关，扫描仪道面扫描时使用的是6 cm×6 cm的平面靶标，放置在10 m位置时约有6 000个点，点间距约为0. 8 mm；放置在120 m位置时就只有200个点，点间距约为4. 2 mm。虽然扫描仪没有对中和整平装置，但是一般情况下，扫描仪在测量时基本保持水平，本文中将靶标识别误差设定为点间距σ_xy ≈ σ_z2 ≈ ±d，其中x、y代表水平方向，即使用随机软件RiSCAN中的点云个数求出各靶标的点间距后，将水平方向的误差、垂直方向的误差设定为点间距。

如图 2所示，x、y、z方向的识别误差分别为:

$ \sigma_{x 2} \approx \pm d \cdot \cos (\alpha-\pi / 2) $

(13)

$ \sigma_{y 2} \approx \pm d \cdot \sin (\alpha-\pi / 2) $

(14)

$ \sigma_{z 2} \approx \pm d $

(15)

测站坐标的3个方向总误差分别为：

$ \sigma_{i}{ }^{2}=\sigma_{i 1}{ }^{2}+\sigma_{i 2}{ }^{2} \quad i=x, y, z $

(16)

目标坐标和测站坐标的误差都计算出来后, 就可以对整体最小二乘的权值进行确定, $p_{i}=\sigma_{0}{ }^{2} / \sigma_{i}{ }^{2}$, 其中$\sigma_{0}{ }^{2}$可以选定为$\sigma_{i}{ }^{2}$的最大值。权矩阵可表示为:

$ \boldsymbol P=\operatorname{diag}\left(p_{i}\right) $

(17)

式中，diag代表生成对角化矩阵的符号，即将向量的成员作为矩阵对角线上的元素。这样，权值即可确定。

2 实例数据采集和分析 2.1 数据的获取

机场道面周边无连续制高点，数据采集时使用专用脚架提高仪器的架设高度。普通脚架只能把仪器高度架设到1. 5~1. 7 m，使用工业级的专用脚架，把仪器高度架设到2. 6~3. 0 m。提高仪器高度，不仅可以测量得更远，提高效率；同时，由于反射角的增大，反射率增加，也提高了测量精度。

如图 3所示，将靶标贴在平整长方体一面，直接将靶标的投影中心对准控制点平面中心，放置在控制点的最高处，靶标面朝向仪器中心，优点在于避免了过大的对中误差和量高误差。需要注意的是，匹配该靶标点的高程应该为控制点高程加上3 cm（平面靶标的一半尺寸），平面坐标不变。

靶标布设时要分布均匀，且避免靶标点共线和共面的情况，防止坐标匹配时无穷多解的情况。待匹配的靶标点布设在道面两边线外的助航灯基座螺母上，航灯的间距约为50 m，跑道宽度一般和机场等级有关，多数为45~60 m，如图 4所示。一个测站可以扫描12个靶标，其中8个在本测站扫描范围内，另外4个在范围外附近。机场道面测量时，将扫描仪架设在跑道中心线上，为便于相邻测站的结果进行检核，使后一测站的点云和前一测站有重叠，使测站重叠区域离相邻测站的最大距离在100 m左右。

2.2 分析和比较

实验1：采用8种坐标匹配方案，比较各种方案的匹配精度，具体匹配方案和匹配精度如表 1、表 2所示。

在表 2中，方案1只有必要的3个点参与参数计算，没有多余点，无法平差，虽然内符合精度很高，但从检核点的精度统计可以看出高程的差值范围很大，标准差最大，证实了采用本文算法的必要性。

由于EIV模型把误差作为变量进行平差，考虑了目标坐标系统中存在误差，所以在参与计算点的内符合精度统计中，方案3和4的高程差值标准差比较大，但是检核点的高程标准差最小。方案3的检核点高程精度最高，平面精度符合实际RTK测量的情况，因此8个临近的点的参数估计最稳定。方案4的整体精度和方案3相近，也符合式（3）的ε_x、ε_y为小角度的情况。

方案5~方案8为匹配点在一侧的情况，另一侧的检核点高程出现较大的偏差，这是由于靶标分布不均匀造成的，这个结论也符合文献[18]的结论。方案6、方案8使用了两个近130 m的靶标，靶标的分辨率过低引起这两个方案的高程整体精度不如方案5和方案7。

实验2：TLS多数没有整平和对中系统，仪器的坐标平面和目标坐标平面产生夹角在所难免，主要体现在以坐标x轴和y轴发生旋转，本文假设产生0~10°的旋转，使用测站周边的8个点作为参数计算点，剩下的4个点作为检核点，对不加权和加权进行比较。

如图 5所示，等权的情况下，参与计算点的内符合高程转换差值的标准差在1. 3~1. 8 mm之间，检核点的高程转换差值的标准差在3. 5~5. 5 mm之间。不论是内符合还是检核点的精度，都随着旋转角的加大而增加，最大的情况出现在x轴和y轴同时旋转10°的情况。可以看出，在2°旋转角的范围内，高程的转换误差还算稳定，也就是旋转角度较小的情况。

如图 6所示，加权的情况下，参与计算点的高程差值的标准差非常平稳，保持在1. 3 mm，检核点的转换差值标准差也较稳定，保持在3. 5~3. 7 mm。对比图 5、图 6可以看出，符合式（3）的推导，即没有整平装置时，目标坐标高程和测站3个坐标分量以及误差、平移参数都相关的结论。加权的情况下保持着等权差值的最小值，并非常稳定，说明了该算法具有通用性和强壮性。

由于两站重叠部分的点云距离两测站都最远，该处的扫描精度最低，外业时使用常规的方法（水准仪）抽样测量了几个测站重叠部分的高程，共计119个点，与TLS点云过滤分类后的结果进行比较得出，差值范围在-12. 1~8. 7 mm之间，大于1 cm差值的点只有3个，且都位于道面边缘处，是源于该处有覆土或者植被的影响，所有差值的均值为2. 2 mm，差值标准差为±4. 6 mm，也就是精度优于1 cm，可以认为本文的方法完全能满足机场道面测量的要求。

3 结束语

本文分析了加权整体最小二乘算法的理论和实现方法，并且对TLS的目标坐标系统和测站坐标系统误差进行了分析确定，并针对某机场实测的数据进行计算比较，得出结论：由于测量效率的需要，匹配点的目标坐标采用不同的设备采集，这样就需要使用加权整体最小二乘法进行转换参数的求解，这样才能得到最优的参数估计，并有效减小三维激光扫描仪不能整平带来的误差。同时在高精度的道面高程测量的应用中，采用此方法可以有效的提高高程精度，最终达到优于1 cm的高程精度要求。

传统方式测量一个机场的板角高程需要4人（一名观测员、一名记录员、两名跑尺员），15 d左右才能完成工作。采用本文方法只需要3人（两名仪器操作员，一名靶标摆设人员），3 d就可完成同样工作。无论从人力需求上还是工作效率上都明显得到了提高，而且避免了传统方式平面和高程分开采集，最后融合易产生的错误。本文方法对于其他一些对单一方向有高精度要求的测量工程，如道路、特殊构件的平整度测量等，都有很好的借鉴作用。