基于c#CP3平面网严密平差数据处理

基于c#CP3平面网严密平差数据处理 using System; using System.Collections.Generic; using System.Collections;//使用动态数组需要添加的语句 using System.ComponentModel; using System.Data; using System.Drawing; using System.Linq; using System.Text; using System.Windows.Forms; using System.IO; using ParaSet; using CSAccelerateMatrix; using Trans; using CSUsualFun; namespace FreeStaAdj { public partial class Form1 : Form { static string Path = "C:\\Documents and Settings\\HB\\桌面\\春儿数据-120907"; string[] KPName;//存储已知点点号 double[,] KPXYZ;//存储已知点三维坐标 string[] StaKName;//存储测站点点号 double[,] StaKXYZ;//存储测站点已知三维坐标 public class StaData { public string StaName;//测站名 public int Num;//观测的CPIII点数 public string[] CP3Name = new string[12];//CPIII自由测站学习点数一般不会超过12个,此处也可用动态数组存储 public double[,] LVS = new double[12, 3]; }; ArrayList StaList = new ArrayList(); public struct DLVS//用于存储方差分量估计各次的方差估值 { public double DL,DV,DS; }; //以下为默认先验平差参数 double PriLm=0.5,PriSa=1.0,PriSb=1.0;//水平方向中误差，测距固定误差和比例误差 int DMaxNum = 10;//最大方差分量估计次数 double DLimitValue = 0.01;//方差分量估计默认收敛阀值 bool IsD=true;//true表示进行方差分量估计 bool SP = true;//true表示A+B*S bool Ma = true;//true表示验后单位权中误差 double LVRatio = 2.0;//水平方向 和 天顶距 权的比值 string AdjMethord = "普通";//三维平差方法 const double Pe = 0.0174532925199433;//Math.PI/180.0,便于后面简化运算 public Form1() { InitializeComponent(); } public double Cal_DRatio(double a,double b,double c) { double x, y; if (a>=b) { x=a; y =b; } else { x = b; y=a; } if (xc) { y=c; } return x / y;//三个数中的最大值:最小值 } public double atan(double y, double x)//Math.Atan2的返回值区间(-PI,PI] { double p = Math.Atan2(y, x); if (p 0) { linecount++; } } A0.Close(); StaKName = new string[linecount]; StaKXYZ = new double[linecount, 3]; int k = 0; StreamReader A1 = new StreamReader(fname); while ((line = A1.ReadLine())!= null) { //原文件格式是YXZ GetLVS(line, ref StaKName[k], ref StaKXYZ[k, 1], ref StaKXYZ[k, 0], ref StaKXYZ[k, 2]); k++; } A1.Close(); } } private void 参数设置ToolStripMenuItem_Click(object sender, EventArgs e) { SetPa f1=new SetPa(); f1.ShowDialog(this); this.PriLm = double.Parse(f1.textBox1.Text); this.PriSa = double.Parse(f1.textBox2.Text); this.PriSb = double.Parse(f1.textBox3.Text); this.DMaxNum = int.Parse(f1.textBox4.Text); this.DLimitValue = double.Parse(f1.textBox5.Text); this.IsD =f1.checkBox1.Checked; this.SP = f1.radioButton1.Checked; this.Ma = f1.radioButton4.Checked; this.AdjMethord = f1.comboBox1.Text; this.LVRatio = double.Parse(f1.comboBox2.Text); } private void 三维平差ToolStripMenuItem_Click(object sender, EventArgs e) { //以下进行测站三维平差-只针对CPIII自由测站的三维解算 double[,] StaAppXYZ=new double[StaList.Count,3];//测站概略坐标数组 double[,] StaAdjXYZ = new double[StaList.Count, 3];//测站概略坐标数组 double[] StaZ = new double[StaList.Count];//测站定向角未知数数组 double[] StaAppZ = new double[StaList.Count];//测站定向角未知数数组 double[] K = new double[StaList.Count];//大气折光系数k数组 double[,] AB = new double[StaList.Count, 2];//固定误差 和 比例误差数组 CalStaAppXY(StaAppXYZ,StaList,KPName,KPXYZ);//计算测站点近似坐标 CopyStaAppXY(StaAppXYZ, StaAdjXYZ);//将近似坐标拷贝一个副本 double z0=0;//z0为某测站定向角未知数近似值,单位为弧度 ArrayList D0List = new ArrayList();//存储各次方差分量估值的数组 //平差方法包括:普通 球气差 距离误差 球气差&距离误差 四种 switch (AdjMethord) { case "普通": { //未知参数改正数顺序为 定向角z,X,Y,Z int i; for (i = 0; i < StaList.Count; i++) { double[] NL = new double[10];//水平方向 法方程系数矩阵 double[] NV = new double[10];//天顶距 法方程系数矩阵 double[] NS = new double[10];//边长 法方程系数矩阵 double[] N = new double[10];//法方程系数矩阵 double[] WL = new double[4];//水平方向 法方程常数项 double[] WV = new double[4];//天顶距 法方程常数项 double[] WS = new double[4];//天顶距 法方程常数项 double[] W = new double[4];//法方程常数项 double[] X = new double[4];//未知参数改正数顺序为z,X,Y,Z StaData D = (StaData)StaList[i];//获取各测站的数据,需要强制类型转换获取StaList中的数据 //以下对边长定权 double[] PS = new double[D.Num]; Cal_SWeight(PS, PriLm, PriSa, PriSb, D, SP); double[,] BL=new double[D.Num,2];//用于v=Bx-l计算v double[,] BV=new double[D.Num,3]; double[,] BS=new double[D.Num,3]; double[] lL=new double[D.Num];//用于v=Bx-l计算l double[] lV=new double[D.Num]; double[] lS=new double[D.Num]; double[] LV=new double[D.Num];//用于存储残差 double[] VV=new double[D.Num]; double[] SV=new double[D.Num]; int IterNum = 0;//记录迭代次数 double DL, DV, DS,RLS; ArrayList TempD0 = new ArrayList();//存储各测站方差分量估计各次各分量的方差值 do { IterNum++; if (IterNum > DMaxNum) { break;//超过自定义迭代次数则跳出循环 } //以下填充N和W,此处的三个函数实际上合并成一个函数可以简化运算！但时间限制,有待继续完善 FillNS(NS,WS,i,D,PS,KPName,KPXYZ,StaAdjXYZ,BS,lS); FillNL(NL, WL, i, D, ref z0, KPName, KPXYZ, StaAdjXYZ, BL, lL); FillNV(NV, WV, i, D, LVRatio, KPName, KPXYZ, StaAdjXYZ, BV, lV); AccelerateMatrix A = new AccelerateMatrix(); A.AddN(N, NL, NS, 10); A.AddN(N, NV, 10); A.AddN(W, WL, WS, 4); A.AddN(W, WV, 4); A.SPDMatrixInverse(ref N, 4); A.Cal_InvNW(N, W, 4, ref X); Cal_LV(LV,BL,lL,X,D.Num); Cal_SV(SV,BS,lS,X,D.Num); Cal_VV(VV,BV,lV,X,D.Num); UsualFun UF = new UsualFun(); StaAppZ[i] = UF.ChangeRadToDMS(z0);//存储定向角未知数近似值 StaZ[i] = UF.ChangeRadToDMS(z0 + X[0] / 206264.8062471);//将定向角从rad转换为dms StaAdjXYZ[i, 0] += X[1] / 1000; StaAdjXYZ[i, 1] += X[2] / 1000; StaAdjXYZ[i, 2] += X[3] / 1000; if(!IsD) { break;//不进行Helmert方差分量估计就跳出循环 } DL=Cal_LD0(LV,D.Num,N,NL,4); DV=Cal_VD0(VV,LVRatio,D.Num,N,NV,4); DS=Cal_SD0(SV,PS,D.Num,N,NS,4); DLVS Di=new DLVS(); Di.DL = DL; Di.DV = DV; Di.DS = DS; TempD0.Add(Di); RLS=DL/ DS; for (int j = 0; j < D.Num;j++) { PS[j] *=RLS; } LVRatio*=DV/DL;//此处不理解可以推导 } while ((Math.Abs(Cal_DRatio(DL, DV, DS)- 1)) > DLimitValue); D0List.Add(TempD0); } PrintResult(StaAppXYZ, StaAdjXYZ, StaZ,StaAppZ,StaList, D0List, IsD, StaKName, StaKXYZ); break; } case "球气差": { //未知参数改正数顺序为 定向角z,k,X,Y,Z int i; for (i = 0; i < StaList.Count; i++) { double[] NL = new double[15];//水平方向 法方程系数矩阵 double[] NV = new double[15];//天顶距 法方程系数矩阵 double[] NS = new double[15];//边长 法方程系数矩阵 double[] N = new double[15];//法方程系数矩阵 double[] WL = new double[5];//水平方向 法方程常数项 double[] WV = new double[5];//天顶距 法方程常数项 double[] WS = new double[5];//天顶距 法方程常数项 double[] W = new double[5];//法方程常数项 double[] X = new double[5];//未知参数改正数顺序为z,X,Y,Z StaData D = (StaData)StaList[i];//获取各测站的数据,需要强制类型转换获取StaList中的数据 //以下对边长定权 double[] PS = new double[D.Num]; Cal_SWeight(PS, PriLm, PriSa, PriSb, D, SP); double[,] BL = new double[D.Num, 2];//用于v=Bx-l计算v double[,] BV = new double[D.Num, 4];//k,X,Y,Z double[,] BS = new double[D.Num, 3]; double[] lL = new double[D.Num];//用于v=Bx-l计算l double[] lV = new double[D.Num]; double[] lS = new double[D.Num]; double[] LV = new double[D.Num];//用于存储残差 double[] VV = new double[D.Num]; double[] SV = new double[D.Num]; int IterNum = 0;//记录迭代次数 double DL, DV, DS, RLS; ArrayList TempD0 = new ArrayList();//存储各测站方差分量估计各次各分量的方差值 do { IterNum++; if (IterNum > DMaxNum) {