非线性最小二乘法在工程测量坐标系统转换中的应用_刘洪志.pdf

第 46 卷 第 4 期2023 年 4 月测绘与空间地理信息GEOMATICS&SPATIAL INFORMATION TECHNOLOGYVol.46,No.4Apr.,2023收稿日期:2022-07-25作者简介:刘洪志(1986-),男,黑龙江哈尔滨人,工程师,硕士,2014 年毕业于电子科技大学软件工程专业,主要从事工程测量工作。非线性最小二乘法在工程测量坐标系统转换中的应用刘洪志1,刘莉鑫2(1.中通大地空间信息技术股份有限公司哈尔滨分公司,黑龙江 哈尔滨 150028;2.黑龙江省测绘地理信息学会,黑龙江 哈尔滨 150081)摘要:针对大型工程测量项目,工程坐标系直接转换为 WGS-84 或 CGCS2000 坐标系的过程中,传统方法效率不高,而且经过多次旋转造成误差传递等问题,本文提出一种基于非线性最小二乘法的算法,可避免误差的传递,减少循环次数,保证转换精度,提高运行速度。首先介绍在大型工程测量项目中坐标系的转换场景,然后介绍传统的坐标系转换方法,并分析其利弊,最后论述利用非线性最小二乘法解决工程测量坐标系转换的可行性,并模拟转换数据,避免数据发散。关键词:最小二乘法;工程测量;坐标系转换中图分类号:P282.2 文献标识码:A 文章编号:1672-5867(2023)04-0186-03Application of Nonlinear Least Square Method in the Transformation of Engineering Surveying Coordinate SystemsLIU Hongzhi1,LIU Lixin2(1.Zhongtong Dadi Space Information Technology Co.,Ltd.,Harbin Branch,Harbin 150028,China;2.Heilongjiang Provincial Society of Surveying and Mapping Geographic Information,Harbin 150081,China)Abstract:In the process of directly transforming engineering coordinate systems into WGS-84 or CGCS2000 coordinate system in large-scale engineering measurement projects,the efficiency of traditional methods is not high,and error propagation is caused by multiple rotations.This paper proposes an algorithm based on nonlinear least square method,which avoids error propagation,reduces the num-ber of cycles,ensures the transformation accuracy,and improves the running speed.This paper first introduces the transformation scene of coordinate system in large-scale engineering survey projects,then introduces the traditional coordinate system transformation method,analyzes its advantages and disadvantages,and finally discusses the feasibility of using nonlinear least square method to solve the transformation of engineering survey coordinate systems,and simulates the transformed data to avoid data divergence.Key words:least square method;engineering survey;coordinate system transformation0 引 言大型工程测量项目,具有一定的特殊性,在测量过程中可能涉及多个坐标系统,有世界级的坐标系统,也有国家级的坐标系统,还有地方坐标系统和工程独立的坐标系统,这些坐标系统各有特点,在不同的场景下发挥着各自的作用1-2。但从项目全局考虑,设计、GPS 定位、测量和放样等工作最好在同一坐标系统下开展,因此需要在各个坐标系统之间进行转换3。1 传统坐标系转换方法工程测量项目中经常用到的坐标系统有北京 54 坐标系、西安 80 坐标系、WGS-84 坐标系、CGCS2000 坐标系、城市坐标系、独立坐标系等4,各坐标系之间的转换方法基本相似,如可以利用 3 个公共点求得转换参数的方法,建立数学模型后,利用算法进行坐标系转换,坐标系之间的转换还要满足一定的精度需求5-6。基准点转换是将一个坐标系统中的点坐标转换为另一个坐标系统中对应的点坐标,再利用 GPS 定位得到的地心坐标,通过基准点转换,可以得到对应的地面坐标。地面坐标系是在一定区域范围内,利用最小二乘法,获得最佳拟合椭球,椭球中心看作坐标系统的原点,投影后得到的平面坐标。一般地,使用 GPS,联测 3 个地面点,即可计算得到地面坐标系统的坐标7。大型工程测量项目坐标转换技术路线如图 1 所示。图 1 大型工程测量项目坐标转换技术路线Fig.1 Technical route of coordinate transformation for large-scale engineering survey projects在实际应用过程中,从工程坐标系统转换到世界坐标系统,往往要经历 1 次二维三参数转换到区域坐标系统,再经历 2 次三维七参数转换才能转换到 WGS-84 或CGCS2000 坐标系统,中间要进行多次旋转,导致误差的传递,造成精度不理想8。工程坐标系转换为 WGS-84或 CGCS2000 坐标系统流程如图 2 所示。图 2 工程坐标系转换为 WGS-84 或 CGCS2000 坐标系统流程Fig.2 Process of transforming engineering coordinate system into WGS-84 or CGCS2000 coordinate system2 非线性最小二乘法利用非线性最小二乘法转换坐标系统的步骤,一般是从一组选定的初始值开始,按照一定的规律沿二次型VTPV 减小而变化,得到区域空间中使得 VTPV 较小的一个点的坐标,然后再以这个点作为新起点,进行下一次循环,直到 VTPV 不再减少为止9。设空间中直角坐标系统的转换函数如式(1)所示。y=f(t;x,y,z,x,y,z,m)(1)式中,t 为自变量;x、y、z 为平移参数;x、y、z为 x、y、z 方向上的旋转欧拉角;m 为尺度比变化参数。由于函数 y=f 的参数是非线性的,所以用最小二乘法确定七参数的数学模型,就是求如式(2)所示的极小值。(x,y,z,x,y,z,m)=ki=1fti3 1;x,y,z,x,y,z,m()-yi|2(2)由式(2)可知,当 k 3 时,可以建立如式(3)所示的方程组,式(3)中 k 3,1 i k。U(x)=i1(x,y,z,x,y,z,m)=(1+m)R11xi+R12yi+R13zi+z-xii2(x,y,z,x,y,z,m)=(1+m)R21xi+R22yi+R23zi+z-yii3(x,y,z,x,y,z,m)=(1+m)R31xi+R32yi+R33zi+z-zik1(x,y,z,x,y,z,m)=(1+m)R11xk+R12yk+R13zk+z-xkk2(x,y,z,x,y,z,m)=(1+m)R21xk+R22yk+R23zk+z-ykk3(x,y,z,x,y,z,m)=(1+m)R31xk+R32yk+R33zk+z-zk|=0(3)设 U:D R7 R3k,U=(u1,um)T,式(3)可以简化为式(4)。(x)=12U(x)TPU(x)(4)式中,:D R,因此方程组的最小二乘法得出的解就是(x)的极小值,如式(5)所示。(x)=min(x)=12minU(x)TU(x)(5)设(x)梯度函数为 g(x)=(x)=0,多元函数极值存在的必要条件是 U 在上 D 可微,可以得到式(6)。g(x)=(x)=U(x)TxU(x)=0(6)令 G(x)=U(x)TxU(x)x,则式(6)可以写作如式(7)所示。xk+1=xk-G(xk)-1g(xk)(7)对于非线性函数 U(x),G(x)为对称非负矩阵,为了避免 G(xk)发散,可以设G(xk)=G(xk)+kI,当 kI 0时,G(xk)是对称正矩阵,因此得出式(8)。pk(k)=-G(xk)-1(xk)=-G(xk)+kI-1(xk)(8)可以推导出循环算法公式如式(9)所示。xk+1=xk-G(xk)-1(xk)=xk-G(xk)+kI-1(xk)(9)当 kI 0 时,(xk+1)0,0=10-3,=2,0 k2计算 U(xk),U(xk)xk,(xk),G(xk),(xk),0 j3解 G(xk)+kIpk=-(xk)方程组,求 pk(k)=pk4计算 xk+1=xk+pk(k)和(xk+1)5如果(xk+1)(xk),并且 j=0,取 k=k,1j,跳转到步骤 3,如果 j 0,跳转到步骤 76如果(xk+1)(xk),取 k=k,1j,跳转到步骤 37如果满足 pk ,则 xk+1为极小点 x的近似值,此时循环终止,否则将 xk+1 xk,k+1 k,跳转到步骤 2781第 4 期刘洪志等:非线性最小二乘法在工程测量坐标系统转换中的应用3 实例分析在某大型工程测量项目中,从工程独立坐标系统转换为 WGS-84 坐标系或 CGCS2000 坐标系,如果旋转角过大,则传统的线性坐标转换模型会导致误差传递。为了避免在一些工程项目中采用高斯-牛顿法造成的矩阵发散情况,本文采用非线性最小二乘法将工程坐标系直接转换为 WGS-84 或 CGCS2000 坐标,能够有效避免误差的传递,具体实例分析如下。利用 Matlab 编写算法程序,设定模拟初始值,x、y、z方向上的旋转欧拉角x,y,z=45,60,75,平移参数x,y,z=10 000 000,200 000,10,尺度比变化参数 m=0.9,根据七参数模型,得出转换后坐标,模拟计算数据见表 2。从表 2 可以看出,通过 1 次循环,即可得到近似的结果,4 次循环后,得到收敛的结果,程序效率较高。误差分析见表 3。表 2 计算结果Tab.2 Calculation results循环次数xyzmxmymzm147.608 860 357 073.431 638 6887.719 044 382 00.832 555 672 01 000 703.860 28199 730.327 1280.874 69244.596 638 845 460.144 195 5875.515 337 793 80.899 846 250999 996.280 88200 004.677 522.950 267345.000 000 008 960.000 000 2374.999 999 980 10.900 000 004 1999 999.999 95200 000.000 019.999 942445.000 000 000 060.000 000 0075.000 000 000 00.900 000 000 01