基于空间直线拟合的古塔变形趋势研究_王玥

第２８卷第３期２０１４年５月兰州文理学院学报（自然科学版）

）ＪｏｕｒｎａｌｏｆＬａｎｚｈｏｕＵｎｉｖｅｒｓｉｔｏｆＡｒｔｓａｎｄＳｃｉｅｎｃｅ（ＮａｔｕｒａｌＳｃｉｅｎｃｅｓ　　　　　　　ｙ　

Ｖｏｌ．２８Ｎｏ．３

Ｍａ２０１４ｙ　

（）０９５１２０１４０２２６９９００３０５－－－　　文章编号：２

基于空间直线拟合的古塔变形趋势研究

赵振学２，王江荣２王　玥１，

（兰州石化职业技术学院石油化学工程系，甘肃兰州７１．３００６０；）兰州石化职业技术学院信息处理与控制工程系，甘肃兰州７３０２．０６０

摘　要：采用粒子群算法求出古塔各层的中心点坐标，并对塔层中心进行空间直线拟合．根据拟合直线的方向向量算出塔体每年偏转角度的近似值．对塔层中心的投影点作直线拟合，以投影点在拟合直线上的疏密程度判断塔体局部弯曲程度；以投影点偏离拟合直线的远近及分布判断塔体局部的扭曲程度．研究结果为古塔变形趋势研究提供了新方法．

关键词：几何中心；粒子群算法；直线拟合；倾斜；弯曲；扭曲；偏转角中图分类号：ＯＡ２４　　　文献标识码：

DOI:10.13804/j.cnki.2095-6991.2014.03.009

０　引言

影响古塔变形的因素通常有自重、气温、风力以及地震、飓风等．这些因素会使古塔产生各种变

］２１－

为了保护古塔，需形，诸如倾斜、弯曲、扭曲等［．

，，其中，ｘａｂｃｄｚｅｙ０∈［０，０］０∈［０，０］０∈［０，，区间为中心坐标值的搜索空间．因为塔每层ｆ０］

的观测点（数据点）近似为正八边形的八个顶点，所以这样定义中心是非常合理的．１．１　粒子群算法

［］

ＰＳＯ算法３如下：

要适时对古塔进行观测，了解各种变形量，以便制定必要的保护措施．某古塔已有上千年历史，是我国重点保护文物．管理部门委托测绘公司先后于６年７月、１９９６年８月、２００９年３月和２０１１１９８

年３月对该塔进行了４次观测，得到４次观测数据（见２３年度全国大学生数学建模Ｃ题附０１件）本文先给出确定古塔各层中心位置的通用方．法，并利用这些观测数据和粒子群算法求出各次测量的古塔各层中心坐标，以所得中心坐标和空间直线拟合来分析该塔倾斜、弯曲及扭曲等变形情况，同时分析了该塔的变形趋势．

设在一个Ｄ维搜索空间中，由ｎ个粒子组成…，，其中第ｉ个粒子表的种群Ｘ＝（Ｘ１，Ｘ２，Ｘｎ）

Ｔ

…，代示为一个Ｄ维的向量Ｘｉ＝（ｘｘｘｉ１，ｉ２，ｉＤ），表第ｉ个粒子在Ｄ维搜索空间中的位置，亦代表问题的一个潜在解．根据目标函数即可计算出每第ｉ个粒子的速个粒子位置Ｘｉ对应的适应度值．

Ｔ

…，…，记第ｉ度为Ｖｉ＝（ｖｖｖｉ＝１，ｎ．２，ｉ１，ｉ２，ｉＤ），

为个粒子搜索到的最优位置（即个体最优值ｐｂｅｓｔ）

Ｔ

…，整个粒子群搜索到的最Ｐｐｐｐｉ＝（ｉ１，ｉ２，ｉＤ）；

…，优位置（全局最优值ｇｐｐｔ）为Ｐｂｅｓ１，２，ｇ＝（ｇｇ

１　中心点确定

凸多边形的几何中心有多种确定方法（例如，，坐标平均值法，分形迭代法等）本文按以下方法来确定凸多边形的几何中心：

空间ｎ个点的几何中心：空间中一点到ｎ个点距离的平方和的唯一最小值点为几何中心．，，…，设Ａ１（ｘｚＡ２（ｘｚＡｎ（ｘｙｙ１，１，１）２，２，２）ｎ，

，，中心Ｐ（则ｚｘｚｙｙｎ，ｎ）０，０，０）

ｎ

２

（Ｓ＝∑［ｘｍｉｎｉ－ｘ０）＋

当两个最优解都找到后，每个粒子根据下ｐＤ）．ｇ

式来更新自己的状态．

Ｔ

粒子状态更新操作如下：

）＝ｗ·）ｖｔ＋１ｖｔ＋ｉｉｄ（ｄ（

）］（）［ｒａｎｔｄ＋ｐ１ｉｄ－ｘｉｄ（η

）］，（）［ｒａｎｔｄ　　　　ｐｄ－ｘ２ｉｄ（ｇη

①

））＝ｘ）ｔ＋１ｔｔ＋１．②＋ｖ　　ｘｉｄ（ｉｄ（ｉｄ（

）表示第ｉ个粒子在ｔ＋１次迭其中，ｖｔ＋１ｉｄ（…，代中第ｄ（２，ｄ＝１，Ｄ）维上的速度；ｗ为惯性（）为０，权重，学习因子）ｄｒａｎ１，２为加速常数（ηη

ｉ＝１

２２

ｚ．　　　（ｙｉ－ｙ０）＋（ｉ－ｚ０）］（）１

此外，为使粒子速度不致过～１之间的随机数．

收稿日期：２０１４０１０３．－－

，作者简介：女，甘肃静宁人，王玥（在校学生，主要从事石油化工生产技术研究．１９９４－）

第３期基于空间直线拟合的古塔变形趋势研究王玥等：

３３

，大，可设置速度上下限，即ｖｖｍａ－ｖｍａｉｄ∈［ｘ，ｘ］，边界值）ｖｍａｔ为当前迭ｘ是之前设定的最大速率（代次数．Ｐｔ是粒子迄今为止搜索到的最优位置；

好，则将其作为当前的Ｇｔ．

）按式①和②更新粒子的速度和位置．４

）如果没有满足终止条件（通常为预设的最５

）；，否大迭代次数和适应值下限值）则返回步骤２则，退出算法，得到最优解．

１．２　参数设定及模型求解

，学习因子：维惯性权重ｗ＝０．６；１＝２＝２ηη数为３，种群数为３适应函数最大迭代数为１００，０，

最小值ｍＶｍａｘ＝１；Ｖｍｉｎ＝－１．ｆｉｔ为０．１；ｉｎ

用粒子群算法分别求得１９８６年、１９９６年、这里２００９年和２０１１年所观测的各层的中心坐标（

，利用四舍五入保留四位小数）数据结果见表１．

Ｇｔ是整个粒子群迄今为止搜索到的最优位置．

用于参数估计的ＰＳＯ算法具体实现过程如下：

）初始化粒子群，随机产生所有粒子群的位１

置和速度，并确定粒子的Ｐｔ和Ｇｔ．）对每个粒子，将其适应值与该粒子所经历２

过的最优位置Ｐ如较好，则ｔ的适应值进行比较，将其作为当前的．

）对每个粒子，将其适应值与整个粒子群所３

经历过的最优位置Ｇ如较ｔ的适应值进行比较，

表１　４次观测各层的中心坐标

层数１２３４５６７８９

１９８６年

１９９６年

２００９年

２０１１年

（）（）（）（）５６６．６６４７，５２２．７１０５，１．７８７４５６６．６６５０，５２２．７１０２，１．７８３０５６６．７２６８，５２２．７０１４，１．７６４４５６６．７２７０，５２２．７０１４，１．７６３３）（）（）（）（１９６，５２２．６６８３，７．３２０２５６６．７２０５，５２２．６６７４，７．３１４６５６６．７６４０，５２２．６６９３，７．３０９０５６６．７６４２，５２２．６６９０，７．２９０５５６６．７（）（）（）（）５６６．７７３５，５２２．６２７３，１２．７５５２５６６．７７５１，５２２．６２５６，１２．７５０８５６６．８００１，５２２．６３８４，１２．７３２２５６６．８００４，５２２．６３８８，１２．７２６９（）（）（）（）５６６．８１６１，５２２．５９４４，１７．０７８３５６６．８１８３，５２２．５９２２，１７．０７５１５６６．８２９３，５２２．６１３２，１７．０６９７５６６．８２９７，５２２．６１２７，１７．０５２０（）（）（）（）５６６．８６２１，５２２．５５９１，２１．７２０５５６６．８６４９，５２２．５５６２，２１．７１５９５６６．８６０３，５２２．５８６６，２１．７０９４５６６．８６１０，５２２．５８６０，２１．７０３９（）（）（）（）５６６．９０８４，５２２．５２４４，２６．２３５１５６６．９１１８，５２２．５２１０，２６．２２９５５６６．９４７１，５２２．５３４２，２６．２１１０５６６．９４７８，５２２．５３３４，２６．２０４５（）（）（）（）５６６．９４６８，５２２．５０８１，２９．８３６９５６６．９５０５，５２２．５０４２，２９．８３２２５６６．９７９２，５２２．５１２３，２９．８２４６５６６．９８００，５２２．５１１５，２９．８１７０（）（）（）（）５６６．９８４３，５２２．４９２４，３３．５３０９５６６．９８８４，５２２．４８８０，３３．３４５４５６７．０３０５，５２２．４７９７，３３．３３９９５６７．０３１３，５２２．４７８８，３３．３３６６（）（）（）（）５６７．０２１７，５２２．４７６４，３６．８５４９５６７．０２６５，５２２．４７１４，３６．８４８３５６７．０８１６，５２２．４４６６，３６．８４３７５６７．０８２５，５２２．４４５７，３６．８２２２

）（）（）（）１０（５６７．０５６９，５２２．４６２４，４０．１７２１５６７．０６２０，５２２．４５７２，４０．１６７７５６７．１３７０，５２２．３９３７，４０．１６１１５６７．１３８０，５２２．３９２６，４０．１４４１）（）（）（）１１（５６７．１０４５，５２２．４２３０，４４．４４０９５６７．１１０２，５２２．４１７３，４４．４３５３５６７．１７９９，５２２．３５４６，４４．４３２６５６７．１８１０，５２２．３５３５，４４．４２４９）（）５）（）１２（５６７．１５１８，５２２．３８３６，４８．７１１８５６７．１５７８，５２２．３７７５，４８．７０７４６７．２２２５，５２２．３１６０，４８．６９９７５６７．２２３８，５２２．３１４７，４８．６８３９）（）（）（）１３（５６７．４４５６，５２２．７８６７，５２．８２７９５６７．４５１０，５２２．７８１３，５２．８２４６５６７．２７１２，５２２．２７１５，５２．８１８４５６７．２７２５，５２２．２７０１，５２．８１３１））（）（）（塔尖（５６７．２８０４，５２２．２４８４，５５．０９０８５６７．１７１２，５２２．１３４０，５５．１８２１５６７．３３６，５２２．２１４８，５５．０９１５６７．３３５７，５２２．２１３５，５５．０８７０对于１其中１６年及１９９６年第１３层所缺失的数据可通过拟合插值的办法填补．６年的补充数据为（５６９．９７０１，５２３．９８９８　　说明一点，

），）１１１９，５２．７８３１９９６年的补充数据为（５６９．９７０３，５２３．１１２３，５２．７８７．

２　空间直线拟合

４］

利用表１中的数据采用空间直线拟法［对各

层中心进行直线拟合，结果见图１－图４

．

图２　１６年中心拟合图９９

图１　１９８６

年中心拟合图

兰州文理学院学报（自然科学版）８卷　　　　　　　　　　　　　　　　第２　　　　　　　　　　　　　　　　３４

（）５

００２２．７ｚ＋５９０５５８．ｙ＝－０．

２．２　方向角计算

）与ｘ轴，向量ＯＰ＝（ａ，ｂ，ｃｚ轴的夹解ｙ轴，　　　　则分别为α，γ，β，ｃｏｓα＝

｛

ｘ＝０．ｚ＋５０１１６６．６７５５６，

，ｃｏｓ＝β２２２２２２

＋ｂ＋ｃ＋ｂ＋ｃｃｏｓγ＝

图３　１９８６

年中心拟合图

．２２２

＋ｂ＋ｃ））得拟合直线的方向向量及方向由（２５－（

角，结果如表２所示．

表２　中心拟合直线的方向向量及方向角

年份

方向向量

向量的方向角（αγ）β

））（，，１００７４，１８９．３７５４６°９０．４２４°０．８１０２９°６（０．０１０９，９８－０．）（，，）１９９００８４，１８９．４２１３０°９０．４８１３°０．８１０２９°６（０．０１０１，－０．）（，，）００９０，１８９．３２９６°９０．５１５７°０．８１０２９°２００９（０．０１１７，－０．）（），，００９０，１２０１１（０．０１１７，８９．３２９６°９０．５１５７°０．８１０２９°－０．

约定ｘ轴正方向为正东方向，　　结论：ｙ轴正方向为正北方向．从方向角可以看出古塔从１６年９８

图４　１６年中心拟合图９９

以后出现偏东方向倾斜．平均每年的偏角为，，）（０６５５°０．００２４６°０°．０．０

２．１　空间拟合直线方程

６年的中心直线拟合方程：１９８

ｘ＝０．ｚ＋５０１０６６．６９３０１，

　　　　

００２２．７ｚ＋５７４２８４；ｙ＝－０．

｛｛｛

（）２

３　倾斜角计算

各层中心点与底层中心点连线的倾斜角（线面角，面为ＸＯＹ坐标面）及斜率．

古塔无任何变形时，则各层中心及塔尖在塔底（光滑平面）的投影为同一点（底层中心）．

以塔底中心为基准点来考查各层中心的偏移中心差Δ各层的倾角量，并求各层的高差ΔＨ、Ｄ、如表３．α及倾斜角．

１９９６年的中心直线拟合方程：　　　　

ｘ＝０．ｚ＋５０１０６６．６１４８９，

ｚ＋５００２２．７８４４５５；ｙ＝－０．

（）３

２００９年的中心直线拟合方程：０１１６６．６ｘ＝０．ｚ＋５７５５５，

　　　　

ｚ＋５００２２．７９０５５８；ｙ＝－０．

２０１１年的中心拟合直线方程：

（）４

表３　４个时间段的古塔各层倾斜角及斜率

１９８６年倾斜度

塔层

倾斜率

２　

３　４　５　６　７　８　９　

８０．０６５１８４　８０．０６５１８４　９５．１４７６１４　８０．１５０８８８　７９．７５０７６０　８０．８０７２１４　８１．８０６４５１　８２．０１６１７３　

倾斜角

（单位：度）８９．２８４４２３　８９．２８４９４７　８９．３９７８４４　８９．２８５１８８　８９．２８１６０２　８９．２９０９９３　８９．２９９６５３　８９．３０１４４３　８９．３０７４３５　８９．２９４３４２　８９．２８３８５６　８９．１１９４１４　８９．１９８２３７　

１９９６年倾斜度倾斜率７８．９２５９０　７８．９９０２２　７９．０４７４２　７９．０１１４０　７８．６１１８４　７９．６５３４３　８０．４５０３４　８０．９３４９１　８１．５３７０２　８０．０３６６８　７９．０３０１６　７９．９９７５５　７４．９１８７４　

倾斜角

（单位：度）８９．２７４０９　８９．２７４６８　８９．２７５２１　８９．２７４８８　８９．２７１１９　８９．２８０７２　８９．２８７８８　８９．２９２１１　８９．２９７３３　８９．２８７１０　８９．２７５０５　８９．２８３８　８９．１９８２３　

２００９年倾斜度倾斜率

倾斜角

（单位：度）

２０１１年倾斜度倾斜率１１２．０６１５２　

倾斜角

（单位：度）４８８８９．７２５

１１２．６９２１５９９．４９１５８６　８　１１３．３９８４８５９．４９４７５３　８　１１３．１２９９５９９．４９３５５３　８　１１３．１８６２９４９．４９３８０５　８　８８．３７４１５９　８８．７６４８９７　８３．９６１３５９　８０．２９６４２２　７４．８７０５８１　７４．０７２５６８　７４．７４２８９６　７３．５９１２３４　７２．３８７１３０　

８９．３５１６９６　８９．３５４５４９　８９．３１７６２６　８９．２８６４８４　８９．２３４７８１　８９．２２６５３８　８９．２３３４７４　８９．２２１４８０　８９．２０８５３１　

１１３．５５８５２９９．４９５４６４　８１１２．６６４５６９９．４９１４６１　８１１２．７８３７５４９．４９１９９８　８８８．０８８０６７　８８．６７３８６２　８３．７３８８５２　８０．０５０５８１　７４．６４９１７０　７４．５９２４０９　７４．５２５２７８　７３．４０７１７４　７３．４０４７１４　

８９．３４９５０４８９．３５３８８６８９．３１５８１２８９．２４８２８８８９．２３２５１１８９．２３１９２７８９．２３１２３６８９．２１９５２８２１９８９．５０１

１０２．７２５７６２　８　１１１．１９０６８６　８　１２０．００１８２５　８　１３５．０６０４０８　６　

塔尖７１．４５７６１６　

．ｒｃｔａｎ　　其中：α＝ａ

ＤΔ

数据分析：

垂直方向分析：

所用ＭＡＴ利用Ｋ均值聚类法（ＬＡＢ函数为

［５］）对斜度进行聚类分析．ａｎｓ９年与２０１１ｋｍｅ２００

年的斜度具有相同的聚类结果：第一层至第五层为第一类，第六层至第八层为第二类，第九层至塔从图５可看出聚类结尖为第三类，轮廓图如图５，果非常好．聚类结果表明古塔出现了较明显的弯曲变形，有两明显的拐点（第五层与第六层间有一，拐点，第八层与第九层间有一个拐点）弯曲度较大（从斜度变化可看出）

．

图６　１６

年投影及直线拟合图９８

图７　１９９６

年投影及直线拟合图

图５　２０１１年斜度Ｋ均值聚类

第一层至第十二层１９８６年的斜度聚类结果：为一类，第十三层至塔尖为一类．数据表明古塔整体没有出现较明显的弯曲，在第十二层与第十三层间出现一个拐点，即从第十三层后发生了弯曲．第一层至第七层为第一６年的斜度聚类结果：１９９

，类（含第十二层）第八层至十三层为第二类（除第，十二层）塔尖为第三类．分类表明在第十一层与第十三层间发生了扭曲现象，在第十三层与塔尖间也发生较大弯曲．

图８　２００９

年投影及直线拟合图

４　投影点直线拟合

将各层中心投影到水平面（即ＸＯＹ坐标面）

６－７］

，结果见图６－图９．上，并进行直线拟合［

水平方向分析：各层中心在水平面上的投影点的分布情况能够反映古塔的变形情况．从测量数据及所求中心坐标来看，相邻中心距近似相等．如果无任何变形，则各层的中心投影点重合，如果无扭曲的倾斜，则各中心的投影响点必在一条直线上，

并且可根据投影点在直线上的疏密程度判

图９　２０１年投影及直线拟合图１

断出古塔的弯曲程度．投影点越密，则弯曲程度越，高（即曲率越大）反之，弯曲程度越小（即曲率越小）如果古塔出现扭曲或扭曲性倾斜等，则各层．中心的投影点必不在一条直线上．根据投影点的拟合直线可判断，如果投影点偏离拟合直线较大，则扭曲程度较大．

６年的中心投影拟合图表明古塔在第六１９８

层至第九层处发生了扭曲，但整体上还是相对完好．６年中心投影拟合图显示古塔１０年后的１９９在六层以上发生了较大的扭曲和弯曲，特别是１２层以上，说明在这一阶段，古塔受自然力或人为因素的影响较大，塔的破坏程度严重．０９年及２０古塔受到文物的保护局的抢２０１１年拟合图表明，救，超好的方向发展．

并以投影点在拟合直线上的疏密程度判断塔体局部弯曲程度；以投影点偏离拟合直线的远近及分布判断塔体局部的扭曲程度．该文方法为古塔变形趋势研究提供了新思路和新方法，值得有关部门借鉴．参考文献：

］［］梁海奎．古塔变形测量方法探讨［城市勘测，１Ｊ．２０１１

（）：３１１４．１１３－

［］黄强古．］塔变形监测探讨［测绘与空间地理信息，２Ｊ．

（）：２２２０．３，３６６２１７０１－

［］王江荣．］基于粒子群优化算法的直线拟合及应用［３Ｊ．

（）：工业仪表与自动化装置，５．３４７３２０１７－

［］袭杨．］空间直线拟合的一种方法［齐齐哈尔大学学４Ｊ．

（）：报，２００９，２５２６４６７．－

［］谢中华．５ＴＬＡＢ统计分析与应用：ＭＡ４０个案例分析

［北京：北京航空航天大学出版社，Ｍ］．２０１０．［］王江荣，刘建清．数学建模与数学实验［北京：高６Ｍ］．

等教育出版社，２１．０１

［］张德丰．北京：北京７ＴＬＡＢ数值分析与应用［Ｍ］．ＭＡ

国防工业出版社，２００９．

５　结语

该文采用空间直线拟合算法对塔层中心点进行拟合，根据方向向量求出拟合直线的方向角，进而得到塔体每年偏转角度的近似值．通过对塔层倾斜度的聚类分析，得出塔体发生弯曲的具体位置．另外，对塔层中心的投影点也作了直线拟合，

ＴｏｗｅｒＤｅｆｏｒｍａｔｉｏｎＴｒｅｎｄＲｅｓｅａｒｃｈＢａｓｅｄｏｎＳａｃｅＳｔｒａｉｈｔＬｉｎｅＦｉｔｔｉｎ　　　　　　　　　ｐｇｇ

１２２

，，ＷＡＮＧ　ＹｕｅＺＨＡＺｈｅｘｕｅＷＡＮＧ　ＪｉａｒｏｎＯ　ｎ－ｎｇｇ－

，；（ａｎｚｈｏｕＰｅｔｒｏｃｈｅｍｉｃａｌＣｏｌｌｅｅｏｆＴｅｃｈｎｏｌｏＬａｎｚｈｏｕ７３００６０，Ｃｈｉｎａ１．Ｌ　　　　　ｇｇｙ

，，），Ｌａｎｚｈｏｕ７３００６０ＣｈｉｎａｆｏｒＩｎｆｏｒｍａｔｉｏｎＰｒｏｃｅｓｓｉｎＬａｎｚｈｏｕｅｔｒｏｃｈｅｍｉｃａｌＣｏｌｌｅｅｏｆｔｅｃｈｎｏｌｏ２．ＤｅａｒｔｍｅｎｔｏｆＣｏｎｔｒｏｌＥｎｉｎｅｅｒｉｎ　　　　　　　　　　ｐｇｇｙｇｐｇｇ　

：ＡｂｓｏｉｎｔａｏｄａｅａｃｈｃｏｏｒｄｉｎａｔｅｏｆｔｈｅｃｅｎｔｅｒｔｒａｃｔＵｓｉｎａｒｔｉｃｌｅｓｗａｒｍｏｔｉｍｉｚａｔｉｏｎａｌｏｒｉｔｈｍ，ｔｈｅ　　　　　　　　　　ｐｐｇｇｐｐｇ　

，ｗｅｒｅｓｏｌｖｅｄａｎｄｓａｃｅｌｉｎｅａｒｆｉｔｔｉｎｏｎｆｌｏｏｒｓｏｆｔｈｅｔｏｗｅｒｃｅｎｔｒｅｗａｓｍａｄｅ．Ａｃｃｏｒｄｉｎｔｏｔｈｅｄｉｒｅｃ－　　　　　　　　　　　　　ｐｇｇ　　，ｔｉｏｎｖｅｃｔｏｒｔｏｗｅｒｄｅｆｌｅｃｔｉｏｎａｎｌｅａｒｏｘｉｍａｔｉｏｎｏｆａｅａｒｗａｓｃａｌｃｕｌａｔｅｄ．Ｔｈｅｒｏｅｃｔｉｏｎｏｉｎｔｏｆ　　　　　　　　　　　　ｇｐｐｙｐｊｐ

，ｔｏｗｅｒｃｅｎｔｅｒａｖｅａｓｔｒａｉｈｔｌｉｎｅｆｉｔｔｉｎｔｈｅｌｏｃａｌｃｕｒｖａｔｕｒｅｏｆｔｏｗｅｒｍａｂｅｕｓｄｅｄｂｒｏｅｃｔｅｄｏｉｎｔ　　　　　　　　　　　　　　ｙｇｇｇｊｇｙｐｊｐ　　，ｄｅｎｓｉｔｉｎｆｉｔｔｉｎｌｉｎｅ．Ｂａｓｅｄｏｎｒｏｅｃｔｅｄｏｉｎｔｄｉｓｔａｎｃｅｏｆｄｅｖｉａｔｉｎｆｉｔｔｉｎｌｉｎｅａｎｄｌｏｃａｌｄｉｓｔｒｉｂｕｔｉｏｎ　　　　　　　　　　ｙｇｐｊｐｇｇ　　　　ｔｈｅｔｗｉｓｔｄｅｒｅｅｏｆｔｏｗｅｒｌｏｃａｌｗａｓｄｅｔｅｍｉｎｅｄ．Ｔｈｅｓｔｕｄｒｅｓｕｌｔｓｒｏｖｉｄｅｄａｎｅｗ　ｍｅｔｈｏｄｆｏｒｓｔｕｄｉｎ　　　　　　　　　　　　　ｇｙｐｙｇ　ｔｏｗｅｒｄｅｆｏｒｍａｔｉｏｎｔｅｎｄｅｎｃ．　　ｙ

：；；；；；Ｋｅｗｏｒｄｓｅｏｍｅｔｒｉｃｃｅｎｔｒｅａｒｔｉｃｌｅｓｗａｒｍａｌｏｒｉｔｈｍ；ｌｉｎｅａｒｆｉｔｔｉｎｓｌｏｅｂｅｎｔｄｉｓｔｏｒｔｉｏｎｔｈｅｄｅｆｌｅｃ－　　　　　ｇｐｐｇｇｙ　ｔｉｏｎａｎｌｅ　ｇ