工程测量中一类参数曲线的拟合

书书书第２７卷第１期２００７年２月大地测量与地球动力学ＪＯＵＲＮＡＬＯＦＧＥＯＤＥＳＹＡＮＤＧＥＯＤＹＮＡＭＩＣＳＶｏｌ．２７Ｎｏ．１　Ｆｅｂ．，２００７　　文章编号：１６７１５９４２（２００７）０１０１０００４工程测量中一类参数曲线的拟合陈基伟（上海市地质调查研究院，上海　２０００７２）摘　要　测量中经常遇到一些特殊参数曲线，已知它的类型，但需要了解它的参数。不同于以往基于最小二乘距离作为拟合函数（不符合度）进行拟合，提出了一种基于欧氏距离的不符合度拟合算法，按照非线性方程解法求取了各项参数，并通过数据进行了验证。结果表明，这种拟合算法容易得到稳定解且具有较好的拟合精度。由于这种拟合函数较易构造，因此对于其它类参数曲线拟合亦便于采用。关键词　工程测量　参数曲线　非线性方程　最小二乘　拟合中图分类号：ＴＢ２２　　　　文献标识码：ＡＡＫＩＮＤＯＦＰＡＲＡＭＥＴＲＩＣＣＵＲＶＥＦＩＴＴＩＮＧＩＮＥＮＧＩＮＥＥＲＩＮＧＳＵＲＶＥＹＣｈｅｎＪｉｗｅｉ（ＳｈａｎｇｈａｉＩｎｓｔｉｔｕｔｅｏｆＧｅｏｌｏｇｉｃａｌＳｕｒｖｅｙ，Ｓｈａｎｇｈａｉ　２０００７２）Ａｂｓｔｒａｃｔ　Ｓｏｍｅｐａｒｔｉｃｕｌａｒｐａｒａｍｅｔｒｉｃｃｕｒｖｅｓａｒｅｕｓｕａｌｌｙｓｅｅｎｉｎｅｎｇｉｎｅｅｒｉｎｇｓｕｒｖｅｙ，ｗｈｏｓｅｔｙｐｅｓａｒｅｋｎｏｗｎｂｕｔｗｈｏｓｅｐａｒａｍｅｔｅｒｓａｒｅｕｎｋｎｏｗｎ．Ｄｉｆｆｅｒｉｎｇｆｒｏｍｔｒａｄｉｔｉｏｎａｌｌｅａｓｔｓｑｕａｒｅｄｉｓｔａｎｃｅｆｉｔｔｉｎｇ，ｔｈｉｓｐａｐｅｒｐｒｅｓｅｎｔｓａｍｏｄｅｌｏｆｐａｒａｍｅｔｒｉｃｃｕｒｖｅｆｉｔｔｉｎｇ，ｗｈｉｃｈｉｓｂａｓｅｄｏｎＥｕｃｌｉｄｄｉｓｔａｎｃｅ．Ｐａｒａｍｅｔｒｉｃｒｏｏｔｓｏｆｔｈｅｓｅｃｕｒｖｅｓａｒｅｄｅｒｉｖｅｄｂｙｅｘｔｒａｃｔｉｎｇｎｏｎｌｉｎｅａｒｅｑｕａｔｉｏｎ．Ｔｈｅｄｉｇｉｔａｌｒｅｓｕｌｔｓｓｈｏｗｔｈａｔｔｈｉｓｍｅｔｈｏｄｃａｎｇｅｔｓｔａｂｌｅｓｏｌｕｔｉｏｎｗｉｔｈｈｉｇｈｐｒｅｃｉｓｉｏｎａｎｄｉｔｃａｎｂｅｅａｓｉｌｙａｄｏｐｔｅｄｆｏｒｏｔｈｅｒｐａｒａｍｅｔｒｉｃｃｕｒｖｅｓｂｅｃａｕｓｅｏｆｉｔｓｓｉｍｐｌｙｄｅｄｕｃｅｄｆｏｒｍｕｌａｓ．Ｋｅｙｗｏｒｄｓ：ｅｎｇｉｎｅｅｒｉｎｇｓｕｒｖｅｙ，ｐａｒａｍｅｔｒｉｃｃｕｒｖｅ，ｎｏｎｌｉｎｅａｒｅｑｕａｔｉｏｎ，ｌｅａｓｔｓｑｕａｒｅ，ｆｉｔｔｉｎｇ１　前言在许多工业部件的设计当中，除了一些我们经常遇到的用代数方程表示的二次曲线和二次曲面类型以外，还有许多曲线方程是用参数方程来表示的。对于某些曲线如螺旋线，当用普通代数方程表示时往往非常繁琐也难以让人理解，而用参数方程来表示往往让人易于理解同时便于进行模拟。由于参数方程利用若干独立参数对直角坐标系的二维或三维坐标进行描述，因此，在拟合它们的参数方程时，首先需要利用最小二乘法构建参数与直角坐标之间的函数关系，才能够利用直角坐标来对这些参数进行求解，从而完成对这些独立参数所表示的参数方程的求解。通常建立平面数据的二次曲线拟合目标函数有两种方法：基于代数距离和基于垂直距离。常用的多项式拟合就是基于代数距离的拟合问题。而基于垂直距离的拟合会产生四次方程或者二次方程组，而四次方程（二次方程组）的求解属于非线性问题，往往使得计算太复杂，在实际中一般不被采用［１］。ＴａｕｂｉｎＧ研究了通过增加二次约束条件，最终可将二次曲线拟合问题转变成条件极值问题，并且得到收稿日期：２００６０６２９作者简介：陈基伟，男，１９６９年生，博士，高级工程师，主要从事大地测量与遥感应用研究．Ｅ－ｍａｉｌ：[e-mail]　第１期陈基伟：工程测量中一类参数曲线的拟合的是非线性方程［２］。对于一般二次曲线，基于欧氏距离的拟合由于没有限定形式的解，同时考虑计算效率等原因，以往研究对基于欧氏距离的最小二乘拟合法研究较少①②。本文提出一种基于欧氏距离的不符合度进行拟合的算法，按照非线性方程解法求取了各项参数。２　参数曲线方程的拟合２．１　方法简介以心形线的拟合为例。如图１，心形线的标准参数方程表示为ｘ′（ｔ）＝ａｃｏｓｔ（１＋ｃｏｓｔ）ｙ′（ｔ）＝ａｓｉｎｔ（１＋ｃｏｓｔ{）（１）其中ａ≠０，ｔ∈［０，２π］。而心形线的代数方程表示为（ｘ′２＋ｙ′２－ａｘ′）２＝ａ２（ｘ′２＋ｙ′２）（２）在心形线边缘上测定了一定数量的点坐标（ｘｉ，ｙｉ，ｚｉ）Ｔ。如果测定的点没有准确位于心形线上，就需要对参数ａ进行估计。根据最小二乘法，若观测了ｎ个测点，第ｉ个测点到心形线的欧氏距离平方可以表示为：Ｄ２（ａ）＝∑Ｎｉ＝１｜（ｘｉ，ｙｉ）－（ａｃｏｓｔ（１＋ｃｏｓｔ），ａｓｉｎｔ（１＋ｃｏｓｔ））｜２＝∑Ｎｉ＝１（ｘｉ－ａｃｏｓｔ（１＋ｃｏｓｔ））２＋∑Ｎｉ＝１（ｙｉ－ａｓｉｎｔ（１＋ｃｏｓｔ））２（３）则Ｄ２（ａ）对ａ的偏导数为１２Ｄ２ａ＝∑ｉ－ｃｏｓｔ（１＋ｃｏｓｔ）［－ａｃｏｓｔ（１＋ｃｏｓｔ）＋ｘｉ］＋∑ｉ－ｓｉｎｔ（１＋ｃｏｓｔ）［－ａｓｉｎｔ（１＋ｃｏｓｔ）＋ｙｉ］（４）由一般函数的泰勒展开可以得出，要使第ｉ个测点到心形线的欧氏距离平方Ｄ２（ａ）趋于０，只有上式的偏导数为０。于是有∑ｉ－ｃｏｓｔ（１＋ｃｏｓｔ）［－ａｃｏｓｔ（１＋ｃｏｓｔ）＋ｘｉ］＋∑ｉ－ｓｉｎｔ（１＋ｃｏｓｔ）［－ａｓｉｎｔ（１＋ｃｏｓｔ）＋ｙｉ］＝０（５）由式（１）得到ｔ＝ａｒｃｔｇｘｉｙｉ（６）将式（５）与式（６）联立，按照单变量函数的标准数值寻根方法（如牛顿法）来迭代解求出，则心形线的方程得到拟合。２．２　各类参数曲线的限制条件方程同心形线的拟合方法相似，可以根据这种测点到参数曲线的欧氏距离平方为最小原则，构造出不同参数曲线的拟合目标函数。１）星形线如图１，星形线的参数方程表示为图１　心形线（ａ）、星形线（ｂ）、圆渐开线（ｃ）、三尖内摆线（ｄ）和螺旋线（ｅ）的拟合Ｆｉｇ．１　Ｆｉｔｔｉｎｇｏｆｈｅａｒｔｓｈａｐｅｌｉｎｅ（ａ），ｓｔａｒｓｈａｐｅｌｉｎｅ（ｂ），ｉｎｖｏｌｕｔｅｏｆｔｈｅｃｉｒｃｌｅ（ｃ），ｔｒｉｃｕｓｐｉｎｎｅｒｃｙｃｌｏｉｄ（ｄ）ａｎｄｈｅｌｉｃａｌｌｉｎｅ（ｅ）ｘ′（φ）＝３ｂｃｏｓφ＋ｂｃｏｓ（３φ）＝４ｂｃｏｓ３φ＝ａｃｏｓ３φｙ′（φ）＝３ｂｓｉｎφ－ｂｓｉｎ（３φ）＝４ｂｓｉｎ３φ＝ａｓｉｎ３{φ（７）其中ａ≠０，φ∈［０，２π］。其代数方程表达为ｘ′２／３＋ｙ′２／３＝ａ２／３（８）可以得出星形线的拟合目标函数为∑ｉｃｏｓ３φ（ｘｉ－ａｃｏｓ３φ）＋∑ｉｓｉｎ３φ（ｙｉ－ａｓｉｎ３φ）＝０（９）其中φ为φ＝ａｒｃｔｇｙｉｘ()ｉ１／３（１０）式（９）和式（１０）联立，按照单变量函数的标准数值寻根牛顿法迭代解求出ａ。２）圆渐开线如图１，圆渐开线的参数方程表示为ｘ′（ｔ）＝ａ（ｃｏｓｔ＋ｔｓｉｎｔ）ｙ′（ｔ）＝ａ（ｓｉｎｔ－ｔｃｏｓｔ{）（１１）上式中ａ≠０。可以得出圆渐开线的拟合目标函数为∑ｉ（ｃｏｓｔ＋ｔｓｉｎｔ）（ａｃｏｓｔ＋ａｔｓｉｎｔ－ｘｉ）＋∑ｉ（ｔｃｏｓｔ－ｓｉｎｔ）（ａｔｃｏｓｔ－ａｓｉｎｔ＋ｙｉ）＝０（１２）其中ｔ＝ｘ２ｉ＋ｙ２ｉ－ａ２ａ槡２（１３）１０１①②ＰｅｔｋｏＦａｂｅｒａｎｄＢｏｂＦｉｓｈｅｒ．ＡＢｕｙｅｒ’ｓｇｕｉｄｅｔｏｅｕｃｌｉｄｅａｎｅｌｌｉｐｔｉｃａｌｃｙｌｉｎｄｒｉｃａｌａｎｄｃｏｎｉｃａｌｓｕｒｆａｃｅｆｉｔｔｉｎｇ［ＥＢ／ＯＬ］．ｗｗｗ．ｂｍｖａ．ａｃ．ｕｋ／ｂｍｖｃ／２００１／ｐａｐｅｒｓ／２１／ａｃｃｅｐｔｅｄ＿２１．ｐｄｆ，２００１：５２１－５３０．ＲｏｂｅｒｔＢＦｉｓｈｅｒ．Ａｐｐｌｙｉｎｇｋｎｏｗｌｅｄｇｅｔｏｒｅｖｅｒｓｅｅｎｇｉｎｅｅｒｉｎｇｐｒｏｂｌｅｍｓ［ＥＢ／ＯＬ］．ｈｏｍｅｐａｇｅｓ．ｉｎｆ．ｅｄ．ａｃ．ｕｋ／ｒｂｆ／ＣＶｏｎｌｉｎｅ／ｌｏｃａｌ＿ｃｏｐｉｅｓ／ｆｉｓｈｅｒ／ｃａｄ０３／ｃｖｏｎ．ｈｔｍｌ，２００１．大地测量与地球动力学２７卷式（１２）与式（１３）联立，利用牛顿法迭代解求出ａ值。３）三尖内摆线如图１，三尖内摆线的参数方程为ｘ′（φ）＝２ｂｃｏｓφ＋ｂｃｏｓ（２φ）ｙ′（φ）＝２ｂｓｉｎφ－ｂｓｉｎ（２φ{）（１４）上式中φ∈［０，２π］，ｂ＞０，ｂ∈Ｒ。可以得出三尖内摆线的拟合目标函数为∑Ｎｉ＝１（２ｃｏｓφ＋ｃｏｓ２φ）（２ｂｃｏｓφ＋ｂｃｏｓ２φ－ｘｉ）＋∑Ｎｉ＝１（２ｓｉｎφ＋ｓｉｎ２φ）（２ｂｓｉｎφ＋ｂｓｉｎ２φ－ｙｉ）＝０（１５）其中φ＝１３ａｒｃｃｏｓｘ２＋ｙ２－５ｂ２４ｂ()２（１６）式（１５）与式（１６）联立，通过牛顿法迭代求解出ｂ值。４）螺旋线如图１，螺旋线的参数方程为ｘ′（ｔ）＝ａｃｏｓ（ωｔ＋φ）ｙ′（ｔ）＝ｂｓｉｎ（ωｔ＋φ）ｚ′（ｔ）＝{ｔ（１７）其中ａ＞０，ｂ＞０，ω＞０，φ∈［０，２π］，ｔ∈Ｒ。如果椭圆的长、短轴ａ，ｂ相等，即成为圆。这时，式（１７）方程退化为一般的数值方程，则Ｄ２（ｒ，ω，φ）对ｒ，ω，φ的偏导数为１２Ｄ２ｒ＝Ｎｒ－∑ｉｘｉｃｏｓψｉ－∑ｉｙｉｓｉｎψｉ１２ｒＤ２ω＝∑ｉｘｉｚｉｓｉｎψｉ－ｙｉｚｉｃｏｓψｉ１２ｒＤ２φ＝∑ｉｘｉｓｉｎψｉ－ｙｉｃｏｓψｉ（１８）其中ψｉ＝ωｚｉ＋φ，根据上式中的第一个方程可以解得ｒｒ＝１Ｎ∑ｉｘｉｃｏｓψｉ＋ｙｉｓｉｎψｉ（１９）通过三角函数和差公式对后二个方程进行化解，再根据二个矢量平行的条件可以得到关于ω的非线性方程；再利用牛顿法迭代解求出ω；求得ω后，再逐项回代可求解出φ和ａ。２．３　曲线误差的评定根据测点坐标和各类曲线的参数方程，可以容易得到测点到曲线的偏差Δｉ以及度量误差ＬＳ。１）心形线通过计算得到心形线上某点的偏差为Δｉ＝１ｘ２ｉ＋ｙ２ｉ（ｘ２ｉ＋ｙ２ｉ）２－２ａ（ｘ２ｉ＋ｙ２ｉ）（ｘｉ＋ｘ２ｉ＋ｙ２槡ｉ([）＋ａ２（ｙ２ｉ＋２ｘｉ（ｘｉ＋ｘ２ｉ＋ｙ２槡ｉ)]））１２（２０）心形线的度量误差表示为ＬＳ＝ｍａｘ（Δｉ）－ｍｉｎ（Δｉ）　ｉ＝１，２，…，ｎ（２１）２）星形线星形线的误差评定公式为Δｉ＝（ｘｉ－ａｃｏｓ（ｙｉｘｉ）１／()３２＋（ｙｉ－ａｓｉｎ（ｙｉｘｉ）１／()３槡２ｉ＝１，２，…，ｎ（２２）３）圆渐开线、三尖内摆线根据这两类曲线的参数方程，可以很容易得到它们的误差评定公式。４）螺旋线很容易得到某点到圆螺旋线的偏差为Δｉ＝（ｘｉ－ａｃｏｓ（ωｚｉ））２＋（ｙｉ－ｂｓｉｎ（ωｚｉ））槡２ｉ＝１，２，…，ｎ（２３）螺旋线的度量误差如式（２１）。３　算例为验证本文介绍的方法，对某个螺旋线测量了一些离散点，其坐标及偏差（由式（２３）计算得到）如表１所示。由式（１７）计算ω的初值为２．０２２，ａ的初值为０．２５３。利用Ｍａｔｈｅｍａｔｉｃａ数值分析软件中ＦｉｎｄＲｏｏｔ函数对ω，ａ进行迭代计算，得到稳定迭代结果ω＝表１　一组离散点的坐标Ｔａｂ．１　Ｃｏｏｒｄｉｎａｔｅｓｏｆａｓｅｔｏｆｄｉｓｃｒｅｔｅｐｏｉｎｔｓ离散点的坐标ｘ３．９９１７３．９３００３．８０５１３．６１９５３．３７７０３．０８１１２．７３９０２．３５２０１．９２７５ｙ０．２５１７０．７４９０１．２３５６１．７０２９２．１４２４２．５４９９２．９１６７３．２３６３３．５０４６ｚ０．０３１１０．０９３８０．１５７７０．２１９７０．２８３７０．３４５３０．４０７５０．４７０９０．５３３６Δ０．００２９０．００３２０．００５７０．００１８０．００８４０．００２６０．００６９０．００２４０．００３２ｘ１．４７２３０．９９４４０．５０１０－０．０００２－０．５００８－０．９９５６－１．４７１５－１．９２６１－２．３５１９ｙ３．７１９１３．８７３６３．９６９３３．９９９６３．９６９１３．８７５０３．７１９１３．５０５９３．２３６０ｚ０．５９７７０．６６０１０．７２３１０．７８６３０．８４９００．９１１５０．９７４０１．０３６２１．１００４Δ０．００６２０．００３１０．００３８０．００７１０．００６８０．００３１０．００１５０．００３