关于古塔变形的数学模型

PAGEPAGE1关于古塔变形的数学模型摘要本文主要研究古塔在自重、气温、风力等因素的影响下产生变形的问题。采用中垂线求解外切圆圆心的模型以及多次平均除误差的方法，找到了确定古塔中心的通用方法，并用多元线性回归模型及插值拟合等方法对倾斜、弯曲、扭曲等变形情况进行分析，从而通过残差拟合得出预测数据对古塔变形趋势进行描述。针对问题一：论文采用古塔八个角点中任意三个角点构成的两两连线，取其中垂线的交点得到外接圆圆心，已知正八边形的中心与外接圆圆心一致，但古塔八角点构成的八边形存在轻微不规则，所以我们采用多次取点求外接圆圆心，并用其平均值消除误差，最后对不同取点方式进行了精度分析（答案详见表一）。针对问题二：首先是古塔倾斜分析，根据测量学本文取塔尖和塔底的中心连线作为倾斜角计算的倾斜方程，算出塔顶在水平面投影与塔底中心的间距，引入实测高程数据，可以得到古塔四次测量的倾斜角（），对其倾斜情况经行描述；然后是弯曲情况分析，根据问题一中古塔各层中点坐标，本文对其进行多元回归分析及多项式拟合，得出函数曲线，并将其和倾斜方程进行比较得到最大差值即挠度（材料力学中对弯曲的描述量）；最后是扭曲分析，本文分垂直和水平两个方向进行讨论，垂直方向上涉及高程Z，即对各层中心点多元线性回归得到的拟合值与实测值进行残差分析，得到扭曲描述量（）。水平方向，本文参考材料力学中扭转角的计算，对古塔各层间的轴向扭转进行分析，得到扭转角对古塔扭曲情况进行描述。针对问题三：在分析了四次观测值中倾斜、弯曲，扭曲的情况下，本文采用加权平均的方法各产生影响数据进行处理后，进行残差拟合，得到下一次观测的模拟数据，对古塔的变形进行变形趋势描述关键词多边形中心确定多元回归分析多项式拟合残差分析一、问题重述由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。某古塔已有上千年历史，是我国重点保护文物。管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。请你们根据附件1提供的4次观测数据，讨论以下问题：1.给出确定古塔各层中心位置的通用方法，列表给出各次测量的古塔各层中心坐标。2.分析该塔倾斜、弯曲、扭曲等变形情况。3.分析该塔的变形趋势。二、模型假设根据测量学对倾斜的描述，倾斜角的计算为塔顶与塔底中心连线的倾角，假设各层中心的偏移不影响塔尖和底层的倾斜计算。2、假设古塔的各层面是一个稳定整体，平面转动是同时进行的。3、假设附件给出的数据是八角古塔各角点的坐标。三、符号说明1、：圆上任意三点构成的弦的斜率2、：塔顶在底面的投影与塔底中心的距离3、：塔顶的高程、塔层间距离4、：位移变化量、残差值5、：单位扭转角6、：多元线性回归系数7、：i次项变量8、：参与计算的古塔层数9、：扭曲描述量四、问题分析对于问题一，求解该层的中心，已知数据是八边形古塔的各角点，根据正多边形中心的定义，各对角线的交点以及外接圆的圆心是正多边形的中心，但是对于不规则的多边形，其对角线的交点并不重合，同时外切圆的圆心也会偏离。于是在不同选点情况下，圆心位置可能存在的偏差，我们采取多次取点计算平均值的方式来消除。由于八边形存在的种选点情况，于是我们对不同角点选取对圆心精度变化的情况进行了验证。对于问题二，倾斜分析，在第一问得出的古塔各层中心的基础上，首先按照测量学对倾斜观测的方法，根据塔顶三维坐标与底层中心，可以得出各层高差以及在平面投影上的中心差，由此可得倾斜角：；弯曲分析，在倾斜分析中我们已经得出塔尖和底层中心的连线，可以用各层中心点与该连线的偏离程度来描述古塔的弯曲程度，根据问题一得出的古塔各层中点坐标，通过多项式插值和多元回归分析可以得到中心点连线的回归方程，将此方程分别投影到，平面后，和塔的倾斜方程进行比较得到最大差值即挠度（材料力学中对弯曲的描述量）和差值总量可以得到古塔的弯曲描述；扭曲分析，按照结构力学中扭曲是空间概念的说明，我们分垂直和水平两个方向对古塔的扭曲状况进行讨论，垂直方向上涉及高程Z，即对各层中心点多元线性回归得到的拟合值与实际值进行残差分析，得到扭曲描述量（）。水平方向，参考材料力学中扭转角的计算，对古塔各层间的轴向扭转进行分析，计算得到平面内位移量后通过计算古塔各层高差得到扭转角，对古塔平面扭曲情况进行描述。对于问题三，古塔的变形趋势，首先对前四次测得数据中对倾斜，弯曲，扭曲有影响的变量根据观测年份间距进行加权平均，然后通过残差拟合得出下一次测量的预测值，再对预测数据进行倾斜、弯曲、扭曲的分析，比对前四次测量数据结果后对古塔变形趋势进行描述。五、模型的建立与求解5.1各层中心位置的求解模型对于该层的中心，已知数据是八边形古塔的各角点，根据正多边形中心的定义，各对角线的交点以及外接圆的圆心是正多边形的中心，但是对于不规则的多边形，其对角线的交点并不重合，同时外切圆的圆心也会偏离。于是不同选点情况下，圆心位置可能存在的偏差，我们采取多次取点计算平均值的方式来消除。首先，引入圆上三点坐标求解中心的模型已知不在同一直线上的三点可以唯一确定一个圆。（如图一）（图一）测量圆形构筑物上的任意3点坐标，设为A（），B（），C（），圆心的坐标为O（），AB的中点为M（），BC的中点为N(),AB与BC的垂直平分线的交点即为圆心由图可知设线段AB的斜率，BC的斜率，其AB的垂直平分线的斜率，其CB的垂直平分线的斜率，其CB的垂直平分线的斜率，则且由此可列出如下方程整理以上两式以矩阵方式表示则求上式解可得（1）公式（1）即为三点坐标法求解中心位置坐标的计算公式。以1986年第一层测量值为例，在多边形上任取三点，两两连线的中垂线角点为圆心，得出外接圆。（如图二）（图二）由于三点的选择存在种情况，所以为消除可能存在的偏差，采用了多次取平均值的办法，对不同角点进行选择，多次平均得到结果。（如图三）（图三）中点坐标为：根据此方法，可以得到该次观察的各层古塔中心和其他三次观察结果：表一：古塔各层中心坐标1986年1996年1566.67,522.70,1.78566.65,522.72，1.782566.72,522.67,7.32566.72,522.68,7.313566.78,522.66,12.75566.78,522.66,12.754566.83,522.64,17.07566.80,522.57,17.075566.89,522.62,21.72566.86,522.65,21.896566.94,522.60,26.23566.95,522.59,26.227566.97,522.58,29.83566.96,522.45,29.838567.16,522.70,33.35566.98,522.57,33.349567.34,522.75,36.85567.00,522.56,36.8410567.06,522.52,40.17567.05,522.36,40.1611567.11,522.48,44.44567.11,522.47,44.4312567.15,522.43,48.71566.84,522.28,48.7013567.18,522.35,52.83567.18,522.38,52.83塔尖2009年2011年1566.77,522.71,1.76566.74,522.72,1.762566.85,522.69,7.30566.79,522.68,7.293566.85,522.72,12.73566.83,522.52,12.724566.86,522.64,17.06566.89,522.72,17.055566.89,522.62,21.70566.89,522.62,21.706567.00,522.04,26.21566.80,522.45,26.207567.04,522.61,29.82567.04,522.60,29.818567.00,522.45,33.33567.03,522.57,33.339567.06,522.45,36.84567.09,522.51,36.8210567.12,522.46,40.16567.16,522.44,40.1411567.13,522.41,44.43567.23,522.38,44.4212567.22,522.36,48.69567.22,522.35,48.6813567.26,522.31,52.81567.35,522.28,52.81塔尖5.2.1，古塔的倾斜分析模型根据测量学定义，取其中三个点使用前方交会的方法，拟合出塔底中心的坐标，塔尖到塔底中心的连线构成倾角。取其与塔顶坐标于底层坐标面上的投影，算出距离，即。（图四）根据题目提供的塔顶坐标数据可以得出塔顶高程，以及塔尖在水平面上的投影与塔底中心间距，即可得出倾斜角。以此类推得到4次观测的塔倾斜情况表二：各年塔尖投影与塔底中心距离及倾斜角1986年1996年2009年2011年塔顶高程55.1255.1255.1255.12中心距离0.740.770.750.74倾斜角0.77°0.80°078°0.77°5.2.2各层中心点和轴线的线性回归分析模型上一问中，我们对古塔的倾斜程度进行了分析，根据测量学的标准，其塔尖到塔基中点的连线构成倾角，然而实际情况中塔存在弯曲，其每一层的中心都偏离该连线，故我们这一问要分析的是各层中点的连线构成一个怎样的线性回归方程。一般称为高斯—马尔柯夫线性模型(k元线性回归模型)，并简记为，，，线性模型考虑的主要问题是：（1）用试验值（样本值）对未知参数和作点估计和假设检验，从而建立y与之间的数量关系；（2）在处对y的值作预测与控制，即对y作区间估计.多项式回归模型设变量X、Y的回归模型为：其中p是已知的，是未知参数，服从正态分布.称为回归多项式.令，多项式回归模型变为多元线性回归模型.以古塔1986年记录的数据为例，将古塔各层中点投影到和坐标面上。各中心点数据投影在坐标：表三：中心点在投影面上坐标层数1234567X566.67566.72566.78566.84566.89566.95566.98Z1.797.3212.7617.0821.7226.2429.84层数8910111213塔尖X566.17567.34567.07567.11567.15567.19567.25Z33.3536.8540.1744.4448.7152.8355.12我们使用matlab拟合得到多元回归系数及图像（matlab源程序附件贴出）。（图五）由该图可以看出1986年古塔的各层中心是向着X轴正向倾斜（即有向东的倾斜），同时我们可以看出在投影面上各层中心基本保持一次线性，弯曲程度轻微，有些许下凹。多项式插值得到的多元线性回归方程为：通过该方程与塔倾斜方程的比较，可以得到塔弯曲的在XOZ投影面上的最大挠度，以及拟合点与倾斜方程的误差平均值表四：回归函数在面上拟合值层数1234567Z实测1.797.3212.7617.0821.7226.2429.84Z拟合2.566.9811.7616.7221.7527.2130.36层数8910111213Z实测33.3536.8540.1744.4448.7152.83Z拟合33.2737.0339.7744.1548.6252.93matlab拟合得到平面上多元回归系数函数图像：（图六）由该图可以看出1986年古塔的各层中心是向着Y轴负向倾斜（即有向南的倾斜），同时我们可以看出在yoz投影面上各层中心存在上凸的弯曲。得到的多元线性回归方程为：通过该方程与塔倾斜方程的比较，可以得到塔弯曲的在YOZ投影面上的最大挠度，以及拟合点与倾斜方程的偏差平均值表五:回归函数在面上拟合值层数1234567Z实测1.797.3212.7617.0821.7226.2429.84Z拟合1.938.0711.9216.7422.0426.6630.30层数8910111213Z实测33.3536.8540.1744.4448.7152.83Z拟合32.3036.8239.2845.2449.4352.37以此，类设变量X、Y的回归模型为：其中p是已知的，是未知参数，服从正态分布.称为回归多项式.令，多项式回归模型变为多元线性回归模型得到四次测量的中心连线多元回归函数表格：表六：各年中心点函数方程1986199620092011在用上述回归方程算得四次测量数据在和面上投影拟合点坐标与倾斜方程得到的残差，筛选最大值（即挠度）与平均值，得到下表：表七：各年挠度与残差平均值年份1986199620092011平均0.390.972.282.8最大0.992.435.095.61平均0.555.0313.383.14最大1.0511.127.069.255.2.3扭转程度分析按照结构力学中扭曲是空间概念的说明，我们分垂直和水平两个方向对古塔的扭曲状况进行讨论，垂直方向上涉及高程Z，即对各层中心点多元线性回归得到的拟合值与实际值进行残差分析，得到扭曲描述量（）。水平方向，参考材料力学中扭转角的计算，对古塔各层间的轴向扭转进行分析，计算得到平面内位移量后通过计算古塔各层高差得到扭转角，对古塔平面扭曲情况进行描述。5.2.3.1垂直方向上的扭曲程度分析通过多项式插值得到回归函数，可以拟合出相应的点坐标，但是实测值必在该曲线呈现离散分布，其程度大小会影响相关系数的大小，同时也说明了塔在垂直方向上的扭曲程度，如下图：（图七）所以我们找到一个对拟合值与实测值的差异评判量：通过对该评判量的统计分析来描述古塔垂直方向上的扭曲程度，统计表格如下：表八：各年扭曲差异量1986199620092011XOZ面0.140.370.780.93YOZ面0.171.681.150.645.2.3.2角点坐标位移残差分析模型通过对古塔八个角点上的坐标分析我们发现，除了偏离中轴回归曲线外，还有古塔单个层面内的扭转变形。以1986年的第一次观察数据为例（图八）我们看投影到水平面上的投影：（图九）可以看到其相对于塔尖中心存在不同程度的转角（图十）计算得到平面内位移量后通过计算古塔各层高差得到扭转角由此可以得到古塔各层相对于下一层的转角，以扭转角来描述该古塔的扭转情况，得出结果如下表表九：各层扭转角（单位：度）年份层数19861996200920111—26.946.947.027.042—36.946.947.027.013—46.936.947.027.044—56.946.687.027.015—66.967.2615.5415.546—713.5713.577.227.227—813.5413.5411.8411.828—913.5613.5711.8411.909—1013.5613.5510.3510.3410—119.769.7610.0910.0711—129.769.7610.0710.0912—1310.0611.5110.1010.075.3.1分析该塔的变形趋势。根据观测年份间距进行加权平均，然后通过残差拟合得出下一次测量的预测值，再用预测数据对倾斜、弯曲、扭曲的情况进行比对描述出各自的变化趋势，并分析各种变形之间的联系，最终得出该塔的变形趋势预测。加权平均公式：若在一组数字中，出现次，出现次，…，出现次，那么叫做、、…、的加权平均数。记作。其中，、、…、分别是、、…、的频数，同时也是它们的权。在预测值中，根据位移残差公式：可计算出位移变化量。预测的下一次测量古塔各层的中心点坐标如下：（五年后观测为例）表十：各层预测中心坐标（单位：m）XYZ1566.79522.711.762566.78522.667.293566.85522.7212.724567.81522.7417.055566.88522.5321.706567.05522.6126.207566.94522.6529.818567.07522.5533.33679567.06522.4136.8210567.13522.4940.1411567.19522.4044.4212567.23522.3548.6813567.27522.3152.81塔顶567.34522.2155.09塔顶坐标数据可以得出塔顶高程，以及塔尖在水平面上的投影与塔底中心间距，即可得出倾斜角所以，塔的倾斜角=0.77°将古塔各层中点投影到和坐标面上，使用matlab拟合得到多元回归系数及图像::用上述回归方程算得四次测量数据在和面上投影拟合点坐标与倾斜方程得到残差，筛选最大值（即挠度）与平均值表十一：预测挠度与差异平均值XOZYOZ平均2.543.26最大5.358.155根据第二问所得的扭转角，依据加权平均数的计算方法，可得出预测值对应的塔的扭转程度。表十二：各层预测扭转角层1—22—33—44—55—66—7扭转角7.037.027.037.0115.547.22层7—88—99—1010—1111—1212—13扭转角11.8311.8710.3410.0810.0810.08六、结果分析与建议6.1问题一：对于角点的的选择，存在种情况，对于不同选产生外界圆圆心的精度影响，我们采用精度公式的转换为了分析三个监测点构形对所求圆心精度的影响，对公式（2）进行转换（如图十一)A、B、C为圆上任意三点，其所对应的边长分别为SA、SB、SC，三边所对圆心角分别为α、β、γ。（图十一）则公式（2）可化简为由余弦定律可得：同理：，因此公式（2）可变换为：三点构形理论分析如图十一，设A、B两点固定，则C点位于圆弧上的什么位置，求解圆心位置的精度最高。如A、B两点固定，则其边长所对应的圆心角为一定值设（3）为了求式（3）的最小值，也即是求的最小值，在γ一定的情况下，分别对α、β求偏导数并令其为0，即可求出极值点。则（同理可得求解以上两式可得α=β或α=β+180°此种情况下取得极小值，即两点一定，第三点位于大小圆弧的中点精度最高。6.2、问题二：汇总古塔倾斜、弯曲、扭曲的四次数据，对古塔的变形做出如下描述6.2.1倾斜：表十二：各年塔尖投影与塔底中心距离及倾斜角1986年1996年2009年2011年塔顶高程55.1255.1255.1255.12中心距离0.740.770.750.74倾斜角0.77°0.80°078°0.77°由上表可以看出古塔的倾斜角的角度值较小，且稳定在之间，并通过查看塔尖投影坐标与塔底坐标，可以得知，塔的倾斜偏向为南偏东之间6.2.2弯曲：表十三：各年挠度与残差平均值年份1986199620092011平均0.390.972.282.8最大0.992.435.095.61平均0.555.0313.383.14最大1.0511.127.069.25由上表可以看出古塔的弯曲状况在1986年份：古塔在平面，古塔的弯曲程度较弱，最大差值（即挠度）为0.99m，且偏移的残差总量较小，所以可以得出结论在面投影上古塔基本保持一次线性。在平面上，古塔中心函数曲线与倾斜方程的最大差值及偏移总量基本与平面保持一致，所以得出结论，此次测量的古塔弯曲程度较弱。在1996年份：古塔在平面上的弯曲程度较1986年的情况严重，无论是平均差值还是最大差值都有所增长，但相较于古塔的体积在可接受范围内，由此可以认定，其偏移的残差总量较小，古塔基本保持一次线性，但在平面上与1986年的数据比较古塔弯曲数据的最大值与平均值都约十倍增长，因此可以认定古塔在平面弯曲程度与平面相比较强古塔弯曲程度较强。在2009年古塔在平面与1996相比偏差较大，最大差值为5.09m，1996年与1986年平均差值为0.58m而2009年增长到1.31m,增长幅度较小。而在平面最大值数据偏大但平均值在可接受范围内。可以判断2009年古塔弯曲程度较弱。2011年观测时间距离2009年仅有两年但是根据观察数据可以看出古塔在平面平均弯曲值增长了约0.52m相当于1986至1996年10年的增长量。而在平面上平均值有所下降但最大值依然超出范围可以判断2011古塔倾斜程度较强。6.3问题三：古塔的变形趋势预测根据前四次测得的各数据对倾斜，弯曲，扭曲有影响的变量由观测年份间距进行加权平均，然后通过残差拟合得出下一次测量的预测值，预测再对数据进行倾斜、弯曲、扭曲的分析。预测得到塔的倾斜角、中心点回归方式、挠度及偏移残差：倾斜角：=0.77°投影面上回归方程：投影面上回归方程：上最大差值（即挠度）为5.35m，平均差值为2.54m，上最大差值（即挠度）为3.26m，平均差值为5.15m扭转角稳定在结合预测数据，古塔的变形趋势为倾斜角度基本保持稳定，即塔顶与塔底中心的连线有一个基本不变的斜率，但是塔身的弯曲情况却在加剧，挠度数据在不断增大，层中心曲线函数与倾斜直线的平均差值也随之增大，而扭曲变形趋势，随着弯曲的加剧，各层中心点反而有趋近回归函数的趋势，所以，扭曲描述量（）呈现出与弯曲变形程度负相关，水平扭转方面塔的轻微扭转依旧在继续，并未产生过多变化。七、模型评价与推广7.1各层中心位置的求解模型的改进在各层中心位置的求解模型中，我们采取的是将古塔的八个角点左边投影到水平面上的方式，用多次取点获得外接圆圆心平均值的方式逼近该层的中心，但是这样忽略了塔的角点坐标是存在高程差异的，故八个角点可以构成的是空间六面体，在该六面体中的体心才是古塔该层的中心位置。两种模型求解得到的中点对比：（图十二）空间六面体体心积分公式：对于该思路的求解可以将八角点够构成的空间六面体进行三重积分，从而得到体心求解。但是古塔每层八角点的高程坐标无稳定变化规律，每次构成同意体面的各角点不能确定，故无法代入积分公式，此思路还有推广改进余地。7.2角点坐标位移残差分析模型的改进在平面扭转分析中对角点坐标位移残差计算的时候，我们假定古塔的各层是稳定刚体，各角点平面转动同时同步进行，但是实际过程中，除了均匀圆杆，任何棱柱在扭转过程中都会产生形变，且变形程度不均匀故对于水平方向上的扭曲分析，可以进一步的对八个角点各自的对每一条棱产生的扭曲进行分析，同时考虑弯曲与