台风年鉴十表元数据

模型算法说明文档

算法名称：拟合多边形为椭圆算法

作者：白舒扬

测试人：徐宏

完成日期：2010年6月19日

签收人：

签收日期：2010年月日

1算法简述 本算法将多边形拟合为椭圆。算法利用带约束的最小二乘法，由给定的若干顶点得到椭圆在平面直角坐标系的一般方程，继而求得椭圆各特征参数。算法特点为稳定和快速。

2算法输入

◆ x ：多边形顶点横坐标组成的列数组 ◆ y: 多边形顶点纵坐标组成的列数组

若给待拟合多边形顶点编上序号，则x 和y 的第i 个元素分别是第i 个顶点的横和纵坐 标。输入的点顺序不一定按多边形周边的绕行顺序。一个顶点仅输入一次，不能重复。

3数学原理

◆ 椭圆方程及其变换

在平面直角坐标系上，一般位置的椭圆可表示为关于,x y 的方程：

220ax bxy cy dx ey f +++++=， 132(0,0)I I I <> (1)

的解。其中

1I a c

=+

222

b a

I b c

=

322222

2

b d a b e I

c

d

e f

=

通过旋转和平移，可将一般位置的椭圆方程变换为中心在原点、对称轴重合于坐标轴 标准椭圆方程，从而提取出椭圆的特征参数。

首先假设（1）中的各系数已知。现在把坐标轴逆时针转动θ角度，其中(',')x y 表示点

(,)x y 在坐标轴变动后的新坐标：

'cos 'sin 'sin 'cos x x y y x y θθ

θθ=-⎧⎨

=+⎩

(2) 把(2)代入(1)，设新方程为：

22''''''''''''0a x b x y c y d x e y f +++++= (3)

容易求出：

22'cos cos sin sin a a b c θθθθ=++

()

'sin 2cos 22

c a b b θθ-=

+ 22'sin cos sin cos c a b b θθθθ=-+

'cos sin d d e θθ=+ 'sin cos e d e θθ=-+

'f f =

为了使坐标轴变换后，方程不在出现''x y 一项，即'0b =，则有：

arctan(

)b

a c

θ=- 为了让''x y 项的系数为0，需令

1= arctan(

)2b

a c

θ- (4)

在将(4)代入(3)中各系数确定它们的值。此时再对(3)进行配方得到

22

22'''''['()]'['()]'02'2'4'4'd e d e a x b y f a c a c --+--+--=

令22

'''4'4'

d e F f a c =-- 则有

22

0022

('')('')1x x y y A B

--+= 其中，

A =

B =

0''2'd x a =- 0'

'2'e y c =- （5） （（1）中的限制条件保证了根号下为正）

此时的A 和B 即椭圆的两半轴长，再利用(2)将目前的椭圆中心点坐标00(',')x y 转变为坐标轴变换之前的00(,)x y 。

于是，我们从椭圆一般方程的系数得到了椭圆的两半轴长，中心坐标以及A 半轴关于x

坐标轴的角度（逆时针为正向）。

用现在得到的参数，椭圆在原坐标系下的参数（t ）方程为：

00cos sin cos sin cos sin x x A t y y B t θθθ

θ⎧-⎡⎤⎡⎤⎡⎤⎡⎤

=+⎨⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦

⎣⎦⎩

带约束的最小二乘

首先考虑带约束的方程：

220

ax bxy cy dx ey f +++++=（2241I ac b =-=） (6) 之所以可以把约束20I >换为21I =，是因为若(6)中方程系数都乘以某个非零数，方

程解没有变化，所以无妨预先假定21I =。

而最小二乘拟合，实际上是求函数

2

221

(,,,,,)

()n

i

i i i i i i F a b c d e f ax

bx y cy dx ey f =+++++∑ 的最小值，其中(,)i i x y 为第i

个待拟合顶点坐标，另外还带有41ac b -=的约束条件。为简化，记

22[,,,,,1]i i i i i i i V x x y y x y =,[,,,,,]P a b c d e f =。再令

12V D V ⎡⎤

⎢⎥=⎢⎥

⎢⎥⎣⎦

，于是问题转化为

求2

min P

DP ，1T PCP =

其中，

002000010000200000000000000000000000C ⎡⎤⎢⎥-⎢

⎥⎢⎥

=⎢

⎥⎢⎥⎢⎥

⎢

⎥

⎣⎦

由拉格朗日乘子法，转化为解方程组：

T T T D DP CP λ=

（7）

1T PCP =

（8）

接下来通过矩阵分块和求特征向量的方式能解出上述方程组（详见参考文献1），从而得到所需的系数组,,,,,a b c d e f 。

现在只需要说明得到的此系数组满足约束130I I <。而实际上，方程(6)在实平面上有非单点解当且仅当130I I <（见参考文献2，P158），或者说仅需(6)有非单点解，则解一定为椭圆。我们又有以下结论：

【结论】有上述方法解出的[,,,,,]P a b c d e f =必定使得方程(6)在实平面上有非单点解。

证明：记2

2

(,)T

l V P VP ax bxy cy dx ey f ==+++++

此时，由拉格朗日乘子法，仅考虑(6)最后一行的那个方程有：