毕业设计（论文）-测量控制网平差系统设计.docx

测量控制网平差系统设计摘 要本软件的是以Visual Basic 6.0为开发环境，用ActiveBar制作主界面，软件应用的关键是控制网数据的输入，软件支持数据输入采用表格形式，当软件得到有关控制网的相关信息后，自动生成相匹配的表格，并且表格以不同颜色显示可编辑（也即必须填写）和不可编辑的单元表格。平差数学模型采用间接平差，平差成果以表格形式显示，成果还包括一些重要的网行参数和平差参数，成果表格以不同颜色显示已知值和平差值，由于受时间的限制，平差软件只实现了简易平差的水准网、导线网、单三角锁，严密平差的单导线（符合导线、闭合导线和无定向导线）平差、导线网平差和高斯正反算、坐标换带、交会计算等。操作简单，便于使用，提高软件的直观性是本软件的主要设计宗旨。关键词：严密平差 高斯正反算 换带计算AbstractThe software is based on the Visual Basic 6.0 development environment. ActiveBar production with the main interface, the software applications, the key is to control network data input, software support, data entry using tabular form, the control software is network related information, Automatic Generation of match form in a different form and color indicates the editing (that is required) and an editing module forms. Adjustment mathematical model using indirect adjustment, the adjustment results in tabular form, the results also include some important parameters of the peace line network parameters poor results in a different form known color values indicate the margin of peace, Due to time constraints, software adjustment only a simple adjustment to achieve the standard networks, cable networks, single triangular locks, Adjustment tight single-wire (with wire coil, and non-directional wire) Adjustment. Traverse Network Adjustment and Gaussian positive and negative counting, coordinates exchange zone, intersection calculation. Simple, easy to use and improve software visualization software is the main design objectives. Key words : positive and negative close Adjustment Gaussian operator for the Calculation目 录引 言1第一章 软件设计平差计算的理论基础31.1水准网简易平差31.2导线测量简易平差31.3单三角形锁的平差理论61.4导线网严密平差81.5高斯约当法101.6高斯投影及换带计算111.7坐标换带公式171.8 交会计算18第二章 本平差软件的简介与开发环境202.1主要平差软件简介202.2本软件的开发环境和运行环境20第三章 测量平差软件操作使用说明223.1测量平差软件安装说明223.2测量平差软件的主要菜单介绍24第四章 平差软件的应用264.1简易平差的应用264.2单导线平差实例334.3测量工具的应用374.3.1高斯投影374.3.2交会计算40结 论41致 谢41参 考 文 献4244引 言测量平差的主要任务是对一系列带有观测误差的观测值，运用概率统计的方法来消除它们之间的不符值，求出未知量的最可靠值，并且评定成果的精度。由于测量仪器的精度不完善和人为因素及外界条件的影响，测量误差总是不可避免的。为了提高成果的质量，处理好这些测量中存在的误差问题，观测值的个数往往要多于确定未知量所必须观测的个数，也就是要进行多余观测。有了多余观测，势必在观测结果之间产生矛盾，测量平差的目的就在于消除这些矛盾而求得观测量的最可靠结果并评定测量成果的精度。测量平差与其他学科一样，试由于生产的需要而产生的，并在生产实践过程中，随着科学技术的进步而发展。18世纪末，在测量学、天文测量学等实践中提出如何消除由于观测误差引起的观测值之间矛盾的问题，即如何从带有误差的观测值中找到观测值的最优值。1794年，年仅17岁的高斯（CF Gauss）首先提除了这个问题的解决方法最小二乘法。他是根据偶然误差的四个特性，并以算术平均值为待求量的最或然值出发，导出了偶然误差的概率分布，给出了在最小二乘原理下求待定量最或然值的计算方法。当时高斯没有正式发表。19世纪初（1801年），天文学家对刚发现的谷神星运行轨道的一段弧长作了一系列观测，后来因故中止了。这就需要根据这些带有误差的观测结果求出该星运行的实际轨道。高斯用自己提出的最小二乘法解决了这个当时很大的难题，对谷神星运行轨道进行了预报，使天文学家及时的又找到了这颗彗星。1809年高斯才在天体运动的理论一书中正式发表了他的方法。在此之前，1806年，勒戎德尔（AMLegendre）发表了决定彗星轨道的新方法一文，从代数观点上也独立的提出了最小二乘法，并定名为最小二乘法。所以，后人称它为高斯勒戎德尔方法。 自19世纪初到20世纪五六十年代的一百多年来，测量平差学者在基于最小二乘原理的平差方法上做了许多研究，提出了一系列解决各类测量问题的平差方法，针对这一时期的计算工具的情况，提出了许多分组解算线性方程组的方法，达到了简化计算的目的。自20世纪70年代开始，随着计算机技术的进步和生产实践中的高精度要求，测量平差得到了很大发展，主要表现在：从单纯的偶然误差到包含系统误差和粗差。21947年，铁斯特拉（TMTienstra）提出了相关观测值的平差理论。31969年，克拉鲁（T.Krarup）提出最小二乘虑波，也称为拟和推估法。41962年，迈赛尔（P meissl）提出针对非满秩平差问题的内制约束平差原理。5对系统误差的特性、传播、检验、分析的理论的发展6.测量中粗差理论的研究。第一章 软件设计平差计算的理论基础1.1水准网简易平差一、水准测量原理：是借助仪器提供的水平视线，配合水准尺测定地面上两点间的高差，然后根据已知点的高程来推求未知点的高程。高差闭合差： 1h2h1BM1fh=h测-（H始-H终）2高差闭合差的允许值： .fh=±40Lmm或±10n mmh3式中 L水准路线的长度， 以km计；h4 n测站数。图1 闭合水准路线高差闭合差的调整 在同一条水准路线上，可以认为各测站条件大致相同，故闭合差的调整可以反号按距离往测2（或测站数）成正比分配各测段的高差中，即 i=-fh/L*Li1BM1i=-fh/n*n返测式中 总和的简写符号； BM2n测站数总和；图2 附合水准路线 ni 某测段的测站数；i 某测段的高差改正数； Li 某测段的路线长度；L水准路线总长度。1.2导线测量简易平差导线测量：是在地面上按照一定的要求选定一系列的点（导线点），将相邻点连成直线而构成折线形，依次测定各折线边（导线边）的长度和各转折角（导线角）；根据起算数据，推算各边的方位角，从而求出各导线点的坐标。导线网中导线点和路线的编号规则：与水准路线基本相同。1AB2BA3图3 闭合导线ABCD31DAC2B图4 附合导线1. 角度闭合差的计算与调整 由平面几何学可知，n边形的内角和应为 理=（n-2）*180° 由于测角误差的存在，所观测的n边形内角和测不等于理论上的内角和理，其差值称为角度闭合差，用f表示，即 f=测-理角度闭合差的允许值为：f允=±60n式中 n内角个数2. 导线的坐标方位角的推算相同前进方向的相邻两条导线边，起坐标方位角的关系为 前=后+左±180前=后-右±180式中 前、后相邻导线前、后的坐标方位角； 左、右相邻导线边左、右转折角。3. 相邻导线点之间的坐标增量计算所谓坐标增量，就是两导线点坐标 X Y B X AB A Y 值之差。也就是从一个导线点到另 一个导线点的坐标增加值。坐标增量 有纵坐标增量X与横坐标增量Y 随着导线边所处的象限不同，坐标增 量有正负之分。 坐标增量是用导线边的边长与坐标方 位角计算的。如图所示，导线边AB的 边长与方位角分别为DAB、AB。 从A点到B点的纵、横坐标增量为图5 坐标推算 XAB=DABcosAB YAB=DABsinAB4. 坐标增量闭合差的计算与调整闭合导线的纵、横坐标增量代数和在理论上应该分别等于零。即X理=0，Y理=0，因为坐标增量就是两导线点的坐标之差，而闭合导线是从已知点开始，最后又闭合到同一点上，那么同一点坐标增量应该等于零，其与理论值不相符的数值称为纵、横坐标增量的闭合差： 闭合导线： fx=X fy=Y 附合导线： fx=X-（X终 -Y始） fy=Y-（X终 -Y始） 由于fx、fy的存在，使得闭合导线起算点的起算坐标与通过坐标增量计算出的坐标两者不相等。fx、fy值的大小反映了起算坐标所确定的点B与计算坐标所确定的点B1之间的距离大小。两点之间的距离可用下式计算： fD=f2x+f2y fD称为导线全长绝对闭合差。单从fD的大小还不能说明导线测量的精度如何，还不能判断是否达到规范的要求。在导线测量中，用fD与导线边全长之比fD/D来衡量测量的精度。 用表示其比值的大小，并且要把比值的结果化算为分子为1的分数形式表示： K=fD/D=1/D/fD K值称之为导线的全长相对闭合差。各种等级的导线测量对K值的限差要求不同。 如果K值达到精度要求，应进行闭合差fx、fy的调整，以提高导线点的精度。 坐标增量闭合差fx、fy的调整方法是按照与导线边长成正比例的原则，反符号分配到各坐标增量中。通常是先求出纵、横坐标增量的改正数xij、yij，然后加到相应的坐标增量中。 xij=-fx/D*Dijyij=-fy/D*Dij 式中 xij、yij 导线点i到j的纵、横坐标增量改正数; D 导线边长总和； Dij 导线边ij的边长1.3单三角形锁的平差理论1. 角度闭合差的计算与调整各三角形内角和应该分别等于180°，但是由于所测水平角存在着不可避免的误差，使得各三角形内角和不等于180°，其差值即角度闭合差用f表示 f i =ai+bi+ci-180°式中 i=1，2，.n。测角中误差的计算公式为 m=±fi*fi/3n若f i小于限差要求，则将其反符号平均分配给各内角观测值，并计算出第一次改正后的各内角角值： ai =ai-fi/3 bi =bi-fi/3 ci =ci-fi/3 经过第一次角度闭合差调整后的各三角形，其内角和应等于180°。作为计算较核。 2. 基线闭合差的计算与调整根据起始基线边长D0及经过第一次角度闭合差平差后的各传距角ai、bi，应用正弦定律即可推算出传距边的长度，D1=D0sina1/sinb1D2=D1sina2/sinb2=D0sina1*sina2/sinb1*sinb2=D0/. Dn= D0sina1*sina2sinan/sinb1*sinb2sinbn= D0/ W基=（1-Dn/ D0）基线闭合差的允许值为 W基允=2m 式中 m相应等级规定的测角中误差 所以： Vai =-Vbi= -W基/ 经过第二次角度闭合差调整后的各三角形内角为 ai=ai+ Vai bi=bi+ Vbi ci=ci 调整后 ai+ bi+ ci=180°D0/ Dn=13.各三角形边长的计算 根据起始基线边D0及第二次调整后的各三角形内角，按三角形顺序，应用正弦定律逐一推算各三角形的边长。 D23= D0/sinb1*sina1D13= D0/sinb1*sinc15. 各三角形坐标推算各三角点坐标的计算可根据坐标起算数据、第二次调整后的各三角形内角及推算出的边长，按照闭合导线12467531计算各点坐标。因为，经过平差后的三角锁，其角度闭合差、坐标增量闭合差均已满足理论要求。所以，不存在角度闭合差、坐标增量闭合差的计算与调整的问题。1.4导线网严密平差间接平差法：是通过t个独立参数，将每个观测量分别表达成这t个参数的函数，建立函数模型，按最小二乘原理；用求自由极值的方法解出参数得最或然值，从而求得各观测量的平差值。间接平差函数模型： （1） 间接平差随机模型： （2）用参数的估值和的估值V代入，可以写成： （3）表达了参数估值与观测值改正数V之间的函数关系，称为误差方程。由于误差方程个数为n，待求量和V的总数为n+t，而n<n+t，（3）式具有无穷多组解，但可按最小二乘原理，在VTPV=min下求的唯一解。下面导出间接平差的计算公式。一. 间接平差的基础方程及其解 设平差中有n个观测值，已知其协因数阵Q=P-1，必要观测数为t，选定t个独立量为参数，其估值为，观测值L与改正数V之和,称为观测值平差值。按具体平差问题，可列n个平差值方程为： (i=1,2,，n)令： 则平差值方程的矩阵形式为 （4）令 （5）式中X0为参数的充分近似值，于是得误差方程为： （6） 按最小二乘原理，上式的必须满足VTPV=min的要求，因为个参数为独立量，故可按数学上求函数自由极值的方法，得 转置后得 BTPV=0 （7） 以上所得的(6)和(7)式中的待求量是n个V和t个，而方程个数也是n+t个，有唯一解，称此两式为间接平差的基础方程。 解此基础方程，一般是将(6)代入(7)式，以便先消去V，得 BTPBBTPL=0 令 ，上式可简写成： NbbW=0 （8）式中系数阵Nbb 为满秩，即R(Nbb)=t , 有唯一解，上式称为间接平差的法方程 。解之，得 （9）或 =(BTPB)-1BTPl (10)将解出的代入误差方程(6)，即可求得改正数V,从而平差结果为： ， 特别地，当P为对称阵时，即观测值间互独立，则法方程(8)的纯量形式为.按间接平差法求平差值得计算步骤（1）.根据平差问题的性质，选择t个独立量作为参数；（2）.将每一个观测量的平差值分别表达成所选参数的函数，若函数非线性要将其现性化，列出误差方程；（3）.由误差方程系数B和自由项l组成法方程，法方程个数等于参数的个数t；（4）. 解算法方程，求出参数，计算参数的平差值 ；（5）.由误差方程计算，求出观测量平差值 1.5高斯约当法矩阵求逆采用全选主元高斯约当（Gauss-Jordan）消去法。高斯约当（全选主元）求逆的步骤如下。首先，对于k从0到n-1作如下几步：a) 从第k行、第k列开始的右下角子阵中选取绝对值最大的元素，并记住此元素所在的行列号，再通过行交换和列交换将它交换到主元素位置上。这一步称为全选主元。b)c)d) e)最后，根据在全选主元过程所记录的行、列交换信息进行恢复，恢复原则如下：在全选主元过程中，先交换的行、列后进行恢复；原来的行（列）交换用列（行）交换来恢复。矩阵的行列式值也是用全选主元高斯约当（Gauss-Jordan）消去法。先将矩阵变为上三角矩阵，通过计算对角线元素的乘积可得。矩阵的其它计算比较简单，在这里不再说明。1.6高斯投影及换带计算高斯投影的主要内容就是导出高斯平面坐标（x,y）与大地坐标（L，B）的相互关系式。关系式分两类：第一类称高斯投影正算公式，亦即由L，B求x,y；第二类称高斯投影反算公式，亦即由x,y求L，B。.高斯投影坐标正算公式 高斯投影必须满足以下三个条件： 中央子午线投影后为直线；中央子午线投影后长度不变；投影具有正形性质，即正形投影条件。由第一个条件可知，由于地球椭球体是一个旋转体，所以高斯投影必然有这样一个性质，即中央子午线东西两侧的投影必然对称于中央子午线。具体地说，比如在椭球面上有对称于中央子午线的两点P1，P2 ，它们的大地坐标分别为（l，B）及（-l，B），式中l为椭球面上有对称于中央子午线的经度与中央子午线的经差，P点在中央子午线之东，l为正，在西则为负，则投影后的平面坐标一定为P1（x,y）和P2（x,-y）。这就是说，在所求的投影公式中，当B=常数，l以-l代替时，x值不变号，而y值则变号，。又由于高斯投影是按带投影的，在每带内经差l是不大的，l/是一个微小量，所以可将式中的函数展开为经差l的幂级数，它可以写成如下形式： x=m0+ m2l2+m4l4+ y=m1+m3l3+m5l5 （1）式中m0，m1，m2，是待定系数，它们都是纬度B的函数。由第三个条件：将式（1）分别对l和q求偏导数代入，得m1+3m3l2+5m5l4+=dm0/dq+dm2/dq*l2+dm4/dq*l4+2m2l+4m4l3+=-dm1/dq*l-dm/dq*l3- （2）为使上面两式两边相等，其必要而充实的条件是l的同次幂的系数相等，因而有m1=dm0/dqm2=-1/2dm1/dq （3）m3=-1/3dm2/dq为要最终的求出待定系数，显示矛盾的焦点在于求得导数dm0/dq。为此，首先要确定m0的表达式。 由第二条件可知，位于中央子午线上的点，投影后的纵坐标x应该等于投影前从赤道量至该点的子午线弧长，即在（1）第一式中，当l=0时， x=m0=X （4）式中X为自赤道量起的子午线弧长。顾及中子午线弧长微分公式dX/dB=M和式dB/dq=NcosB/M，于是得 dm0/dq=dm0/dB*dB/dq=dX/dB*NcosB/M=M*NcosB/M=NcosB （5）故 m1= NcosB=c/VcosB （6）其次求dm1/dq，由上式对q求偏导数得 dm1/dq= dm1/dB*dB/dq= d/dB（c/VcosB）dB/dq =（-c/V2 dV/dBcosB- c/VsinB）顾及 dV/dB=-1/V2t于是得出 dm1/dq=-c/V2（-1/V2t）cosB- c/VsinB* c/VcosB/ c/V3 = c/V3sinB（2- t2） V2cosB=- c/VsinBcosB于是 m2=N/2sinBcosB （7）依次求导，并依次代入（3）式右可得m3 ，m4，各值： m3= N/bcos3B(1- t2+2)m4=N/24sinBcos3B(5-t2+92) （8） m5= N/120cos5B(5- 18t2+ t4) 将上面已经求出的各个确定系数m5，带入（1）式，并略去2l5及l6以上各项，最后得出高斯投影坐标正算公式如下： x=X+N/22* sinBcosB l2+ N/244sinBcos3B(5- t2+92) *l4 （9） y= N/cosB l+ N/63cos3B(1- t2+2) l3 + N/1205cos5B(5- 18t2+ t4) l5当l3.5°时，公式换算的 精度为±0.1 m。欲要换算精确至0.001 m的坐标公式，可将上式继续扩充，现直接写出如下：x=X+N/22* sinBcosB l2+ N/244sinBcos3B(5- t2+92+44) *l4+ N/7206* sinBcos5B（61- 58t2+ t4）l6 （10）y= N/cosB l+ N/63cos3B(1- t2+2) l3+ N/1205cos5B(5- 18t2+ t4+142-582 t2) l5 （9），（10） 两式即为高斯投影坐标正算公式，他们就是式中的F1和F2的具体形式。.高斯投影坐标反算公式高斯投影反算时，原面是高斯面，投影面是椭球面，相应地可以写出如下投影方程： B=1（x,y） l=2（x,y） （11）同正算一样，对投影函数1，2提出以下三个条件：x坐标轴投影成中央子午线，是投影对称轴；x轴上的长度投影保持不变；正形投影条件，即高斯面上的角度投影到椭球面上后角度没有变形，仍然相等。反算公式的 推导方法和正算公式相类似，它的基本思想是，首先根据x计算纵坐标在椭球面上的投影的垂足纬度Bf,接着按Bf计算（Bf-B）及经差l，最后得到 B=Bf-（Bf-B） L=L0+l （12）为了简化差值（Bf-B）及l的计算公式的推导，我们还是借助于等量纬度q。由于高斯投影区域不大，其中y值和椭球半径相比也是不大的，因此可将（q，l）展开为的幂级数。顾及到上述第一个条件，即投影对轴子午线对称的要求，q应是y的奇函数，因此q= n0+ n2y2+n4y4+l= n1+ n3y3+n5y5+ (13)由第三个条件：将（13）式分别对x和y求偏导数；并考虑到上述第二个条件，当y=0时，点处在中央子午线上，亦即l=0,x=X,1(x)=qf于是可求出待定点系数n，把它们代入上式，则得q= qf-y2 （d2q/ d X2）f+ y4/24（d4q/ d X4）f-l= y（dq/ d X）f- y3/6（d3q/ d X3）f+ （14）式中qf和（dq/ d X）f均是对应于垂足纬度Bf的数值。至此还只是求得了等量纬度q，为最终求得大地纬度B，还需进一步化算。由式可知，q和B必有固定的函数关系，今设 B=f（q），Bf= f（qf） （15） B=f（q）= f（qf+ dq） （16） dq=q- qf式中按台劳级数展开则有 B= f（qf）+ （dB/dq）dq+1/2（d2B/dq2）+1/6（d3B/dq3）dq3+因此 （17） -（Bf-B）= B- Bf=（dB/dq）（q- qf）+1/2（d2B/dq2）（q- qf）2 +1/6（d3B/dq3）（q- qf）3+ （18）为了求得高斯坐标的实用反算公式，我们必须求出（14）式和（18）式中的各阶导数值dnq/ d Xn和dnB/ d qn。 首先求导数值dnq/ d Xn。由式知一阶导数 dq/ d x= dq/ d X=1/NfcosBf （19）其他各阶导数可依上式逐次对X求导得到，因此有 d 2q/ d X2=tf/ N2fcosBf d3q/ d X3=1/N3fcosBf(1- 2tf2+f2) （20）d4q/ d X4=1/N4fcosBf(5+ 6t2+2f-46f)其次求导数dnB/ d qn。由（5）式易知一阶导数 （dB/ d q）f= （NcosB/M）=V2 fcosBf （21）其他各阶导数可依上式对q逐次求导得到，其结果是 （ d2B/ d q2）f=-sinB f cosBf(1+4f2+3f4) （ d2B/ d q2）f= -cos3Bf(1- tf2+5f2-13f2 tf2)(22) 当将式中有关代入到（14）式第一式后，又可得到（q- qf）=- tf cosBf/2N2f* y2+ tfsec Bf/24 N4f(5+ 6t2+2f-46f) y4- tfsec Bf/720 N6f(61+ 180t 2 f +120t 4 f+464f-482f t 2 f) y6 对上式分别平方和三次方得（q- qf）2= t2fsec Bf/4 N4f* y4-t2fsec 2Bf/24 N6f (5+ 6t2 f +2f-46f)y6 （q- qf）3= t3fsec 3B/8 N6f y6 (23)将（21）（23）式代入（18）式，将（19），（20）式代入（14）式的第二式，并略去4fy4和6fy6次项，经整理得B= Bf-tf/2MfNf* y2+ tf/24MfN3f(5+ 3t2 f +2f-92f t2 f)y4 -tf/720 MfN5f(61+ 90t2 f + 45t2 f) y6L=1/Nf cosBf y-1/6N3f cosBf(1+ 2t2 f +2f)y3+ (24) 1/120N5f cosBf(5+ 28t2 f +24t2 f+62f-82f t2 f)y5当13.5°时，上式换算精度为0.0001。欲要换算精确至0.01时， 可按上式简化成 B= Bf-tf/2MfNf* y2+ tf/24MfN3f(5+ 3t2 f +2f-92f t2 f)y4L=1/Nf cosBf y-1/6N3f cosBf(1+ 2t2 f +2f)y3+ (25) 1/120N5f cosBf(5+ 28t2 f +24t2 f)y5上式即为高斯投影坐标反算公式，它们就是（11）式中的1和2的具体形式。1.7坐标换带公式坐标换带公式坐标的换带计算有两种方法:直接换算法和间接换算法。直接换算法就是推导出相邻带间的换算关系式,并编制出专门的计算用表,然后直接把一点的坐标换算成相邻带坐标系的坐标。直接换算法有几种不同的公式和相应的计算用表,它们的精度有所不同。可见,直接法需要人工查表,效率较低,不适合大量点的坐标换算,在电算相当发达的今天,经常采用的是间接换算法。如已知p点在带的坐标x1、y1,首先用高斯投影反算公式求得其大地坐标B、L(B为大地纬度,L为大地经度),然后计算P点与带的中央子午线L中的经度差,即=L-L中,最后把B、l代入高斯投影正算公式求得P点在带的坐标值xP、y2可见,在这种方法中大地坐标B、L起到了间接的桥梁作用。间接法理上是严密的,精度也最高。下面就介绍间接换算法的有关公式:(1)高斯投影反算公式求P点的大地纬度B(单位:度):B=Bf-1+2ftf90(y/Nf)2-7.5(5+3t2f+2f-92ft2f)(y/Nf)4+0.25(61+90t2f+45t4f)(y/Nf)6 （26）求P点距带中央子午线L中的经差l(单位:度):=1/cosBf(180(y/Nf)-30(1+2t2f+2f)(y/Nf)3+1.5(5+28t2f+24t4f+62f+82ft2f)(y/Nf)5) (27)公式26和公式27中:f=cosBf其中=0.08209446887257f=/ 其中2=0.00673950181947tf=tgBf =6399596.65199()现在,以上所有公式中只有一个未知数Bf,它叫做底点纬度或垂足纬度,是纵坐标x1(以千公里为单位)的一元函数, Bf的单位为“度”。公式如下:当x1<3000公里时:Bf=9.04369066313x1-0.00000049618x21-0.00075325505x31-0.00000084330x41-0.00000426157x51-0.00000010150x61 （28）当3000公里<x1<6000公里时:Bf=27.11162289465+9.02483657729(x1-3)-0.00579850656(x1-3)2-0.00043540029(x1-3)3-0.00004858357(x1-3)4-0.000