毕业设计（论文）-沉降监测网拟稳平差及程序设计.doc

本科毕业论文本科毕业论文 沉降监测网拟稳平差及程序设计 SUBSIDENCE MONITORING NETWORK QUASI-STABLE ADJUSTMENT AND PROGRAM DESIGN 学院（部）： 测绘学院 专业班级： 测绘 11-1 学生姓名： 指导教师： 年 月 日 沉降监测网拟稳平差及程序设计 摘要 本文简单的介绍了沉降监测网的起因、经过、形成，以及对沉降观测网数 据处理的方法。拟稳平差在沉降观测网的数据处理应用，使用 C 语言对拟稳平 差进行编程，以函数的形式实现拟稳平差中的误差方程、组建法方程以及矩阵 求逆等功能。本文还介绍了有关拟稳点的判断和选取的内容，最后详细的说明 了运行程序的一些必要步骤和上机操作问题。 关键词：沉降监测网，关键词：沉降监测网，C C 语言，拟稳平差，程序设计语言，拟稳平差，程序设计 SUBSIDENCE MONITORING NETWORK QUASI-STABLE ADJUSTMENT AND PROGRAM DESIGN ABSTRACTABSTRACT This paper simply introduced the cause, process, formation of subsidence monitoring network, as well as to the settlement observation data processing method.The quasi-stable adjustment in settlement observation data processing applications, the use of C language programming of quasi-stable adjustment, in the form of function implementation error equation, a method of quasi-stable adjustment equation and matrix inversion, and other functions.This paper also introduced the quasi steady point of judgment and choice of content, finally illustrates in detail some steps necessary to run the program and computer operation. KEYWORDS：Settlement Monitoring Network，C Programming Language，Software，Quasi-Stable Adjustment 目录 ABSTRACT.3 1 沉降监测网概述.6 1.1 沉降监测网的形成.6 2 地下开采引起的地表移动和破坏.6 2.1 地表移动的形式.6 2.2 地表移动盆地.7 3 沉降监测网的数据来源.7 3.1 开采沉陷的观测工作.7 3.2 本次实测数据说明.8 4 拟稳平差程序概述.9 4.1 拟稳平差的基本原理.9 4.2 拟稳平差程序编译的相关函数.10 4.3 函数功能说明及其源代码.11 4.4 原始数据文件格式设计.30 4.5 拟稳点信息的读取与存储.31 4.6 拟稳点的判定与选取.31 4.6.1 拟稳点的判定方法.31 4.6.2 拟稳点的选择原则.32 5.拟稳平差程序的运行.33 5.1 原始数据的处理.33 5.2 程序运行的步骤.34 5.3 程序运行结果显示.35 绪论 中国的经济发展随着改革开放而腾飞，而经济的发展离不开对各种能源的 需求。中国经济从 2003 年起进入了重化工阶段，对能源的需求快速增长，尤其 是对煤矿的需求大大增加。因此从 2003 年起到 2013 年这十年，是煤炭行业黄 金的十年。煤矿企业挖走了大部分的煤，但是造成了开采沉陷留下了大量的采 空区。 这些采空区对国民经济的发展有着不可小觑的阻碍。其中尤其对国家的基 础设施影响极大，例如：在交通中的影响，铁路和国家高速公路不得穿过采空 区上方。为了监测采空区和采区对周边建筑、环境和其它地物的影响，这时要 布设沉降观测网，用以监测和评价采区和采空区对环境的影响。本文正是基于 这一理论的基础上，研究拟稳平差在沉降观测网中的应用以及程序设计。 1 沉降监测网概述沉降监测网概述 1.1 沉降监测网的形成 有用矿物被采出以后，开采区域周围的岩体的原始应力平衡状态受到破坏， 应力重新分布，到达新的平衡。在此过程中，使岩层和地表产生连续的移动、 变形和非连续的破坏（开裂、冒落等） ，这种现象称为“开采沉陷” 。 岩体本身是一种非常复杂的介质。它不仅是由各种不同性质的岩层组成， 而且还由于各种地质作用（如褶皱、断层、开裂、火成岩、侵入，陷落柱等） 而产生了大量的不连续面。岩体在受到各种不同开采方法影响时，产生的开采 沉陷是一个在时间上和空间上都是非常复杂的过程。在时间上来说，在移动的 过程中，开采沉陷的形式和大小在不同的时间上是不同的，也就是说，此时的 开采沉陷是“动态的” ；随着时间的推移，开采沉陷的形式和大小逐渐趋向于稳 定，开采沉陷变成“静态的”或“最终的” 。从空间上来说，若地下开采的范围 较小、开采的矿物的埋藏深度较小，则开采沉陷波及的范围往往只局限于开采 区域的周围的岩体；若开采范围较大、开采矿物的埋藏深度较大，则开采沉陷 波及的范围就会从岩体发展到地表，引起“地表移动” 。 沉降监测网是为了监测“地表移动”的表面上指定位置的点的沉降变化而 建立的一种水准网。通过实地多次观测这些指定位置的点的高程和平面坐标， 确定这些点的空间位置坐标，进行数据统计分析，观察这些点的沉降情况、沉 降速度、沉降变化等，判断并预测这些点位在空间中的变化趋势。 2 地下开采引起的地表移动和破坏地下开采引起的地表移动和破坏 2.1 地表移动的形式地表移动的形式 所谓地表移动，是指在采空区面积扩大到一定范围后，岩层移动发展到地 表，是地表产生移动和变形，在地表沉陷的研究中称这一过程和现象为地表移 动。开采引起的地表移动过程，受多种地质采矿因素的影响，因此，随开采深 度、开采厚度、采煤方法及煤层产状等因素的不同，地表移动和破坏的形式也 不完全相同。在采深和采厚的比值比较大时，地表的移动和变形在空间和时间 上是连续的、渐变的，具有明显的规律性。当采深和采厚的比值较小（一般小 于 30）或具有较大的地质构造时，地表的移动和变形在空间和时间上将是不连 续的，移动和变形的分布没有严格的规律性，地表可能出现较大的裂缝或坍塌 坑。地表移动和破坏的形式，归纳起来有：地表移动盆地、裂缝及台阶、坍塌 坑。 本程序是针对地表移动盆地而言。 2.2 地表移动盆地地表移动盆地 在开采影响波及到地表以后，受采动影响的地表从原有标高向下沉降，从 而在采空区上方地表形成一个比采空区面积还大的沉陷区域。这种地表沉陷区 称为地表移动盆地，或称下沉盆地。在地表移动盆地形成的过程中，改变了地 表原有的形态，引起了高低、坡度及水平位置的变化。 地表移动盆地的形成：地表移动盆地是在工作面的推进过程中逐渐形成的。 一般是当回采工作面自开切眼开始向前推进的距离相当于 1/41/2H0（H0 为平 均采深）时，开采影响即波及到地表，引起地表下沉。然后，随着工作面继续 向前推进，地表的影响范围不断扩大，下沉值不断增加，在地表就形成一个比 开采范围大得多的下沉盆地。 3 沉降监测网的数据来源沉降监测网的数据来源 3.1 开采沉陷的观测工作开采沉陷的观测工作 要保护井巷、建筑物、水体及铁路等。是它们免受或减少开采的有害影响， 减少地下资源的损失，必须研究地下开采引起的岩层与地表移动规律。岩层与 地表移动的过程十分复杂，它是许多地质采矿因素综合影响的结果。认识岩层 与地表移动这一复杂过程，目前的主要方法是实地观测。通过观测获得大量的 第一性资料，然后对这些资料进行综合分析，找出各种因素对移动过程的影响 规律。 为了进行实地观测，必须在开采进行以前，在选定的地点设置开采沉陷观 测站。简称观测站。所谓观测站，是指在开采影响范围内的地表、岩层内部或 其它研究对象上，按一定要求设置的一系列互相联系的观测点。在采动过程中， 根据需要定期观测这些测点的空间位置及其相对位置的变化，以确定各测点的 位移和点间的相对移动，从而掌握开采沉陷的规律。 开采沉陷观测站的类型有以下几种： 1 .按观测站设置的地点不同可分为： （1）地表移动观测站：测点布设在地表。主要研究地表移动和变形的规律； （2）岩层内部观测站：测点一般布设在井下巷道或岩层内部的钻孔中，用 于研究岩层内部的移动和变形规律； （3）专门观测站：为了某一个特定的目的所设立的观测站，如建筑物观测 站、铁路观测站、边坡移动观测站等。 2. 普通观测站和短期观测站 普通观测站的观测时间较长（一般在一年以上） ，它是在地表移动的开始到 移动结束的整个过程中定期进行观测，主要用于研究地表移动和变形的规律。 短期观测站观测时间较短（几个月到一年） ，只在移动过程中的某个阶段进行观 测，而对观测资料的处理，求出一些近似的移动参数，如最大下沉速度、走向 移动角等。短期观测站只是在急需开采沉陷资料的情况下才采用。 3. 按步站形式的不同分为： （1）网状观测站：在产状复杂的煤层或在建筑物密集的地区开采时，可考 虑多布设一下测点，组成网格状观测站。网格状观测站可以对整个采动影响范 围进行观测，所得资料比较全面、准确，但测点数目较多，野外观测和室内成 果整理工作量大，且受地形、地物条件的限制，所以只在研究专门问题时采用。 （2）剖面线状观测站：是指在沿移动盆地主断面的方向上，将观测点布设 成直线的观测站。剖面线状观测站通常由两条互相垂直且相交的观测线所组成。 沿走向主断面布设的观测线称为走向观测线，沿倾斜主断面布设的观测线称为 倾斜观测线。 图 3-1 观测站的布置形式示意图 3.2 本次实测数据说明本次实测数据说明 本程序处理的原始数据包括各个点的平面坐标和各个点之间的观测高差， 其中点位的平面坐标用于计算水准路线的长度，用以定权，原始数据以 excel 表格形式如下图（3-2）和（3-3）所示 图 3-2 点位的平面坐标 图 3-3 各点的观测高差 4 拟稳平差程序概述拟稳平差程序概述 4.1 拟稳平差的基本原理拟稳平差的基本原理 自由网拟稳平差公式 设误差方程式为 （式 4-1）lAXV 式中： 为 n 维自由项向量；V 是 n 维残差向量；A 为系数矩阵；X 是 tltn 维参数向量，这里 X 是高程点近似值的平差改正数。由上式（4-1）组成法方程 式； （式 4-2）0)(PlAXPAA TT 式中，P 为观测值的权矩阵，对称正定。 水准网平差中至少有一个已知高程点。假如网中没有已知高程点，误差方 程式的系数阵 A 的秩就会小于 t，法方程的系数阵也会降秩。拟稳平差将全部 高程点分为拟稳点和非拟稳点两部分，仅在拟稳点参数的平方和最小的约束下 求法方程式的解。设水准网中共有 t 个点，其中 s 个点为拟稳点，为拟稳点 S X 参数解向量，约束条件为 （式 4-3）min S T SX X 我们采用虚拟观测值法求解。设 G、G2都是 t 维向量， sss GT 1 . 11 . 1102 t T gggG 其中（i=0,1,t-1）为 i g i 号点位拟稳点（式 4-5） s gi 1 i 号点为非拟稳点 0 i g 拟稳平差参数解及其权逆阵的公式为 （式 4-6）PlAGGPAAX TTT1 22 )( （式 4-7） TTT X GGGGPAAQ 1 22 )( 单位权中误差计算公式为 （式 4-8） 1 tn pvv 4.2 拟稳平差程序编译的相关函数拟稳平差程序编译的相关函数 本程序使用 C 语言作为编程语言，有以下函数组成： 1)MyBreak 信息提示函数 2)ij 对称矩阵下标计算函数 3)PrintM 数组输出函数 4)PrintM2 对称矩阵输出函数 5)PrintEquation 线性方程组输出函数 6)inverse 对称正定矩阵求逆函数 7)Calculate_BQBT 权逆阵传播函数 8)Calculate_q 权倒数计算函数 9)构造函数 10) 析构函数 11) GetStationNumber 点名获取函数 12) Inputdata 原始数据输入函数 13) Printdata 原始数据写至结果文件函数 14) ca_H0 高程近似值计算函数 15) ca_ATPA 组成法方程函数 16) ca_dX 高程平差值计算函数 17) ca_V 残差计算函数 18) PrintResult 平差值输出函数 19) Quasi_Stable 拟稳平差计算函数 20) main 主函数 上述函数全部为在 C 语言中自定义的函数。 4.3 函数功能说明及其源代码函数功能说明及其源代码 1.MyBreak 信息提示函数： 应用程序都需要通过计算机界面动态地显示一些提示性信息，例如，正定 矩阵求逆计算遇到被求逆的矩阵是非正定矩阵，程序就无法继续运行，应该通 过用户界面显示“矩阵求逆失败”的提示信息，然后终止程序运行。 函数源代码如下： void MyBreak(char* fmt, .) char buffer256; va_list argptr; va_start(argptr, fmt); vsprintf(buffer, fmt, argptr); va_end(argptr); #ifdef VC_EXTRALEAN AfxMessageBox(buffer); #else printf(buffer); getchar(); #endif / VC_EXTRALEAN 2.ij 对称矩阵下标计算函数 平差程序经常会处理对称矩阵，例如观测值的权矩阵、权逆阵、法方程系 数阵等都是对称矩阵，为了节省存储空间和避免重复计算，采用仅存下三角矩 阵（含主对角线元素）的存储方案存储对称矩阵，亦即将对称矩阵的下三角元 素按顺序存到数组中。对称矩阵的下标计算，就是当对称矩阵仅存下三角矩阵 时根据矩阵元素的行号和列号确定相应矩阵元素在数组中的存储位置。 函数源代码如下： int ij(int i,int j) return (i=j)? i*(i+1)/2+j :j*(j+1)/2+i; 3.PrintM 数组输出函数 设有双精度型数组 A，数组长度为 size，函数 PrintM 将数组 A 的元素写入 指定文件，函数原型如下： void PrintM(FILE *fp,double A,int size, int t,char* fmt, char* title,bool IsLabel) fp文件指针 A待输出的 double（双精度）型数组 Size数组的长度，即输出元素的总个数 t每行元素的个数，即每 t 个数据进行一次换行 fmt输出格式 title标题字符串地址，默认值为空 IsLabel当值为逻辑真时每行添加行号，为假时不加行号，默认值为真 函数源代码如下： void PrintM(FILE *fp,double A,int size, int t,char* fmt, char* title,bool IsLabel) if(title)fprintf(fp,n %s: ,title); int j=0; for(int i=0;isize;i+) if(i%t=0) j+; if(IsLabel)fprintf(fp,n%3d ,j); else fprintf(fp,n ); fprintf(fp,fmt,Ai); fprintf(fp,n); 4. PrintM2 对称矩阵输出函数 设有阶对称矩阵，程序中仅存矩阵下三角元素，矩阵元素为双精度型，nn 函数 PrintM2 将矩阵输出至指定文件，函数原型如下： void PrintM2(FILE* fp, double M, int n, int t,char *fmt, char* title,bool IsLabel) fp文件指针 M待输出的数组 n矩阵的阶数 t格式控制变量 fmt输出格式 title标题字符串地址，缺省值空 IsLabel值为 true 时每行前添加行号，等于 false 时不加行号，默认值为 true。 函数源代码如下： void PrintM2(FILE* fp, double M, int n, int t,char *fmt, char* title,bool IsLabel) if(title)fprintf(fp,n %s: ,title); int index=0; for(int i=0;in;i+) if(IsLabel)fprintf(fp,n%3d ,i+1); else fprintf(fp,n ); for(int j=0;j0 fprintf(fp,fmt,Mindex+); fprintf(fp,n); 5.线性方程组输出函数 设有线性方程组，A 是矩阵，b 是 n 维向量。函数将此方程0bAXtn 组的系数矩阵和常数项向量输出至指定文件，函数原型如下： void PrintEquation(FILE* fp, double A, double b, int n, int t, char *fmt, char* title) fp文件指针 A方程组系数矩阵，数组长度为tn b方程组常数项向量，数组长度为 n; n方程的个数 t方程未知数个数 fmt输出格式 title标题字符串地址，默认值为空字符串 函数源代码如下： void PrintEquation(FILE* fp, double A, double b, int n, int t, char *fmt, char* title) if(title)fprintf(fp,n %s: ,title); for(int i=0;in;i+) fprintf(fp,n%3d ,i+1); for(int j=0;jt;j+) fprintf(fp,fmt,Ai*t+j); fprintf(fp,fmt,bi); 6.inverse 对称正定矩阵求逆函数 从误差方程到权矩阵，到最后的求其平差值和改正数，都离不开矩阵的求 逆运算，函数原型： bool inverse(double a,int n) a函数调用前为待求逆矩阵的元素； n矩阵的阶数； 函数返回值若计算成功返回 true；若计算失败返回 false。 函数源代码如下： bool inverse(double a,int n) double *a0=new doublen; for(int k=0;kn;k+) double a00=a0; if(a00+1.0=1.0) delete a0; return false; for(int i=1;in;i+) double ai0 = ai*(i+1)/2; if(i=n-k-1)a0i= -ai0/a00; else a0i= ai0/a00; for(int j=1;j=i;j+) a(i-1)*i/2+j-1=ai*(i+1)/2+j+ai0*a0j; for(i=1;in;i+) a(n-1)*n/2+i-1=a0i; an*(n+1)/2-1=1.0/a00; delete a0; return true; 7. Calculate_BQBT 权逆阵传播函数 此函数的功能是类似协方差的传播定律，用于计算矩阵乘积计算，如 ， T BQB 计算公式 （4.3-1） T rn nn nrrr B Q BN 则 N 的第 i 行第 j 列元素的计算公式为 （4.3-2） 1 0 1 0 n k n x jskskiij bqbn 函数原型如下； void Calculate_BQBT(double B,double Q,int r,int n,double N)； B（4.3-1）中的矩阵，数组长度为；nr Q（4.3-1）中的矩阵，数组长度为；2/ ) 1( nn N（4.3-1）中的矩阵，数组长度为；2/ ) 1( rr rB 矩阵的行数； nB 矩阵的列数。 函数源代码如下： void Calculate_BQBT(double B,double Q,int r,int n,double N) for(int i=0;ir;i+) for(int j=0;j=i;j+) double nij=0.0; for(int k=0;kn;k+) for(int s=0;sn;s+) nij+=Bi*n+k*Qij(k,s)*Bj*n+s; Nij(i,j)+=nij; 8. Calculate_q 权倒数计算函数 计算说明：设 X 是 t 维随见向量，其权逆阵为阶对称矩阵 Q，z 是 Xtt 的函数： (4.3-3)FXz 其中 F 为系数向量，则 z 的权倒数为t1 (4.3-4) T z z FQFq p 1 由矩阵乘法公式，得 （4.3-5） 1 0 1 0 t k t s skskz fqfq 式中：是矩阵 Q 的第 k 行第 s 列元素。 ks q 函数原型： double Calculate_q(double Q,double F,int Fin,int n) Q权逆阵数组，仅存下三角矩阵，数组长度等于；2/ ) 1( tt F系数数组，数组长度为 t； t权逆阵的阶数； 返回值返回式（4.3-5）计算的结果。 函数源代码如下： double Calculate_q(double Q,double F,int Fin,int n) double q=0.0; for(int k=0;k0) delete Height; delete ATPA; delete ATPL; delete dX; for(int i=0; im_Pnumber;i+) if(Pnamei!=NULL)delete(Pnamei); delete Pname; 10. GetStationNumber 点名获取函数 工作流程如下： （1）将待查点名 name 与 Pname 数组中已经保存的点名利用系统内置函数 strcmp 逐一比较，检查 Pname 中是否已经保存有 name 这个点名，如果 Pname 中存在 name 这个点名，就返回到 name 在 Pname 中的下标。 （2）如果 Pname 中没有 name 这个点名，则向计算机内申请内存，然后利 用系统内置函数 strcpy 将 name 复制到申请的内存中，再将内存地址存到 Pname 数组中，再将内存地址存到 Pname 数组中，最后返回 name 在 Pname 中的下标. 函数源代码如下： int CLevelingAdjust: GetStationNumber(char *name) for(int i=0; im_Pnumber; i+) if(Pnamei!=NULL) /将待查点名与已经存入点名数组的点名比较 if(strcmp(name,Pnamei)=0)return i; else /待查点名是一个新的点名，将新点名的地址放到 Pname 数组 中 int len = strlen(name); Pnamei = new charlen+1; strcpy(Pnamei, name); return i; return -1; /Pname 数组已经存满，且没有待查点名 11. Inputdata 原始数据输入函数 该函数有一个形参 datafile,定义为 string 型地址变量，其内容是数据文件名 称字符串的地址，文件名包括文件的全路径及扩张名，如文件在程序所在的文 件目录下的路径为“data.txt” 。 为了从数据文件读取数据，函数内首先声明了 FILE(文件指针定义的关键 字)型变量 fp，并调用 fopen 函数打开数据文件，将文件的指针赋给 fp。接着， 按照数据文件中的数据内容的顺序与格式，先读取网的概况数据，在读取网的 已知高程数据，最后读取观测高差数据。 函数院代码如下： void CLevelingAdjust: Inputdata(char *datafile) FILE *fp; if(fp=fopen(datafile,r)=NULL) MyBreak(n 数据文件打不开!); exit(0); fscanf(fp,%d%d%d, int unPnumber=m_Pnumber-m_knPnumber; Height=new double m_Pnumber; dX=new double m_Pnumber; ATPA=new double m_Pnumber*(m_Pnumber+1)/2; ATPL=new double m_Pnumber; StartP=new int m_Lnumber; EndP=new int m_Lnumber; L=new double m_Lnumber; V=new double m_Lnumber; P=new double m_Lnumber; fscanf(fp,%lf, Pname=new char* m_Pnumber; for(int i=0;im_Pnumber;i+) / GetStationNumber 函数根据 Pnamei是否为 NULL / 确定 Pnamei是否为点名地址 Pnamei = NULL; char buffer100; /临时数组，保存从文件中读到的点名 / 读取已知高程数据 for( i=0;i=m_knPnumber-1;i+) fscanf(fp,%s,buffer); int c=GetStationNumber(buffer); fscanf(fp,%lf, / 读取观测数据 for(i=0;im_Lnumber;i+) fscanf(fp,%s,buffer); /读取高程起点名 StartPi=GetStationNumber(buffer); if(StartPi0) fprintf(resultfp,n 数据文件出错：); fprintf(resultfp,n 第%d 个高差的起始点名为%s,i+1,buffer); fclose(resultfp); exit(0); fscanf(fp,%s,buffer);/读取高程终点 EndPi=GetStationNumber(buffer); if(EndPi0) fprintf(resultfp,n 数据文件出错：); fprintf(resultfp,n 第%d 个高差终点的点名为%s,i+1,buffer); fclose(resultfp); exit(0); fscanf(fp,%lf%lf, /读取高差值与路线长度 Pi=1.0/Pi; fclose(fp); 12. Printdata 原始数据写至结果文件函数 此函数的功能就是将网的已知信息反馈到结果文件中。在结果文件中显示 网额信息，例如；网的总点数，总观测值，已知点个数。 函数源代码如下： void CLevelingAdjust: Printdata() int i; fprintf(resultfp,n 观测值总数: %d 总点数: %d 已知点数：%d n, m_Lnumber, m_Pnumber,m_knPnumber); fprintf(resultfp,n 验前单位权中误差：%lf ,m_Sigma); fprintf(resultfp,nn 已知高程:n); for(i=0;i=m_knPnumber-1;i+) fprintf(resultfp,n%5d %8s ,i+1,Pnamei); fprintf(resultfp,%10.4lf ,Heighti); fprintf(resultfp,nn 高差观测值:n); for(i=0;i=m_Lnumber-1;i+) fprintf(resultfp,n%5d%8s%8s,i+1,PnameStartPi,PnameEndPi); fprintf(resultfp,%12.4lf %10.3lf,Li,1.0/Pi); 13. ca_H0 近似高程计算 近似高程计算的思路 第一步，计算位置点标志。各点的高程值用 Height 数组保存，近似高程计 算计算之前，先将 Height 数组中的未知点的高程值赋值为-9999.9，由于正常的 高程值不可能小于-9999.9，根据高程数组中的值是否小于-9999.9 可判断某点是 否需要计算近似高程。当某点的近似高程计算出来之后，数组中-9999.9 就会被 实际的值所取代，直到高程数组中所有的高程值均大于-9999.9，近似高程计算 就此结束。 第二步，计算高程。近似高程计算的基本思路是，遍历观测值，找到每一 个观测值起始点号、终点号，再根据点号从 Height 数组中获得每一段高程起始 点的高程值和终点的高程值，当高差的一端是高程已知点、另一端是高程已知 点时，就由已知点高程和观测高差传算出未知点的高程值，并将计算结果赋给 高程数组，如果差两端同为已知点或者同为未知点，就跳到下一次循环，继续 查找下一个观测值。 函数源代码如下： void CLevelingAdjust: ca_H0() for(int i=m_knPnumber;i(m_Pnumber-m_knPnumber) fprintf(resultfp,nn 下列点无法计算出概略高程:n); for(i=0;im_Pnumber;i+) if(Heighti-9999.0) printf(n%s,Pnamei); MyBreak(近似高程计算失败！); fclose(resultfp); exit(0); 14. ca_ATPA 组成法方程函数 该函数的作用是组成法方程和，计算流程如下：PAATPlAT 1)获得高差的起点号 i 和终点号 j 2)获得起点和终点的高程值； ji HH , 3)计算误差方程的自由项； kjik hHHl 4)将累加到法方程系数阵和法方程自由项的相应单元中。 kkk lpp 、 已经在 InputData 函数中根据路线长度计算出来，保存在权数组 P 中。 k p 函数源代码如下： void CLevelingAdjust: ca_ATPA() int t=m_Pnumber; for(int i=0; it*(t+1)/2; i+) ATPAi=0.0; for(i=0; it; i+) ATPLi=0.0; for(int k=0; km_Lnumber; k+) int i=StartPk; int j=EndPk; double Pk=Pk; double Lk=Lk-(Heightj-Heighti); ATPLi-=Pk*Lk; ATPLj+=Pk*Lk; ATPAij(i,i)+=Pk; ATPAij(j,j)+=Pk; ATPAij(i,j)-=Pk; 15. ca_dX 高程平差值计算函数 该函数首先调用对称正定矩阵求逆函数 inverse，计算法方程系数矩阵的逆 阵，然后以总点数为循环界逐点计算各点高程的改正数，将高程改正数保存 dX 数组中，同时将高程改正数累加到数组 Height 的对应单元中。这样，Height 数 组就从高程近似值变成了高程平差值。 函数源代码如下： void CLevelingAdjust: ca_dX() if(!inverse(ATPA,m_Pnumber) MyBreak(n 法方程系数矩阵降秩!); fclose(resultfp); exit(0); for(int i=0; im_Pnumber; i+) double xi=0.0; for(int j=0; jm_Pnumber; j+) xi+=ATPAij(i,j)*ATPLj; dXi=xi; Heighti+=xi; 16. ca_V 残差计算函数 残差也叫观测值的平差改正数，在参数平差中观测值的平差值都可以直接 用参数平差值计算出来，所以计算残差的目的不是用来改正观测值，而主要是 用残差来进行精度估计。残差计算由函数 ca_V 完成，计算结果存在数组 V 中， 该函数的返回值是pvv。 以观测值序号 k 为循环变量，按观测值循环，按下式计算各观测值的残差 ， k v （k=0,1,2,n-1） kjik hxxv 式中 k 为观测值编号； 是观测高差； k h 是观测值的平差改正数，是称为残差； k v 表示高差两端点的编号（即点号） ；ji, 分别表示观测高差起点和终点的高程平差值，即平差中未知数。 ji xx、 第 k 次循环中所要进行的工作包括： 1)获得高差的起点号 i 和终点 j； 2)获得起点和终点的高程值（已经是平差值） ； ji HH , 3)计算残差，存在 Vk中； kijk hHHv 4)将累加到pvv中。 2 kkv p 结束循环之后，得到pvv，返回pvv. 函数的源代码如下： double CLevelingAdjust: ca_V() double pvv=0.0; for(int i=0;i=m_Lnumber-1;i+) int k1=StartPi; int k2=EndPi; Vi=Heightk2-Heightk1-Li; pvv+=Vi*Vi*Pi; return(pvv); 17. PrintResult 精度估计与平差结果输出 该函数将平差结果以 txt 的文件格式写至结果文件。输出的内容包括高程点 名、高程平差值及其中误差，高差端点点名、高差平差值、高差改正数、高差 平差值的中误差。 函数的源代码如下： void CLevelingAdjust: PrintResult() fprintf(resultfp,nn = 高程平差值及其精度 =n); fprintf(resultfp,n 点名 近似高程 改正数 高程平差值 中误 差n); for(int i=0; im_Pnumber; i+) fprintf(resultfp,n %5s ,Pnamei); double dx=dXi; double qii=ATPAij(i,i); fprintf(resultfp,%12.4lf%9.4lf%12.4lf%9.4lf, Heighti-dx,dx,Heighti,sqrt(qii)*m_mu); fprintf(resultfp,nnn = 观测值平差值及其精度 =n); fprintf(resultfp,n No. 起 点 终 点 观测高差 ); fprintf(resultfp, 高差平差值 观测权 中误差n); for(i=0;im_Pnumber | m_StablePnumber1) MyBreak(拟稳点数错误!); fclose(resultfp); exit(0); IsStable=new intm_Pnumber; for(int i=0;i=m_Pnumber-1;i+) IsStablei=false; fprintf(resultfp,n%sn,n 拟稳点：); for( i=0;i=m_StablePnumber-1;i+) char name100; fscanf(fp,%s,name); int k=GetStationNumber(name); if(k0) fprintf(resultfp,n 拟稳点名文件出错：); fpr