降沉监测网拟稳平差及程序设计.doc

本科毕业论文 沉降监测网拟稳平差及程序设计SUBSIDENCE MONITORING NETWORK QUASI-STABLE ADJUSTMENT AND PROGRAM DESIGN学院（部）： 测绘学院 专业班级： 测绘 11-1 学生姓名： 指导教师： 年 月 日沉降监测网拟稳平差及程序设计摘要本文简单的介绍了沉降监测网的起因、经过、形成，以及对沉降观测网数据处理的方法。拟稳平差在沉降观测网的数据处理应用，使用C语言对拟稳平差进行编程，以函数的形式实现拟稳平差中的误差方程、组建法方程以及矩阵求逆等功能。本文还介绍了有关拟稳点的判断和选取的内容，最后详细的说明了运行程序的一些必要步骤和上机操作问题。关键词：沉降监测网，C语言，拟稳平差，程序设计SUBSIDENCE MONITORING NETWORK QUASI-STABLE ADJUSTMENT AND PROGRAM DESIGNABSTRACTThis 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目录ABSTRACT31沉降监测网概述61.1 沉降监测网的形成62地下开采引起的地表移动和破坏62.1地表移动的形式62.2地表移动盆地73沉降监测网的数据来源73.1开采沉陷的观测工作73.2本次实测数据说明84拟稳平差程序概述94.1拟稳平差的基本原理94.2拟稳平差程序编译的相关函数104.3函数功能说明及其源代码114.4原始数据文件格式设计304.5拟稳点信息的读取与存储314.6拟稳点的判定与选取314.6.1拟稳点的判定方法314.6.2拟稳点的选择原则325.拟稳平差程序的运行335.1原始数据的处理335.2程序运行的步骤345.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）式中：为n维自由项向量；V是n维残差向量；A为系数矩阵；X是t维参数向量，这里X是高程点近似值的平差改正数。由上式（4-1）组成法方程式；（式4-2）式中，P为观测值的权矩阵，对称正定。水准网平差中至少有一个已知高程点。假如网中没有已知高程点，误差方程式的系数阵A的秩就会小于t，法方程的系数阵也会降秩。拟稳平差将全部高程点分为拟稳点和非拟稳点两部分，仅在拟稳点参数的平方和最小的约束下求法方程式的解。设水准网中共有t个点，其中s个点为拟稳点，为拟稳点参数解向量，约束条件为（式4-3）我们采用虚拟观测值法求解。设G、G2都是t维向量， 其中（i=0,1,t-1）为 i号点位拟稳点（式4-5） i号点为非拟稳点 拟稳平差参数解及其权逆阵的公式为（式4-6）（式4-7）单位权中误差计算公式为（式4-8）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_EXTRALEANAfxMessageBox(buffer);#elseprintf(buffer);getchar();#endif / VC_EXTRALEAN2.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对称矩阵输出函数设有阶对称矩阵，程序中仅存矩阵下三角元素，矩阵元素为双精度型，函数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 & j%t=0)fprintf(fp,n );fprintf(fp,fmt,Mindex+);fprintf(fp,n);5.线性方程组输出函数设有线性方程组，A是矩阵，b是n维向量。函数将此方程组的系数矩阵和常数项向量输出至指定文件，函数原型如下：void PrintEquation(FILE* fp, double A, double b,int n, int t, char *fmt, char* title)fp文件指针A方程组系数矩阵，数组长度为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权逆阵传播函数此函数的功能是类似协方差的传播定律，用于计算矩阵乘积计算，如，计算公式（4.3-1）则N的第i行第j列元素的计算公式为（4.3-2）函数原型如下；void Calculate_BQBT(double B,double Q,int r,int n,double N)；B（4.3-1）中的矩阵，数组长度为；Q（4.3-1）中的矩阵，数组长度为；N（4.3-1）中的矩阵，数组长度为；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是X的函数：(4.3-3)其中F为系数向量，则z的权倒数为(4.3-4)由矩阵乘法公式，得（4.3-5）式中：是矩阵Q的第k行第s列元素。函数原型：double Calculate_q(double Q,double F,int Fin,int n)Q权逆阵数组，仅存下三角矩阵，数组长度等于；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;kn;k+)int i=Fink;for(int s=0;s0)delete StartP; delete EndP; delete L; delete P; delete V; if(m_Pnumber0)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,&m_Lnumber,&m_Pnumber,&m_knPnumber); 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,&m_Sigma); 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,&Heightc); / 读取观测数据 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,&Li,&Pi); /读取高差值与路线长度 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就会被实际的值所取代，直到高程数组中所有的