2013全国大学生数学建模山西赛区山西财经大学华商第43队论文----C题-古塔变形

承诺书

我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： C

我们的参赛报名号为（如果赛区设置报名号的话）：

所属学校（请填写完整的全名）：山西财经大学华商学院

参赛队员(打印并签名) ：1. 胡珺怡

2. 赵萌

3. 陈茜

指导教师或指导教师组负责人(打印并签名)：张培强

日期： 2013 年 9 月 16 日赛区评阅编号（由赛区组委会评阅前进行编号）：

2013高教社杯全国大学生数学建模竞赛

编号专用页

赛区评阅编号（由赛区组委会评阅前进行编号）：

全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：

古塔的变形分析

摘要

本文针对古塔的变形这，由于长时间承受自重、气温、风力等各种作用，偶然还要受地震，飓风的影响，古塔会产生各种变形，为保护古塔，文物部门需要对古塔进行观测，了解古塔的变形情况。

因为当3个观测点位于一个圆周上时，可以用图解或解析的方法确定圆心的位置；当观测点数大于3时，可以用最小二乘进行曲线拟合，计算出圆心的坐标。本文用最小二乘法进行曲线的拟合，得到古塔各层中心点坐标。使用线性回归，曲线拟合法，梯度下降法得到各层中心所在的直线方程。对直线进行投影，利用曲率对古塔的变形、弯曲、扭曲进行了研究。

通过对问题一所给数据的处理分析，用matlab对其中已知数据的作图发现所给的数据呈圆形，从而建立圆形塔模型，画出相应图形。发现题设所给的数据所画的图形拟合程度极高。

针对问题二，我们分别建立斜率倾斜模型和曲率弯曲模型，并用matlab进行了有效的验证。

针对问题三，我们发现圆形塔倾斜程度大扭曲程度高，需要有关部门进行维修保养

关键词最小二乘法机器学习法线性回归法圆形塔模型斜率倾斜模型曲率弯曲模型

一、问题的提出

文物部门需要对古塔定期勘测，制定保护方案。由于长时间承受自重、气温、风等外力作用，偶然还要受自然灾害的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。为有效保护古塔，我们在这里以古塔的四次监测数据为例进行探讨。就此现象，建立数学模型，回答以下问题。

1. 找到确定古塔各层中心位置的通用方法，并给出古塔各层中心坐标。

2. 分析古塔的倾斜、弯曲、扭曲情况。

3. 根据古塔的变形趋势来判定古塔的变形程度

二、问题的分析

关于古塔的现状，我们在数据处理研究时以古塔各层中心点和古中心坐标应计量确定为出发点，建立相关模型，在此基础上建立模型解决古塔的倾角、曲率、扭距。

对于问题一，我们观察已知数据进而建立了圆形塔模型，顺利找到塔的中心点。根据观测数据可以看出每一层有八个点，可以把它看成一个八边形，八个点在八边形的边上。我们利用八边形的外接圆寻找其中心点。由于八边形外接圆的圆心和八边形的中心重合，故我们提出的圆形塔模型是合理的

对于问题二，由于长时间承受风力、温度、自然灾害的影响，古塔面临塔结构的变形、扭曲、倾斜等问题。而这些问题可由问题一中的圆形塔模型以及拟合出来的曲线图观测出来。然后利用斜率倾斜模型和曲率弯曲模型可分析得出问题的结果。

对于问题三，通过分析塔中心的变化，对直线方程进行投影，观察塔的倾斜、弯曲、扭曲的情况图表可得出古塔变形趋势，从而得出古塔有受损情况。

三、基本假设

①观测到的数据都具有真实性。

②不考虑风力、湿度、温度对古塔变形的影响

③按Z轴的平均值计算

④题目所给数据真实可靠

四、符号说明及名词定义

五、模型的分析、建立及求解

5.1 问题1模型的提出及分析

使用Matlab对数据进行可视化分析，根据观测数据的二维平面俯瞰分布如图1，提出如下假设：古塔的水平横截面为圆形

图1 观测数据的二维平面俯瞰分布图（参见附件）

5.1.1 问题1模型建立与求解

求古塔各层中心位置的问题可转化为求圆心问题。当3个观测点位于一个圆周上时，可以用图解或解析的方法确定圆心的位置；当观测点数大于3时，可以用最小二乘进行曲线拟合，计算出圆心的坐标。

下面讨论用最小二乘法进行曲线拟合的问题。

（1）模型建立：

从解析几何中可知，圆心在{QUOTE |、半径为的圆方程是：

将方程（1）中的括号展开并移项，得

如果令

则方程（2）变成

只要求出参数就可以求得圆心半径的参数

圆心的Z坐标值C可通过求观测坐标Z分量的均值得到：

样本集中点到圆心的距离为:

点到圆边缘的距离的平方和半径的平方的差为：

令为的平方和：

能够使得的值最小的即为所要求的解。

图2为模型示意图：

图2 模型示意图

（2）模型求解：

平方差大于0，因此函数存在大于或等于0的极小值。对求偏导，令偏导等于0，可求得极值点，比较所有极值点的函数值即可得到最小值。

偏导方程组：

求解方程组：得：

得：

令

可得方程组：

可解得：

进而求得的估计拟合值：

编写matlab函数function[R,A,B,C]=fitted_model(x,y,z,N)计算得到各次测量的古塔各层中心坐标如表1：

表1 古塔各层中心坐标

模型可视化结果如图3：各年份拟合曲线俯瞰图，由图可以看出，拟合曲线完美的拟合了图1中的观测数据，结果准确（参见附件）。

图3 古塔拟合曲线俯瞰图（参见附件）

5.2问题2 模型建立与求解

（1）古塔倾斜分析

线性回归是一种有效的机器学习算法，根据问题1中求得的古塔各年中心的坐标，通过线性回归的方法对中心坐标进行拟合，求得三维空间的一条直线，根据垂直投影的方法求得直线的斜率，从而求出古塔的倾斜度。

（2）数学模型：

假设目标拟合直线方程为：

其中

定义损耗函数：

则满足的即为所求的值。

（3）模型求解：

运用梯度下降法求的值。

梯度下降法：

Repeat{

：＝-

: ＝-

: ＝-

}

重复迭代，直到找到使最小的。

使用matlab仿真（程序见附表2）求得的值如下表：

进而求得拟合直线方程如下表：

（4）结果分析：

取所求直线和面的交点，直线上任意一点，以及点

在面上的投影点

则可求得古塔倾斜角

5.2.1古塔弯曲、扭曲度建模求解

思路

通过将问题1中求得的各层中心坐标投影到和平面内，然后利用曲线拟合的方法，对投影点进行曲线拟合，最后利用曲率的相关数学方法对拟合曲线进行曲率计算。进而求得古塔的弯曲、扭曲程度。

（1）数学模型

对问题1中求得的数据进行从新整理，例如1986年的数据如下表：

以投影到坐标平面为例，假设目标拟合曲线方程为：

定义损耗函数：

则满足的即为所求的值。

（2）模型求解： 运用梯度下降法求的值: 梯度下降法： Repeat{

}

重复迭代，直到找到使最小的。

利用曲率知识定量分析古塔弯曲、扭曲程度： 曲率知识介绍 如图1所示，设曲线上一段弧长为，在点作切线，

在点作切线 ，记切线转过的角度。对于同样的弧长，切线转过的角度。显然，比弯曲程度大。因

此弧长相等时，转角越大，曲率越大。

图1

另一方面，如图2所示，若转角相等时，弧长越短，弯曲程度越大。

M

M

T N a 1

ds M 1N 1T 1

a 2

ds

a

M1

M

N1

N

图2

如上图所示，如果曲线的长度为，两点切线正向所夹的角为，将转角与弧长的比值定义为曲率。

下面给出曲率的计算公式

设曲线的方程为，且具有二阶导数，则曲线的曲率为：

5.3问题3 形变趋势分析

由问题2的建模分析、理论验证以及matlab求解和问题3的仿真作图可得塔的上半部分呈下降趋势，而塔的底部呈上升趋势，由此可知塔身呈一定程度的倾斜。通过塔中心坐标的研究，发现塔身呈一定程度的弯曲和扭曲，需要引起有关部门的重视，加强古塔文物的保护和监测，为古塔提供有力的保障。

图4

八、模型评价

模型评价

优点：模型方法比较具有通用合理性，结果准确，计算精确，科学合理

缺点：本案例所建立的模型对于假设之外的情况存在局限性，有待进一步深化研究。

九、模型推广

该模型可以推广到文物部门对于测量因时间长久及外力作用下发生变形的多种类建筑物。在此模型下,工作人员从时间和空间两个维度下观察构造物的变形，从而简便的得出建筑物的变形趋势。在此基础上制定最优方案对建筑物进行最好的保护。本案例中的这些处理方法可以广泛的应用于古建筑的保护以及已有建筑物的变形监控中。加之，模型本身具有的优势和准确性，不仅有效的利用了测绘公司所测数据得出切实可用的维护方案，而且节约了文物部门的人力、物力、财力，提高了工作效率和工作绩效。而且因此模型有很强的推广性和实用性。

参考文献：

[1]黄强,古塔变形监测的探讨，上海，测绘与空间地理信息，2013.6

[2] 刘珂，周富强，张广军，半径约束最小乘圆拟合方法及误差分析，光电子激光，2006.5

[3] 梁海奎，古塔变形测量方法探讨，城市勘测，2011.6

[4]何凭宗，圆形立柱倾斜度测定的解析方法，矿山测量，2001.12

附件

附件1：

function [R,A,B,C]=fitted_model(x,y,z,N)

%fitted_model.m 最小二乘圆拟合

% x 测量点x坐标值

% y 测量点y坐标值

% N 测量数据集中数据点个数

% R 拟合得到的圆半径

% A 拟合得到的圆心横坐标

% B 拟合得到的圆心纵坐标

x1 = 0;

x2 = 0;

x3 = 0;

y1 = 0;

y2 = 0;

y3 = 0;

x1y1 = 0;

x1y2 = 0;

x2y1 = 0;

for i = 1 : N

x1 = x1 + x(i);

x2 = x2 + x(i)*x(i);

x3 = x3 + x(i)*x(i)*x(i);

y1 = y1 + y(i);

y2 = y2 + y(i)*y(i);

y3 = y3 + y(i)*y(i)*y(i);

x1y1 = x1y1 + x(i)*y(i);

x1y2 = x1y2 + x(i)*y(i)*y(i);

x2y1 = x2y1 + x(i)*x(i)*y(i);

end

C = N * x2 - x1 * x1;

D = N * x1y1 - x1 * y1;

E = N * x3 + N * x1y2 - (x2 + y2) * x1;

G = N * y2 - y1 * y1;

H = N * x2y1 + N * y3 - (x2 + y2) * y1;

a = (H * D - E * G)/(C * G - D * D);

b = (H * C - E * D)/(D * D - G * C);

c = -(a * x1 + b * y1 + x2 + y2)/N;

A = a/(-2); % x 坐标

B = b/(-2); % y 坐标

C = mean(z); % z 坐标

R = sqrt(a * a + b * b - 4 * c)/2;

附件2（线性回归matlab代码）：

function liner_regression

clear ; close all; clc

fprintf('Loading data ...\n');

%% Load Data

data = load('11.txt');

X = data(:, 1:2);

y = data(:, 3);

m = length(y);

% Print out some data points

fprintf('First 10 examples from the dataset: \n'); fprintf(' x = [%.0f %.0f], y = %.0f \n', [X(1:10,:)

y(1:10,:)]');

fprintf('Program paused. Press enter to continue.\n'); pause;

% Scale features and set them to zero mean

fprintf('Normalizing Features ...\n');

[X mu sigma] = featureNormalize(X);

% Add intercept term to X

X = [ones(m, 1) X];

fprintf('Running gradient descent ...\n');

% Choose some alpha value

alpha = 0.01;

num_iters = 400;

% Init Theta and Run Gradient Descent

theta = zeros(3, 1);

[theta] = gradientDescentMulti(X, y, theta, alpha, num_iters); % Display gradient descent's result

fprintf('Theta computed from gradient descent: \n');

fprintf(' %f \n', theta);

fprintf('\n');