古塔变形

古塔变形
摘要
本文主要通过对倾斜、弯曲、扭曲三个角度研究，来分析古塔的变形趋势。

针对问题一，由于观测数据缺失现象，本文先进行数据的预处理，再采用空间圆的拟合方法，建立空间球体和斜平面所组成的方程组来求解。

同时结合运用MATLAB软件，计算出4次测量的古塔各层中心坐标，并以列表形式给出。

针对问题二，本文假设古塔的第一层是水平的，首先过古塔的塔顶的中心点作第一层的垂线，作出射影，构建直角三角形，得到射影和垂线的比值作为衡量古塔倾斜程度的指标。

发现随着年份增加，古塔倾斜程度越大，但1986年—1996年、2009年—2011年的变化幅度不大，1996年—2009年的变化幅度较大。

对于研究弯曲度，将各层的中心连线向X-Z平面及Y-Z平面进行投影，分别比较两者的方差得出向X-Z平面投影的研究价值更大。

通过曲线拟合求得4年曲线对应的方程，利用曲率公式求得曲线上各点的曲率值来分析古塔的弯曲程度，得出从1968年到1996年，各层的曲率有增长趋势，从1996年-2011年各层的曲率逐渐递减。

最后，将相邻两层的中心连线投影到第一层中，得到的夹角大小来确定古塔的扭曲情况，得到各年各层的扭曲率不存在明显的递变规律。

针对问题三，本文从倾斜、弯曲、扭曲三个方面采用灰色预测模型GM(1,1) 进行拟合与预测，并得出各因素相对于各层的时间响应函数，经检验得出扭曲程度最容易预测，其次是倾斜程度，而弯曲程度不易预测。

关键词：古塔变形；空间圆拟合；投影；灰色预测模型GM(1,1)
一．问题重述
塔，作为一种古代高层建筑，不但标志着古代人们征服自然的胜利，同时也记录了当时的历史和工艺。

塔的挺拔高耸的姿态，对佛寺组群甚至城市轮廓面貌都起着一定的美化作用，也成为今天的地标性建筑。

塔采用多种结构，也有变化多端的造型，在中国建筑史上独树一帜，是我国古代优秀的文化遗产。

但是这类建筑物由于长细比大、重心高，且由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响等。

因此对地基要求非常高，因而许多古塔都产生了不同程度的倾斜。

特别是砖石结构的古塔，自重大、强度低，一旦倾斜，往往造成倾斜一侧地基应力集中，又进一步加速了建筑物的倾斜。

虽然塔身的倾斜在美学上营造了一种缺憾美，但在古建筑的保护中应该遵循‘保护第一”的原则。

因此，对于倾斜古塔的纠偏应该慎重，在保证安全的前提下，应尽量保持其原貌。

为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以致定必要的保护措施。

我国目前对于倾斜古塔的保护通常采用扶正纠偏的方法，而对古塔的研究也大多集中于如何扶正以及在非倾斜状态下结构性能的分析上。

某古塔已有上千年历史，是我国重点保护文物。

管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。

根据4次观测数据，进行分析
1. 给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。

2. 分析该塔倾斜、弯曲、扭曲等变形情况。

3. 分析该塔的变形趋势。

二．问题假设
1.假设古塔的塔体是完全刚性的；
2.假设观测点均具有代表性；
3.假设古塔第一层是水平的；
4.假设观测方法正确且测量仪器精度较高；
三．符号说明
()c b a ,,： 空间球体的球心坐标；
i ： 表示塔的层数；
∂∆: 第i 层的中心位置相对第一层水平偏移的距离 H ∆: 第i 层到下一层平面的距离；
1→i α： 第i 层塔与底层塔形成的倾斜度；
),,(i i i z y x ：表示第i 层的中心位置坐标； n ： 相邻两个中心位置点向量；
'z ： 曲线拟合函数一阶求导； ''z ： 曲线拟合函数二阶求导；
1→n θ: 第n 向量与第1向量的角度；
K : 该塔的弯曲度；
四.问题分析
针对问题一，本文首先对提供的观测数据中缺失的数据进行预处理，并对塔尖中心位置点坐标运用求平均值的方法加以确定，然后对每层给出的8组观测数据进行球体拟合和斜平面拟合，分别确定其方程，再利用球心到平面的投影，可得到空间圆的圆心坐标即中心位置点的坐标。

针对问题二，本文要分析该塔倾斜、弯曲、扭曲等变形情况，首先，根据倾斜角建立倾斜程度模型，计算出各年第i 层对底层之间的倾斜角的大小并进行分析。

其次，将各层中心位置点投影到z x -平面，对其z x 、分量数据进行多项式拟合，求得各年中心位置拟合函数方程，建立曲率模型，求出每年各层中心位置点的曲率的大小并进行分析。

然后，将空间各点投影到y x -平面，计算出相邻两点形成的向量，并通过第n 向量和第
1 向量之间的夹角，进行扭曲程度分析。

针对问题三，本文考虑提供数据年数较少，因此采用灰色系统模型GM （1,1）分别对三个影响因素进行拟合、预测，得出相应的响应函数，并对各自时间响应函数进行精度检验，从而分析该塔的变形趋势。

五．模型的建立与求解
5.1 问题一的模型建立与求解 5.1.1 模型准备
1.缺失点的坐标处理
经过观察提供观测数据，可以得到1986年与1996年第13层塔中只有7组数据。

考虑到与12~1层塔所给的数据相同，则通过对所给的点的数据用求平均值的方法确定该点的数据。

得出数据如表1：
表1 确定缺失点的坐标
2.塔尖中心点坐标的处理
经过观察提供观测数据，可以得到各年中的塔顶给出的数据较少。

且给出的点的坐标几乎接近。

再根据古塔的形状，塔顶一般都是尖的，可以近似于一个点，采用求平均值的方法确定塔顶的坐标。

求得数据如表2：
5.1.2 用空间圆的拟合方法确定中心位置的坐标
通过MATLAB 软件分析古塔各层观测点的坐标的空间位置，发现各层的观测点几乎落在空间的一个圆上，则观测点的中心为空间圆的圆心。

要求此圆心，可通过建立空间
球体和空间平面所组成的方程组来求解。

因此，分以下三步进行：（1）球体拟合；（2）斜平面拟合；（3）球心到平面的投影。

下面以求解古塔在1986年时第一层的中心坐标为例，以此类推，同理可得各次测量的古塔各层中心坐标。

1.球体的拟合
球的一般式方程：
0222=+---++D Cz By Ax z y x （1） 现在同层有8个观测点，要拟合这个球体，可写成矩阵的形式为：
⎥⎥⎥⎦⎤
⎢⎢⎢⎣⎡++++=⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--2828282121218
8
8
111
11z y x z y x D C B A z y x z y x
（2） 经过矩阵变换可得:
⎥⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢
⎢⎢⎢
⎣⎡++-++++++⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢
⎢⎢
⎣
⎡------=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∑
∑∑∑∑
∑∑∑∑∑∑∑∑∑∑∑∑∑∑-⨯)()()()(8
2222222222221
442
22
i i i i i i i i i i i i i i i i i i i
i i
i i
i i
i i i i
i i i
i i i i z y x z y x z z y x y z y x x z y x z
z z y z x y z y y y x x z
x y x x D C B A （3） 将8个观测点的数据代入（3）可解得：
32.1133=A ，44.1045=B ，5656.2-=C ，594305=D 则可得到球心坐标：),,(c b a =（566.66，522.72，-1.28179） 2．斜平面拟合
对所测量的斜平面点进行空间平面拟合，以得到空间平面方程。

空间平面的方程通常表示为：
01'''=+++z C y B x A （4） 现在同层有8个观测点，要拟合这个平面，可写成矩阵的形式为：
⎥⎥⎥⎦⎤⎢
⎢⎢⎣⎡---=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡111'''88
8
111
C B A z y x z y x
（5）
两边左乘T
z y x z y x ⎥⎥⎥⎦
⎤
⎢⎢⎢⎣⎡88
8
11
1
得： ⎥⎥⎥
⎦
⎤⎢⎢⎢⎣⎡---=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤
⎢⎢⎢
⎣
⎡∑∑∑∑∑∑∑∑∑∑∑∑i
i
i
i
i
i i i i i i i
i i
i i
i i z y x C B A z z y z x z y y y x z x y x x '
'
'22
2 （6）
可得出：
⎥
⎥⎥⎦⎤
⎢⎢⎢⎣
⎡---⎥⎥⎥⎦⎤
⎢⎢⎢
⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡∑∑∑∑∑∑∑∑∑∑∑∑-⨯i i i i i
i i i i i i
i
i i i i
i i z y x z z
y z x z y y y x z x y x x C B A 1
33222''' (7） 将8个观测点的数据代入可解得系数A '=-0.00101，B '=-0.000337，C '=-0.140567。

即可得到相应的空间平面的方程：
1140567.0000337.000101.0=++z y x （8）
3 球心到平面的投影计算
球心到平面的投影计算可表示为：
'''C
c
z B b y A a x -=-=- (9)
联立（1）（4）（9）式得：
⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=+++-=-=-=+---++0
1''''''0222z C y B x A C c
z B b y A a x D Cz By Ax z y x （10） 由以上方程组联立即可得出空间圆的圆心坐标（566.682,522.727,1.789），即为
第一层中心坐标。

5.1.3 各层中心位置坐标的求解
根据上述算法，利用MATLAB 软件（程序见附件一）求解，得到4次测量的古塔各层中心坐标，见表3。

表3 4次测量的古塔各层中心的坐标
年份 2009年 2011年 塔层 x/m
y/m
z/m x/m y/m z/m 1 566.7441 522.6969 1.7635 566.745 522.698 1.766 2 566.7819 522.6742 7.3116 566.783 522.672 7.290 3 566.8145 522.6488 12.7287 566.815 522.650 12.731 4 566.8364 522.6208 17.0749 566.840 522.621 17.050 5 566.8696 522.6009 21.7137 566.871 522.603 21.711 6 566.9554
522.551
26.2045 566.957 522.550 26.203 7 566.9887 522.5278 29.8161 566.989 522.526 29.801 8 567.0397 522.4932 33.3313 567.042 522.496 33.315 9 567.0923 522.4637 36.8253 567.092 522.460 36.822 10 567.1423 522.4024 40.1507 567.144 522.400 40.142 11 567.1887 522.3667 44.4034 567.186 522.365 44.414 12 567.2289 522.3297 48.7064 567.230 522.327 48.692 13 567.2829 522.2861 52.8102 567.283 522.282 52.813 塔顶
567.336
522.2148
55.091
567.338
522.214
55.087
分析表3数据发现：x 轴上的坐标值随着塔层的增加在逐渐递增,而y 轴上的坐标值随着塔层的增加在逐渐递减；各层中x 轴上的坐标值随着年份的增加在逐渐递增，而各层中y 轴上的坐标值随着年份的增加在逐渐递减。

这些坐标值变化，说明古塔在外界及自身条件下不断地发生不同程度的倾斜、弯曲、扭曲等变形现象。

5.2 问题二的模型建立与求解 5.2.1建立塔倾斜程度模型
由于先前假设古塔的第一层是水平的，要分析塔的倾斜程度，只需把古塔的第一层中心作为塔底的中心，将其与其他塔顶的中心连线，并过塔顶的中心做第一层的垂线，得到射影，即水平位移量，形成直角三角形，利用连线与垂线的夹角α来分析古塔的倾斜程度，如图1所示。

图1 倾斜程度示意图
由两点坐标距离公式得：
()()212212y y x x -+-=
∂∆ （11）
12z z H -=∆ （12）
由三角函数关系可得倾斜角度α为：
=αH
∆∂
∆arctan
（13） 同时，通过MATLAB 软件（程序见附件二）求解每年第i 层与底层之间的倾斜角大小，见表4。

表4 各年每层与低层之间的倾斜度
从表4中分析数据，得出：各年各层的倾斜度变化均不等，且不存在递变规律，这是由于受自重，风力，气温等自然因素影响程度不同所造成的。

然而，2009年与2011年分别在第六层、第十层的倾斜度相对较大，说明这两年的环境对古塔的影响更为恶劣，随着年份的递增，其中塔尖相对于下一层塔的倾斜度在递减，说明塔产生了倾斜现象。

5.2.2建立弯曲程度模型 由表4看到，每一层相对于底层的倾斜角度均不同，说明该古塔不仅出现倾斜现象，
同时出现弯曲现象。

1.对数据进行处理
从表3中的数据可以看出z 轴的变化比较明显，再通过求方差，判断x 、y 轴的偏离程度。

求得方差见表5。

表5 各年x 、y 轴方差表
1986年 1996年 2009年 2011年 x 0.409 0.422 0.498 0.496 y
0.209
0.213
0.303
0.308
方差值越大，偏离程度就越高。

方差值越小，偏离程度就越低。

从表5中数据可以得出x 轴数据偏离程度较大，对塔的弯曲程度影响相对较大。

从而，将塔投影到z x -平面上进行弯曲程度分析。

2.对z x ,分量数据进行多项式拟合
通过MATLAB 软件（程序附件三）对z x ,轴数据进行多项式拟合求得各年曲线拟合函数表达式分别为：
1986年： 1.823794547
7915.289644598.2521-+-=x x z （14） 1996年： 782.77617672911.272869807.2322-+-=x x z （15） 2009年： 3103.16443629398.579139916.5023-+-=x x z （16） 2011年： 2613.16539273295.582502884.5124-+-=x x z （17）
3.曲率的求解
设曲线C 是光滑的, 曲线上点M 对应于弧s , 在点M 处切线的倾角为α, 曲线上另外一点N 对应于弧s s ∆+，在点N 处切线的倾角为αα∆+。

图2 曲率分析图
通过计算，得到曲率方程为：
2/320)
1(lim
z z s K s '+'
'=∆∆=→∆α （18） x
将各层中心点的坐标值代入到各年份的z''、z'的表达式中，再代入到（18）式，求出相应的曲率（见表6）。

观察其弯曲程度。

分析表6的数据，得出：各年各层的曲率都随着层数的递增在逐渐增加；说明古塔有弯曲现象。

从1968年到1996年，各层的曲率有增长趋势，说明在此期间受到恶劣天气等影响。

而从1996年-2011年各层的曲率在逐渐递减。

说明古塔在此期间受到恶劣天气的影响较小。

5.2.3建立扭曲程度模型
1. 判断塔的扭曲程度的依据
x-平面上，根据投影点存首先将三维空间转化成二维空间。

将空间的点投影到y
在形式确定塔是否扭曲。

x-平面上的点之间的连线不存在一个夹角的时候，说明该塔没
(1).如果投影到y
有发生扭曲，但发生了弯曲。

x-平面上的点之间的连线存在一个夹角的时候，说明该塔发生
(2).如果投影到y
扭曲。

x-平面上的点之间的连线形成一个点的时候，说明该塔既没有
(3).如果投影到y
发生扭曲，又没有弯曲。

2. 确定向量，并求出向量之间的夹角
(1).向量的确定
x-平面上记作将在模型一中求出各层之间的中心位置点的坐标，投影到y
()0,,i i y x ，则可以得到()11,----=i i i i y y x x n ，()2121,1------=-i i i i y y x x n
(2).求解向量之间的夹角
2
212
212
12
1211211)
()()()())(())((1cos -------------+-•-+---+--=
-=
i i i i i i i i i i i i i i i i y y x x y y x x y y y y x x x x n n θ （19）
3.求解层与层形成向量的夹角
将求得n 、1-n 向量坐标代入（19）式 ，求出相应的夹角（见表7）。

观察其扭曲程度。

分析表7中的数据，得出：各年各层的扭曲率不存在递变规律，而存在很大的极差，第十二层中扭曲率变化量在各年中的变化相对较大，而第十三层中扭曲率变化量在各年中的变化相对较小，说明各年自然环境恶劣程度对各层扭曲率大小影响不同。

5.3 预测古塔的变形趋势
对该塔的变形趋势进行预测，从三个方面考虑（倾斜、弯曲、扭曲），由于只有四年，数据较少，且变化趋势也不明显，难以单纯的从拟合的角度进行发展趋势的预测，
而灰色模型系统可以解决这类问题，因此，利用灰色系统模型GM(1,1)的方法来研究相应问题。

5.3.1 建立灰色系统模型GM(1,1)
以第一层的倾斜度为例，其他各层以此类推。

1.给定观测数据列：
}478.0,459.0,706.0,696.0{)0(=x
2.经一次累加序列得：
}339.2,861.1,402.1,696.0{)1(=x
3.建立矩阵：y B ,
[][][]⎥⎥⎥⎥⎥
⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡+-+-+-=1861.1339.2211402.1861.1211696.0402.121B []T y 478.0,459.0,706.0=
4.由()y B B B U
T T 1
ˆ-=，求估值a ˆ和u ˆ. ()
⎥⎦
⎤
⎢⎣⎡-==⎥⎦⎤⎢⎣⎡=-7978.00092.0ˆˆˆ1
y B B B u
a U T T 把a
ˆ和u ˆ代入时间响应方程，即时间响应方程为： 1)
--0.2252(t 3.3294e 4.0254+=y （20）
5.3.2 从倾斜程度分析变形趋势
1．通过灰色系统模型对第i 层与第1层之间的倾斜度，通过MATLAB 软件（程序见附件四）对GM （1，1）进行拟合得到表8。

表8各年第
层与第1层之间的倾斜度拟合
分析表8数据，发现：随着时间的变化8
~
2层相对第1层的倾斜度逐年减少，则说明塔相对稳定，而从10
~
9层，倾斜度逐年增大，从而可以得到倾斜度与塔的高度有一定的关系。

2．精度检验与预测
通过查阅资料得出最小误差概率与差比值精度等级表（见表9）。

由最小误差概率与差比值的值检验GM（1,1）模型的预测精度，已提供决策依据。

精度等级越小越好，精度一级，表示预测具有较高的精度，四级为不通过。

通过MATLAB软件对上述结果进行检验
根据表9判断标准，对表10中的数据进行分析得到：关于古塔各层倾斜程度的大部分预测模型等级合格及以上水平，而6层及11层的预测效果不好。

5.3.3 从弯曲程度分析变形趋势
对于弯曲程度的分析，采用同样的灰色模型GM（1,1）进行拟合、预测，得到相应的响应函数，并将各响应函数进行精度检验与预测，得到表11。

表11 弯曲度精度检验与预测
根据表9判断标准，对表11中的数据进行分析得到：关于古塔各层弯曲程度的预测模型的拟合、预测效果不太理想。

5.3.4 从扭曲角度分析变形趋势
对于扭曲程度的分析，采用同样的灰色模型GM（1,1）进行拟合、预测，得到相应的响应函数，并将各响应函数进行精度检验与预测，得到表12。

表12 扭曲率精度检验与预测
根据表9判断标准，对表12中的数据进行分析得到：关于古塔各层扭曲程度的预测模型大部分处于一级精度，表明预测、拟合效果较好。

综上分析，引起该古塔变形的三个因素中，扭曲程度最容易被预测，其次是倾斜程
度，弯曲程度不易被预测，而扭曲程度、倾斜程度，弯曲程度之间又是相互影响的。

有关部门应对古塔进行周期性地观测，做到及时监测，有效监测，以便掌握古塔随时间的变化产生的各引起变形的因素，有利于维修工作的展开。

六．模型的评价
6.1 模型的优缺点
1.本模型优点：
（1）利用空间圆的拟合方法确定中心位置的坐标，相对较精确。

（2）利用灰色模型GM(1,1)对影响古塔变形的三个因素进行拟合、预测，得到了相应的响应函数，并对响应函数做精确检验与预测，得出该塔的变形趋势。

2.本模型缺点：
（1）在假设中假设第一层是水平面所得到的数据存在一定的偏差。

（2）在研究扭曲程度的过程中未考虑古塔各中心的竖坐标分量，可能存在一定的偏差。

6.2 模型的推广
本模型可用于诸如上海明珠电视塔的变形情况研究。

七．参考文献
[1] 周永正詹棠森方成鸿邱望仁，数学建模，上海市四平路1239号：同济大学出版社，2010年8月第1版.
[2] 张晶黄琴兰红军孙小东杨勇，工程测量中空间圆的拟合方法与研究，万方数据，2013/9/13.
[3] 司守奎孙玺菁，数学建模算法与应用，北京：国防工业出版社，2011.8
[4] 张德丰雷晓平周燕，MATLAB基础与工程应用，北京：清华大学出版社，2011.12
[5] 陈刚于丹吴迪，MATLAB基础与实例进阶，北京：清华大学出版社，2012.1
[6] 唐丽芳贾冬青孟庆鹏，用MATLAB实现灰色预测GM（1，1）模型，沧州师范专科学校学报，第24卷第2期，2008.6
八．附件
一 Matlab求解中心点坐标程序（以计算1986年第一层为例）：
1.斜平面方程：
x=[565.454 562.058 561.39 563.782 567.941 571.255 571.938 569.5 ];
y=[528.012 525.544 521.447 518.108 517.407 519.857 523.953 527.356];
z=[1.792 1.818 1.783 1.769 1.772 1.77 1.794 1.801];
a=-sum(x);
b=-sum(y);
c=-sum(z);
d=sum(x.*x);
e=sum(x.*y);
f=sum(x.*z);
g=sum(y.*y);
h=sum(y.*z);
i=sum(z.*z);
j=[a;b;c];
k=[d e f;e g h;f h i];
m=inv(k);
w=m*j;
p=vpa(w,6);
P
2.球方程：
x=[565.454 562.058 561.39 563.782 567.941 571.255 571.938 569.5 ];
y=[528.012 525.544 521.447 518.108 517.407 519.857 523.953 527.356];
z=[1.792 1.818 1.783 1.769 1.772 1.77 1.794 1.801];
a=-sum(x);
b=-sum(y);
c=-sum(z);
d=sum(x.*x);
e=sum(x.*y);
f=sum(x.*z);
g=sum(y.*y);
h=sum(y.*z);
i=sum(z.*z);
j=sum(x.*x.*x+x.*y.*y+x.*z.*z);
n=sum(y.*x.*x+y.*y.*y+y.*z.*z);
l=sum(z.*x.*x+z.*y.*y+z.*z.*z);
o=-sum(x.*x+y.*y+z.*z);
p=[j;n;l;o];
k=[d e f a;e g h b;f h i c;a b c 8];
m=inv(k);
w=m*p;
v=vpa(w,6);
V
3.中心坐标：
[x,y,z]=solve('(x-A)/a-(y-B)/b=0','(x-A)/a-(z-C)/c=0','a*x+b*y+c*z+1=0', 'x,y,z');
A=1133.32/2;
B=1045.44/2;
C=-2.56558/2;
a=-0.00101043;
b=-0.000337003;
c=-0.1405672;
x=eval(x)
y=eval(y)
z=eval(z)
运行结果：
1、斜平面方程系数
p =
-0.001010
-0.000337
-0.140567
2、球方程系数
v =
1133.32
1045.44
-2.56558
594305.0
（3）中心坐标
x =
566.682
y =
522.727
z =
1.789
二倾斜角求解程序（以1986年第一层为例）：
x1=[566.682];x2=[566.722];y1=[522.727];y2=[522.673];z1=[1.789];z2=[7.318];
h=z2-z1;
v=sqrt((x2-x1).^2+(y2-y1).^2);
p=atan(v./h);
p*180/pi
运行结果：
ans =
0.6964
三z
x、轴坐标曲线拟合程序(以1986年为例）：
d=[566.682 566.722 566.777 566.819 566.867 566.924 566.959 566.990 567.024 567.050 567.106 567.155 567.200 567.247];
z=[1.789 7.318 12.757 17.084 21.721 26.245 29.837 33.341 36.860
40.169 44.461 48.696 52.838 55.123];
p=polyfit(d,z,2);
px=poly2str(p,'x');
pv=polyval(p,d);
p,px
plot(d,z,d,pv)
运行结果：
p =
1.0e+006 *
-0.0000 0.0290 -8.2379
px =
-25.4598 x^2 + 28964.7915 x - 8237945.4701
曲线拟合图
四倾斜、弯曲、扭曲变形趋势拟合程序：
function [y,p,e]=huise_1_1(X,k)
if nargout>3;
error('Too many output argument.');
end
if nargin==1,k=1;x_orig=X;
elseif nargin==0|nargin>2
error('Wrong number of input arguments.');
end
x_orig=X;
predict=k;
x=cumsum(x_orig);
n=length(x_orig);
for i=1:(n-1);
B(i)=-(x(i)+x(i+1))/2;
end
B=[B' ones(n-1,1)];
for i=1:(n-1);
y(i)=x_orig(i+1);
end
Y=y';
au=(inv(B'*B))*(B'*Y); coef1=au(2)/au(1);
coef2=x_orig(1)-coef1;
coef3=0-au(1);
costr1=num2str(coef1);
costr2=num2str(abs(coef2));
costr3=num2str(coef3);
eq=strcat(costr1,'+',costr2,'e^',costr3,'*(t-1))'); for t=1:(n+predict)
mcv(t)=coef1+coef2*exp(coef3*(t-1));
end
x_mcv0=diff(mcv);
x_mcve=[x_orig(1) x_mcv0]
x_mcv=diff(mcv(1:end-predict));
x_orig_n=x_orig(2:end);
x_c_error=x_orig_n-x_mcv;
x_error=mean(abs(x_c_error./x_orig_n));。