提取空间拟合椭圆模型的中心坐标算法

提取空间拟合椭圆模型的中心坐标算法范昌胜-李斌'，陈新庄2(1.陕西工商职业学院工程与建筑学院，陕西西安710119；2,西北工业大学理学院应用数学系，陕西西安 710072)

摘要：采用古塔每一层8个观测点的坐标，给出了解决各层中心位置坐标的通用方法.先用一个平面对原始观测数

据进行空间平面拟合，把原本三维的问题降低到二维来考虑，使得原本的空间曲面问题简化成了平面问题.然后在该

平面上建立局部坐标系，通过坐标变换得到观测点在局部坐标系下的坐标.再利用椭圆曲线拟合这8个观测点，得到

椭圆的方程，椭圆的中心即为该塔层时中心.最后利用坐标逆变换得到在原始坐标系下的中心坐标.仿真实验和实验

表明，该算法可以用来解决空间建筑物的中心位置进行标定的问题.算法比其它的算法具有速度快和推广价值高等优 点.关键词：坐标变换;拟合方程；最小二乘法；空间椭圆曲线；中心提取中图分类号:TP391.9

文献标志码:A

文章编号:2095^271(2020)01-0041-08Algorithm for extracting center position of fitted space ellipse modelFAN Chang - sheng1, LI Bin* ,CHEN Xin - zhuang2(1. School of Engineering and Architecture, Shaanxi Vocational and Technical College, Xi'an 710119, P. R. C.；2. Department of Applied Mathematics, School of Science, Northwestern Polytechnic University, Xi'an 710072, P. R. C.)Abstract: The coordinates of 8 observation points in each layer of the ancient tower are used to give a general method for solving

the coordinates of the center position of each layer. Firstly, the spatial plane fitting of the original observation data is performed

by a plane, and the original three - dimensional problem is reduced to two - dimensional, so that the original spatial surface

problem is simplified into a planar problem. Then, a local coordinate system is established on the plane, and coordinates of the observation point in the local coordinate system are obtained by coordinate transformation. Then the elliptic curve is used to fit the eight observation points to obtain the elliptic equation. The center of the ellipse is the center of the tower layer. Finally, the coordinate of the original coordinate system is obtained by inverse coordinate transformation. Simulation experiments and experi­

ments show that the algorithm can be used to solve the problem of calibration of the center position of space buildings. The algo

rithm of this paper has the advantages of faster speed and higher promotion value than other algorithms.Key words: coordinate transformation ； fitting equation ； least squares method ； spatial elliptic curve ； center extraction1研究现状及分析作为古代高层建筑代表之一的砖石古塔,为研究古代建筑技术和古建筑史提供了丰富的物质素材•受古代收稿日期:2018-06-10作者简介:范昌胜(1979-)，男，山东汶上人，讲师，研究方向:最优化方法及科学工程计算基金项目：陕西省教育厅2017年专项科学研究计划课题(17JK0059):陕西省教育厅2017年重点科学研究计划(17JZ011)

42西南民族大学学报（自然科学版）第46卷建筑技术水平的限制以及地质条件恶化、载荷、风力等外力作用下，大多数砖石古塔均存在不同程度的结构损 坏和倾斜、弯扭变形问题，“十塔九斜”已成为公认的事实⑴.为保护古塔，首先需要对古塔的变形形式和变形量

进行观测,为古塔的纠偏等保护措施提供依据•其中，确定古塔各层的中心位置是研究的重要课题⑵.由于古塔 的内底层与顶层之间通常无通视条件，而塔内各层的中心点坐标也无法直接获取，目前只能采用间接方法测量 其中心点坐标卩7 •通常会在古塔的每一层选取了若干个参考点，根据这些观测数据，再确定古塔各层中心位 置.当前经常用的椭圆检测预定位算法中，国内外比较被认可的有:胡志晓⑶使用了间接测量古塔底层中心和

顶层中心的方法，但是未对坐标体系下曲线和曲面的进行转换，中心位置点的坐标存在准确性不足等缺陷.闫 蓿等⑷人采用随机理论的思想，先随机选取6个点拟合椭圆，然后计算与此椭圆匹配的所有样本点个数.重复

此过程一定次数，采用投票机制，匹配样本点多的椭圆即为最优椭圆，构造了一种快速准确剔除误差较大样本

点的改进椭圆拟合算法•李良福㈤等提出了一种有效的基于随机Hough变换的椭圆检测算法（RED）,实现了利 用最小二乘法进行椭圆拟合的检测方法，解决了二维对图像中的多个椭圆进行实时而且高精度检测的难题.安

新源等人e刃提出了一种非线性最小二乘来拟合椭圆，用来提取人眼图像中瞳孔技术，但在参数计算较为复杂

和运算效率不高.ZHOU J,SHIZ, FU Q［9］提出了一种使用射击和弹跳射线技术,提取复杂目标的三维（3 -D） 中心模型的技术•首先基于单视ISAR算法生成目标的3-D逆合成孔径雷达图像•然后使用图像处理算法从三

维ISAR图像中提取三维中心位置.Kyung-Tae Kim,Hyo-Tae Kim［10］提出了一种使用雷达数据来估计二维图

形中心的新技术•该技术利用多个模块网络优化来找到目标中心的二维位置•该方法可以有效地估算目标上的 二维图形的中心.Nair等提出的基于Hough变换的椭圆检测方法以及对它的改进⑴］.这种算法在极坐标下进行

Hough,能在复杂的环境下提取椭圆,但复杂度大,运算效率低氏-131.问题的本质在于如何利用已知的m个参考点的坐标求出古塔中心位置的坐标.由于塔的截面外围形状一

般为正多边形，参考点的位置一般取在正多边形的顶点位置附近，故而存在一个外接圆同时经过这m个参考 点.可以考虑用一个空间圆曲线来拟合这m个参考点,从而获得中心点坐标.事实上，由于塔身受到外载荷，塔

的截面外围形状往往发生变形，截面外接圆也压缩成一个椭圆，因此需要用一个空间椭圆曲线来拟合实际测量 到的m个参考点.难点就在于如何建立空间椭圆曲线的拟合方程，以及如何求解拟合参数.由上述分析，我们设计了古塔中心提取算法.2基于坐标变换的空间椭圆模型中心提取算法2.1基于最小二乘法的空间平面拟合方法考虑一般情况下，假设原始观测点为= 1,2,…,m），用一空间平面眄对原始观测点P,进行最小二乘 拟合•不妨设平面叭的方程为⑶:Ax + By + Cz + 1 = 0 .设原始观测点的坐标为x = （%! ,%2,…，％）'，y =（X』2,…，九），，Z =（Z1 ,Z2,…,Z”），，则该拟合表示成

如下矩阵形式：%、(x y z)B=1)一 1JC丿两边同时左乘（％ y z）7■，得:< -1 >Zy-、X xiziB =-工久工z,丿AC第1期求得拟合系数A、B、C：范昌胜，等:提取空间拟合椭圆模型的中心坐标算法43(Z^Y'X X\£ y^iB =、X xiz<2.2空间平面法向量及空间坐标变换沪丿-Zr,平面771的法向量为n =(A,B,C)r，平面77］与平面7T2之间的夹角”2为法向量n =(A,B,C)r与m =CO,。」)7\"之间的夹角，故n |C|6 = arccos ― ■ ” .y/A2 +B2 +C2旋转原始坐标系O-^z，使得z轴与平面眄的法向量；重合•此坐标变换可以看作两次旋转轴变换，第一次绕%轴旋转0,=arccos系的坐标变换门为:CVb2 + C2，第二次绕y轴旋转e2=arctan —,

Vb2 + c2，从 0 - xyz 系变换到 Ox - x y n‘COS/01_ sin。？、00 )(%)_ sing]0y=、sin&2\\z丿00COS&]sin0]0COS/ ><0ycos% 丿\\z丿，cos&2 - sin0]COS020

一 cos&]Sin02、cos%sin0]Cos&2、sin&2

yCOS&]COS&2 >\\z丿-1-sin%反之，逆变换为:N、‘COS02一 sin°]Cos&2COS0]sin°]Cos&2一 cos&]Sin02、-sin&iCOS&]COS02 丿0y=\\z丿、siri02y\\z丿2.3椭圆方程的拟合方法及中心提取由于塔身在变形前观测点处于同一水平面上，假设7T2就是这个平面，并且观测点P,在水平面772上对应的位置分别为P： ( i = 1 ,2,m),那么Pi ( i = 1,2,…,m)必在一个正八边形的外接圆上 而

Pj( i = 1,2, ,/n)是匕(i二l,2,・・・,m)在平面771上的投影,P：( i = l,2,・・・,m)在一个椭圆上，且椭圆的短轴与长轴长度之比为COS0 .设椭圆的方程为：x2 + cos2^ • y2 + ax + cos2 0 • by + c = 0.通过空间坐标变换求出Pi( i = 1,2,…,m)在0x -x y 〃坐标系中的坐标为兀=(%, ,x2 *• 9xm , y = (71 '』2 '，…，九V，则拟合该椭圆方程可表示为如下矩阵形式：((%I )2 + cos2^ •(y；)2、(-x - cos2 0 • y - 1)b=\\c丿两边同时左乘(-兀-cos2^ • y - IL,得:<(%8)2 + COS% ・(yg)2?44西南民族大学学报（自然科学版）第46卷X cos% • xiyiX cos'e • xiyiX cos% • y-X cos'。•求得拟合系数(一》务(兀? + cos'&・y?)X cos2<9 •儿b=一 丫 cos% • y. (+ cos2^ ・ yn丿\\c丿< _ 丫分 + cos2^ ・ y-丿工 cos2^ ・ xiyi_ i-工省(彳 + COS2^ • y -)b=\\c丿Xcos2^ • xiyiX cos\"e • yi丫 cos% ・ y-cos2^ • y-一工 cos2 3 • Yi (xj + cos2^ • yj)In丿< _ 》球 + cos2^ ・ y- >拟合椭圆中心坐标为0」=（一号，-£），设椭圆中心的z坐标为z'，通过空间坐标逆变换，求得关于z的

在原始坐标系下的中心坐标0 =（x0,y0,z0），再将（勺小凤）代入平面眄方程，解得z和, （%0,y0,z0）即是 该层塔中心位置.3模型的求解3.1实例分析及数学模型的建立(a)变形前(a) before deformation(b) after deformation图1塔身受载前后变形对比示意图Fig. 1 Schematic diagram of the deformation of the tower观察图1可以推测，塔身的截面形状大致为正八边形，其外接圆的直径大约为10 ~ 12m,而观测的参考点 恰选取在正八边形的各个顶点附近•同时，由于受到复杂的外载荷以及恶劣地质条件的影响，塔身已发生明显

的倾斜•如图1所示，这种倾斜应当理解为弯曲，即塔身受到弯矩作用,其外表面有拉压应力的存在，而内应力 必然导致结构的弹性变形「⑷•因此，塔身的截面不再保持水平，并且截面的外接圆也不再是完美的圆形.考虑一个空间椭圆曲线与水平面呈一定夹角8，用这样的椭圆曲线方程来拟合8个空间参考点，以此椭圆

曲线的中心点作为塔层的中心位置点.利用塔的每一层已知的8个观测点P,（i = 1,2,-,8）的坐标确定古塔 各层中心位置的坐标•根据8个观测点的空间位置，合理推测塔的横截面为正八边形，因此其截面外围存在一

个外接圆C .首先用一个空间平面77,来拟合这8个观测点，得到拟合平面的方程及法向量n .考虑到塔身受载 变形后，横截面挤压变形,其截面外围的外接圆C变成椭圆C'，且该椭圆位于77,平面上.然后在77,平面上建立

局部坐标系，通过坐标变换得到p；（i = 1,2,-,8）在局部坐标系下的坐标•然后利用椭圆曲线拟合这8个观测

点，得到椭圆的方程，椭圆的中心即为该塔层的中心，利用坐标逆变换得到在原始坐标系下的中心坐标.第1期范昌胜，等:提取空间拟合椭圆模型的中心坐标算法45Fig. 2 Spatial ellipse fitting model如图2示，以塔身的某一层为例,建立空间椭圆曲线拟合模型，得到该层的中心点坐标，其他层的中心坐标

同理可得.3.2算法的实现以1986年观测的第一层数据为例，计算其层中心位置坐标.观测的P, ~ Ps的坐标如表1所示（已做坐标

平移处理）.表1第一层观测点坐标(1986年观测)Table 1 Coordinates of the first layer of observation points ( observed in 1986)坐标点x/my/mz/m115.45428.0121.792212.05825.5441.818311.3921.4471.783413.78218.1081.769517.94117.4071.772621.25519.8571.77721.93823.9531.794819.527.3561.801首先设拟合平面77!方程为 4% + By + Cz + 1 = 0，原始观测点的坐标为X =\"（,兀2，° ° •，兀8）'（X』2，8）T，Z =（ Z],勺，…,Zg）t ,该拟合表7K成如下形式:71y2MAn-1I- 1Jy8zz丿两边同时左乘（x y

，得:46西南民族大学学报（自然科学版）第46卷<23393028238 VA )(-133.3、- 181.71- 14. 30 丿3028<238求得拟合系数A、B、C：424432532525. 6 >B=(2339B=所以平面叭的方程为：3028238 V1(-133.3]'-5. 1105 x 10\"、3028<238424432532525. 6 J-181.7=14. 30 丿、

1.9397 x IO--0.57935

>-5. 1105 x 10+ 1. 9397 x 10\"3y - 0. 57935z + 1=0.其法向量n =（A,B,C）t =（ -5. 1105 x IO-4,1.9397 x IO-3, - 0.57935）1

平面77）与平面兀之间的夹角为0= 0. 003462 « 0. 1984°.为了使原始坐标系的z轴与法向量；重合，需要依次绕x轴和y轴旋转的角度分别为:厲=0. 003348 « 0. 19183°

02 = 0. 000882 = 0. 0505° '因此从0 - xyz系变换到0} - x'yn系的坐标变换为：(Xs(1=- 0. 00335-0. 00088 Ay0<0. 000881- 0. 00335\\z丿其逆变换为:0. 00335y1 >\\z丿r 10y=\\z丿<一 0. 000880. 003351一 0. 003350. 00088 A0.00335y1丿丿通过空间坐标变换求出= 1 ,2,---,8）在0 - x'yn坐标系中的坐标如表2所示.表2坐标变换后的第一层观测点坐标(1986年观测)Table 2 Coordinates of the first layer of observations after coordinate transformation ( observed in 1986)点123坐标x/my/m28.0058425.53777z/m1.89940915.3586311.970871.914151.8648441.8417741.84609611.3166213.7198121.440914518.1019817.4009719.8509623.9468627.3498217.8811521.186956781.85522221.8562119.406811.8935381.909782设椭圆的方程为异+ COS20 •/+«% + cos20 ■ by+ c = 0，则拟合该方程可表示为如下矩阵形式g :第1期范昌胜，等：提取空间拟合椭圆模型的中心坐标算法47—%!' - cos2 3 ・ y'x-1)fay气珀）2 + cos2^ ・（y；）2、b=l - %8' - cos2 n0 • y，8_ 1丿\\c丿J%8)2 + cos?。・(y；)2丿两边同时左乘（-%—cos20 • y - 1)7 ,得:<2318. 433012. 69132. 70)(-112711A3012. 694242. 17181.63b=-154321I 132. 70181.638 >\\c丿< 一 6560. 6 >求得拟合系数'Q、<2318.43 3012.69 132.70]-i(-11271H33. 17497)b=3012.69 4242. 17 181.63-154321=-45. 40562\\C丿、132.70

181.63 8

>、一 6560. 6><761. 09762 >所拟合椭圆中心坐标为0.（-号，-纟）=（16.58749,22.70281），设椭圆中心的z坐标为z'，对(16.58749,22.70281,z ')进行逆变换:仏、( 10. 003350. 00088 )(16.58749 ATo=010.0033522. 70281vzo丿、-0. 00088-0. 003351 >\\ z J'16.66354 +0. OOO88z'、22. 70281 +0.00335/ •、

-0. 09065 + z' >将(％o,yo,Zo)代入平面771方程-5. 1105 x 10「3 + 1. 9397 x 10_3y - 0. 57935z + 1 = 0 ,解得z' = 1.90745，从而(x0,y0,z0) =(16. 6652,22. 7028,1.7874) •将坐标平移回到观测坐标系，即该层塔的中 心点坐标为(566. 6652,522. 7028,1. 7874)•各次测量的古塔各层中心坐标如表3表3古塔各层中心坐标（1986年观测）Table 3 Center coordinates of the various layers of the ancient tower ( observed in 1986)坐标层 一x/my/mz/m1566.5875522.70281.7873752566.6346522.64477.320253566.6832522.587112.755254566.7198522.540117.078255566.7618522.48921.72056566.3341522.401526.235137566.2622522.35829.83688488西南民族大学学报(自然科学版)566.1815566.0856566.1523522.3139522.268第46卷33.35088936.8548840.1721344.4408848.7118852.834291011522.4464566.0685565.9374522.3588522.29591213565.638522.03831996年、2009年、2011年对应各层中心的坐标,可以类似得出，这里省略.4模型评价1) 在求解各层中心位置点坐标时,先用一个平面对原始观测数据进行空间平面拟合,把原本三维的问题降

低到二维来考虑，使得原本的空间曲面问题简化成了平面问题.2) 考虑到了塔身受载变形后截面形状的变化,采用椭圆模型对观测数据点进行拟合,更符合实际情况,拟

合结果更加准确•同时对原始观测的坐标点进行坐标变换，使得数据点可以用一个标准椭圆方程来拟合，简化

了构造椭圆方程的形式.3) 本算法可以用来解决空间建筑物的中心位置进行标定的问题.不仅解决了如何采用常规测量仪器和测

量方法进行古塔倾斜观测的难题，而且也适用于其他三维空间的建筑变形观测(如圆形构筑物、高层建筑等)•仿真实验和实验表明，本文算法比其它的算法具有速度快和推广价值高等优点，算法总体时间复杂度为0

(nP).在求取椭圆的过程中有效剔出杂质点，准确度得到提高的同时,系统的鲁棒性也得到增强.参考文献[1] 张炜，徐磊.陕西省古塔现状调査及研究[J].文博，2012 (2)：

64 - 70.[9] ZHOU J, SHI Z, FU Q. Three - dimensional scattering center extrac­

tion based on wide aperture data at a single elevation [ J ]. IEEE Trans

Geoscience and Remote Sensing, 2015 , 53(3) : 1638 - 1655.[2] 卢俊龙.湿陷性黄土地区砖石古塔纠偏技术研究[D].西安:西安

[10] KIM K T, KIM H T. Two - dimensional scattering center extraction

based on multiple elastic modules network [J]. IEEE transactions on

建筑科技大学,2005.[3] 胡志晓.古塔倾斜观测和数据分析[J].江苏建筑,2011, 6： 34-

37.antennas and propagation, 2003 , 51(4) : 848 -861.[11 ]NAIR P S, SAUNDERS JR A T. Hough transform based ellipse de­

[4] 黄人堂，卫启云.GPS在应县木塔变形监测网中的应用研究[J].

tection algorithm [ J ]. Pattern Recognition Letters, 1996, 17(7)： 777 -784.北京测绘，2000 ( 2)： 28 -33.[5] 李良福，冯祖仁，贺凯良.一种基于随机Hough变换的椭圆检测算

[12] ZHOU F, CHEN Q, ZHANG G. Composite image processing for cen­

ter extraction of structured light stripe [ J ]. Journal of Optoelectronics.

法研究[J] •模式识别与人工智能,2005,18(4):459 - 461.[6] 闫蓿王斌李媛.基于最小二乘法的椭圆拟合改进算法[J].北京

Laser, 2008, 11 : 024.航空航天大学学报,2008,34(3) ：295 - 29&[7] .安新源，周宗潭，胡德文.椭圆拟合的非线性最小二乘方法[J].

[13] CHEN T C, CHUNG K L. An efficient randomized algorithm for de­

tecting circles[ J]. Computer Vision and Image Understanding, 2001 ,

计算机工程与应用，2009, 45(18)： 188 -190.[8 ] BHALLA R, LING H. Three - dimensional scattering center extraction

using the shooting and bouncing ray techniquef J] . IEEE Transactions83(2) : 172 -191.[14] 尤承业.解析几何[D].北京:北京大学出版社,2004.[15] 周建兴，等.Matlab A门到精通[M].北京：人民邮电出版社，

2009 年.(责任编辑：张阳，付强，李建忠，罗敏;英文编辑：周序林)on Antennas and Propagation, 1996, 44( 11 ) : 1445 - 1453.