基于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 (x<c)
{
x=c;
}
if (y>c)
{
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)
{
return p+2*Math.PI;
}
else
{
return p;
}
}
public double SubRad(double y, double x)//返回两个弧度(角度)的差值
{
double p =y-x;
if (p < 0)
{
return p + 6.28318530717959;
}
else
{
return p;
}
}
public void GetXYZ(string PName,ref double x,ref double y,ref double z,string[] KPName, double[,] KPXYZ)
{
for (int i=0;i<KPName.Length;i++)
{
if (PName==KPName[i])
{
x=KPXYZ[i,0];
y=KPXYZ[i,1];
z=KPXYZ[i,2];
break;
}
}
}
public double GetlL(double Ljk,double z0,double ajk)
{
double ljk =Ljk+z0-ajk;
if (ljk>3.14)//出现359 59 59的情况
{
ljk -= 2 * Math.PI;
}
return ljk * 206264.8062471;
}
public void Cal_LVk(double[] LV, double[,] BL, double[] lL, double[] X, int n)
{
for (int i = 0; i < n; i++)
{
LV[i] = -1.0 + BL[i, 0] * X[2] + BL[i, 1] * X[3] - lL[i];
}
}
public void Cal_LV(double[] LV,double[,] BL, double[] lL,double[] X,int n)
{
for(int i=0;i<n;i++)
{
LV[i] = -1.0 + BL[i, 0] * X[1]+BL[i,1] * X[2]-lL[i];
}
}
public void Cal_SVk(double[] SV, double[,] BS, double[] lS, double[] X, int n)
{
for (int i = 0; i < n; i++)
{
SV[i] = BS[i, 0] * X[2] + BS[i, 1] * X[3] + BS[i, 2] * X[4] - lS[i];
}
}
public void Cal_SVab(double[] SV, double[,] BS, double[] lS, double[] X, int n)
{
for (int i = 0; i < n; i++)
{
SV[i] = -( X[4]+BS[i, 3] * X[5])+BS[i, 0] * X[1] + BS[i, 1] * X[2] + BS[i, 2] * X[3] - lS[i];
}
}
public void Cal_SV(double[] SV, double[,] BS, double[] lS, double[] X, int n)
{
for (int i = 0; i < n; i++)
{
SV[i] = BS[i, 0] * X[1] + BS[i, 1] * X[2] + BS[i, 2] * X[3] - lS[i];
}
}
public void Cal_VVk(double[] VV, double[,] BV, double[] lV, double[] X, int n)
{
for (int i = 0; i < n; i++)
{
VV[i] = BV[i, 0] * X[1]+BV[i, 1] * X[2] + BV[i, 2] * X[3] + BV[i, 3] * X[4] - lV[i];
}
}
public void Cal_VVabk(double[] VV, double[,] BV, double[] lV, double[] X, int n)
{
for (int i = 0; i < n; i++)
{
VV[i] = BV[i, 0] * X[1] + BV[i, 1] * X[2] + BV[i, 2] * X[3] + BV[i, 3] * X[6] - lV[i];
}
}
public void Cal_VV(double[] VV, double[,] BV, double[] lV, double[] X, int n)
{
for (int i = 0; i < n; i++)
{
VV[i] = BV[i, 0] * X[1] + BV[i, 1] * X[2] + BV[i, 2] * X[3] - lV[i];
}
}
public double Cal_LD0(double[] LV,int n,double[] InvN,double[] NL,int r)
{
double v2 = 0;
for (int i = 0; i < n;i++ )
{
v2 +=LV[i]*LV[i];
}
AccelerateMatrix A = new AccelerateMatrix();
return v2/(n-A.Tr_TwoPSMatrix(InvN, NL, r));
}
public double Cal_VD0(double[] VV, double LVRatio, int n, double[] InvN, double[] NV, int r)
{
double v2 = 0;
for (int i = 0; i < n; i++)
{
v2 += VV[i] * VV[i] / LVRatio;
}
AccelerateMatrix A = new AccelerateMatrix();
return v2 / (n - A.Tr_TwoPSMatrix(InvN, NV, r));
}
public double Cal_SD0(double[] SV, double[] PS, int n, double[] InvN, double[] NS, int r)
{
double v2 = 0;
for (int i = 0; i < n; i++)
{
v2 += SV[i] * SV[i]*PS[i];
}
AccelerateMatrix A = new AccelerateMatrix();
return v2 / (n - A.Tr_TwoPSMatrix(InvN, NS, r));
}
public void FillNVk(double[] NV, double[] WV, int Index, StaData D,//Index表示第几个测站
double LVRatio, string[] KPName, double[,] KPXYZ, double[,] StaAppXYZ,
double[,] BV, double[] lV)
{
//此函数考虑大气折光影响
//先计算天顶距的初始权,水平方向权的默认值为1
double PV = 1.0 / LVRatio;
//用天顶距观测值填充NV和WV矩阵
double Stax = StaAppXYZ[Index, 0];
double Stay = StaAppXYZ[Index, 1];
double Staz = StaAppXYZ[Index, 2];
double x = 0, y = 0, z = 0, dx, dy, dz, a, b, c,d, Sij, S2, lv;
for (int i = 0; i < D.Num; i++)
{
GetXYZ(D.CP3Name[i], ref x, ref y, ref z, KPName, KPXYZ);
dx = x - Stax;
dy = y - Stay;
dz = z - Staz;
Sij = Math.Pow(dx, 2) + Math.Pow(dy, 2);
S2 = Sij + Math.Pow(dz, 2);
Sij = Math.Sqrt(Sij);
lv = 206264.8062471 * (D.LVS[i, 1] * Pe - Math.Atan2(Sij, dz));//常数项
a = 206.2648062471 * dz / (Sij * S2);
b = -a * dy;
a *= -dx;
c = 206.2648062471 * Sij / S2;
d = 206.2648062471 * Sij / 6371;//大气折光系数k在误差方程式中的系数
BV[i, 0] = d;//未知数顺序k,X,Y,Z
BV[i, 1] = a;
BV[i, 2] = b;
BV[i, 3] = c;
lV[i] = lv;
NV[2] += PV * d * d;
NV[4] += PV * d * a;
NV[7] += PV * d * b;
NV[11] +=PV * d * c;
NV[5] += PV * a * a;
NV[8] += PV * a * b;
NV[9] += PV * b * b;
NV[12] += PV * a * c;
NV[13] += PV * b * c;
NV[14] += PV * c * c;
lv *= PV;
WV[1] += d * lv;
WV[2] += a * lv;
WV[3] += b * lv;
WV[4] += c * lv;
}
}
public void FillNVabk(double[] NV, double[] WV, int Index, StaData D,//Index表示第几个测站
double LVRatio, string[] KPName, double[,] KPXYZ, double[,] StaAppXYZ,
double[,] BV, double[] lV)
{
//此函数考虑大气折光影响,并估计距离的 固定误差 和 比例误差
//先计算天顶距的初始权
double PV = 1.0 / LVRatio;
//用天顶距观测值填充NV和WV矩阵
double Stax = StaAppXYZ[Index, 0];
double Stay = StaAppXYZ[Index, 1];
double Staz = StaAppXYZ[Index, 2];
double x = 0, y = 0, z = 0, dx, dy, dz, a, b, c,d, Sij, S2, lv;
for (int i = 0; i < D.Num; i++)
{
GetXYZ(D.CP3Name[i], ref x, ref y, ref z, KPName, KPXYZ);
dx = x - Stax;
dy = y - Stay;
dz = z - Staz;
Sij = Math.Pow(dx, 2) + Math.Pow(dy, 2);
S2 = Sij + Math.Pow(dz, 2);
Sij = Math.Sqrt(Sij);
lv = 206264.8062471 * (D.LVS[i, 1] * Pe - Math.Atan2(Sij, dz));//常数项
a = 206.2648062471 * dz / (Sij * S2);
b = -a * dy;
a *= -dx;
c = 206.2648062471 * Sij / S2;
d = 206.2648062471 * Sij / 6371;
BV[i, 0] = a;
BV[i, 1] = b;
BV[i, 2] = c;
BV[i, 3] = d;
lV[i] = lv;
NV[2] += PV * a * a;
NV[4] += PV * a * b;
NV[5] += PV * b * b;
NV[7] += PV * a * c;
NV[8] += PV * b * c;
NV[9] += PV * c * c;
NV[22] += PV * a * d;
NV[23] += PV * b * d;
NV[24] += PV * c * d;
NV[27] += PV * d * d;
lv *= PV;
WV[1] += a * lv;
WV[2] += b * lv;
WV[3] += c * lv;
WV[6] += d * lv;
}
}
public void FillNV(double[] NV, double[] WV, int Index, StaData D,//Index表示第几个测站
double LVRatio,string[] KPName, double[,] KPXYZ, double[,] StaAppXYZ,
double[,] BV, double[] lV)
{
//此函数不考虑大气折光影响
//先计算天顶距的初始权
double PV = 1.0 / LVRatio;
//用天顶距观测值填充NV和WV矩阵
double Stax = StaAppXYZ[Index, 0];
double Stay = StaAppXYZ[Index, 1];
double Staz = StaAppXYZ[Index, 2];
double x=0,y=0,z=0,dx,dy,dz,a,b,c, Sij,S2, lv;//dx,dy方向的坐标差,近似距离,水平方向常数项
for (int i = 0; i < D.Num; i++)
{
GetXYZ(D.CP3Name[i], ref x, ref y, ref z, KPName, KPXYZ);
dx = x - Stax;
dy = y - Stay;
dz = z - Staz;
Sij = Math.Pow(dx, 2) + Math.Pow(dy, 2);
S2 = Sij + Math.Pow(dz, 2);
Sij = Math.Sqrt(Sij);
lv = 206264.8062471*(D.LVS[i, 1] * Pe - Math.Atan2(Sij, dz));//常数项
a = 206.2648062471 * dz / (Sij * S2);
b = -a *dy;
a *= -dx;
c = 206.2648062471 * Sij / S2;
BV[i, 0] = a;
BV[i, 1] = b;
BV[i, 2] = c;
lV[i] =lv;
NV[2] += PV * a * a;
NV[4] += PV * a * b;
NV[5] += PV * b * b;
NV[7] += PV * a * c;
NV[8] += PV * b * c;
NV[9] += PV * c * c;
lv *= PV;
WV[1] += a * lv;
WV[2] += b * lv;
WV[3] += c * lv;
}
}
public void FillNLk(double[] NL, double[] WL, int Index, StaData D,//Index表示第几个测站
ref double z0, string[] KPName, double[,] KPXYZ, double[,] StaAppXYZ,
double[,] BL, double[] lL)
{
//z0定向角未知数
//用方向观测值填充NL和WL矩阵
double Stax = StaAppXYZ[Index, 0];//测站近似坐标
double Stay = StaAppXYZ[Index, 1];
double Staz = StaAppXYZ[Index, 2];
//先计算测站的定向角未知数
double x = 0, y = 0, z = 0;
z0 = 0;//初始化为0
for (int i = 0; i < D.Num; i++)
{
GetXYZ(D.CP3Name[i], ref x, ref y, ref z, KPName, KPXYZ);
z0 += SubRad(atan(y - Stay, x - Stax), D.LVS[i, 0] * Pe);
}
z0 /= D.Num;//求平均值,单位为弧度[0,2*PI]
double dx, dy, a, b, S02, ljk;//dx,dy方向的坐标差,近似距离,水平方向常数项
for (int i = 0; i < D.Num; i++)
{
GetXYZ(D.CP3Name[i], ref x, ref y, ref z, KPName, KPXYZ);
dx = x - Stax;
dy = y - Stay;
S02 = Math.Pow(dx, 2) + Math.Pow(dy, 2);
ljk = GetlL(D.LVS[i, 0] * Pe, z0, atan(dy, dx));
a = 206.2648062471 * dy / S02;
b = -206.2648062471 * dx / S02;
BL[i, 0] = a;
BL[i, 1] = b;
lL[i] = ljk;
//水平方向观测值的权为1.0
NL[0] += 1.0;
NL[3] -= a;
NL[5] += a * a;
NL[6] -= b;
NL[8] += a * b;
NL[9] += b * b;
WL[0] -= ljk;
WL[2] += a * ljk;
WL[3] += b * ljk;
}
}
public void FillNL(double[] NL, double[] WL, int Index, StaData D,//Index表示第几个测站
ref double z0,string[] KPName, double[,] KPXYZ, double[,] StaAppXYZ,
double[,] BL,double[] lL)
{
//z0定向角未知数
//用方向观测值填充NL和WL矩阵
double Stax = StaAppXYZ[Index, 0];//测站近似坐标
double Stay = StaAppXYZ[Index, 1];
double Staz = StaAppXYZ[Index, 2];
//先计算测站的定向角未知数
double x = 0, y = 0, z = 0;
z0 = 0;//初始化为0
for (int i = 0; i < D.Num; i++)
{
GetXYZ(D.CP3Name[i], ref x, ref y, ref z, KPName, KPXYZ);
z0+=SubRad(atan(y-Stay,x-Stax),D.LVS[i,0]*Pe);
}
z0 /= D.Num;//求平均值,单位为弧度[0,2*PI]
double dx,dy,a, b,S02, ljk;//dx,dy方向的坐标差,近似距离,水平方向常数项
for (int i = 0; i < D.Num; i++)
{
GetXYZ(D.CP3Name[i], ref x, ref y, ref z, KPName, KPXYZ);
dx =x-Stax;
dy =y-Stay;
S02 =Math.Pow(dx, 2) + Math.Pow(dy, 2);
ljk=GetlL(D.LVS[i, 0] * Pe,z0,atan(dy, dx));
a = 206.2648062471 * dy / S02;
b = -206.2648062471 * dx / S02;
BL[i, 0] = a;
BL[i, 1] = b;
lL[i] = ljk;
//水平方向观测值的权为1.0
NL[0] += 1.0;
NL[1] -= a;
NL[2] += a * a;
NL[3] -= b;
NL[4] += a * b;
NL[5] += b * b;
WL[0] -= ljk;
WL[1] += a * ljk;
WL[2] += b * ljk;
}
}
public void PrintResultk(
double[,] AppXYZ, double[,] AdjXYZ, double[] StaZ, double[] StaAppZ, double[] K,
ArrayList StaList, ArrayList D0List, bool IsD0,
string[] StaKName, double[,] StaKXYZ)
{
richTextBox1.Clear();
int i, j;
//先输出定向角未知数计算结果
richTextBox1.Text += "测站 定向角近似值(°′″) 定向角平差值(°′″) 改正数(″)\n";
for (i = 0; i < StaList.Count; i++)
{
StaData D = (StaData)StaList[i];//获取各测站的数据,需要强制类型转换获取StaList中的数据
richTextBox1.Text +=
string.Format("{0,4}", D.StaName) +
string.Format("{0,17:f4}", StaAppZ[i]) +
string.Format("{0,23:f4}", StaZ[i]) +
string.Format("{0,17:f2}", (StaZ[i] - StaAppZ[i]) * 10000) + "\n";
}
richTextBox1.Text += "\n-------------------------------------------------\n\n";
//输出大气折光系数估值
richTextBox1.Text += "测站 大气折光系数估值\n";
for (i = 0; i < StaList.Count; i++)
{
StaData D = (StaData)StaList[i];//获取各测站的数据,需要强制类型转换获取StaList中的数据
richTextBox1.Text +=
string.Format("{0,4}", D.StaName) +
string.Format("{0,14:f2}", K[i]) +"\n";
}
richTextBox1.Text += "\n-------------------------------------------------\n\n";
richTextBox1.Text += " 既有坐标(0) | 近似坐标(1) | 平差坐标(2) | 坐标较差[(1)-(0)] | 坐标较差[(2)-(0)]\n";
richTextBox1.Text += "点名 | X(m) Y(m) Z(m) | X(m) Y(m) Z(m) | X(m) Y(m) Z(m) | dx(mm) dy(mm) dz(mm) | dx(mm) dy(mm) dz(mm)\n";
for (i = 0; i < StaKName.Length; i++)
{
richTextBox1.Text +=
string.Format("{0,4}", StaKName[i]) + " | " +
string.Format("{0:0.0000}", StaKXYZ[i, 0]) + " " +
string.Format("{0:0.0000}", StaKXYZ[i, 1]) + " " +
string.Format("{0:0.0000}", StaKXYZ[i, 2]) + " | ";
for (j = 0; j < StaList.Count; j++)
{
StaData D = (StaData)StaList[j];
if (D.StaName == StaKName[i])
{
richTextBox1.Text +=
string.Format("{0:0.0000}", AppXYZ[j, 0]) + " " +
string.Format("{0:0.0000}", AppXYZ[j, 1]) + " " +
string.Format("{0:0.0000}", AppXYZ[j, 2]) + " | ";
richTextBox1.Text +=
string.Format("{0:0.0000}", AdjXYZ[j, 0]) + " " +
string.Format("{0:0.0000}", AdjXYZ[j, 1]) + " " +
string.Format("{0:0.0000}", AdjXYZ[j, 2]) + " |";
richTextBox1.Text +=
string.Format("{0,6:f1}", (AppXYZ[j, 0] - StaKXYZ[i, 0]) * 1000) +
string.Format("{0,7:f1}", (AppXYZ[j, 1] - StaKXYZ[i, 1]) * 1000) +
string.Format("{0,7:f1}", (AppXYZ[j, 2] - StaKXYZ[i, 2]) * 1000) + " |";
richTextBox1.Text +=
string.Format("{0,6:f1}", (AdjXYZ[j, 0] - StaKXYZ[i, 0]) * 1000) +
string.Format("{0,7:f1}", (AdjXYZ[j, 1] - StaKXYZ[i, 1]) * 1000) +
string.Format("{0,7:f1}", (AdjXYZ[j, 2] - StaKXYZ[i, 2]) * 1000) + "\n";
break;
}
}
}
//以下输出方差分量估计信息
richTextBox1.Text += "\n-------------------------------------------------------------\n";
if (IsD0)
{
for (i = 0; i < D0List.Count; i++)
{
StaData D = (StaData)StaList[i];
richTextBox1.Text += "测站 " + D.StaName + " 各次方差分量估计 各类观测值(L V S)的方差估值如下:\n";
richTextBox1.Text += "第*次 水平方向L 天顶距V 斜距S\n";
ArrayList Di = (ArrayList)D0List[i];
for (j = 0; j < Di.Count; j++)
{
DLVS dvs = (DLVS)Di[j];
richTextBox1.Text += string.Format("{0,3}", j + 1);
richTextBox1.Text += string.Format("{0,12:f3}", dvs.DL);
richTextBox1.Text += string.Format("{0,12:f3}", dvs.DV);
richTextBox1.Text += string.Format("{0,12:f3}", dvs.DS) + "\n";
}
richTextBox1.Text += "\n------------------------------------------------------\n";
}
}
else
{
richTextBox1.Text += "未进行方差分量估计!\n";
}
}
public void PrintResultabk(
double[,] AppXYZ, double[,] AdjXYZ, double[] StaZ, double[] StaAppZ, double[,] AB, double[] K,
ArrayList StaList, ArrayList D0List, bool IsD0,
string[] StaKName, double[,] StaKXYZ)
{
richTextBox1.Clear();
int i, j;
//先输出定向角未知数计算结果
richTextBox1.Text += "测站 定向角近似值(°′″) 定向角平差值(°′″) 改正数(″)\n";
for (i = 0; i < StaList.Count; i++)
{
StaData D = (StaData)StaList[i];//获取各测站的数据,需要强制类型转换获取StaList中的数据
richTextBox1.Text +=
string.Format("{0,4}", D.StaName) +
string.Format("{0,17:f4}", StaAppZ[i]) +
string.Format("{0,23:f4}", StaZ[i]) +
string.Format("{0,17:f2}", (StaZ[i] - StaAppZ[i]) * 10000) + "\n";
}
richTextBox1.Text += "\n-------------------------------------------------\n\n";
//输出 固定误差 和 比例误差
richTextBox1.Text += "测站 斜距固定误差(mm) 斜距比例误差(mm/km)\n";
for (i = 0; i < StaList.Count; i++)
{
StaData D = (StaData)StaList[i];//获取各测站的数据,需要强制类型转换获取StaList中的数据
richTextBox1.Text +=
string.Format("{0,4}", D.StaName) +
string.Format("{0,13:f2}", AB[i, 0]) +
string.Format("{0,19:f2}", AB[i, 1]) + "\n";
}
richTextBox1.Text += "\n-------------------------------------------------\n\n";
//输出大气折光系数估值
richTextBox1.Text += "测站 大气折光系数估值\n";
for (i = 0; i < StaList.Count; i++)
{
StaData D = (StaData)StaList[i];//获取各测站的数据,需要强制类型转换获取StaList中的数据
richTextBox1.Text +=
string.Format("{0,4}", D.StaName) +
string.Format("{0,14:f2}", K[i]) + "\n";
}
richTextBox1.Text += "\n-------------------------------------------------\n\n";
richTextBox1.Text += " 既有坐标(0) | 近似坐标(1) | 平差坐标(2) | 坐标较差[(1)-(0)] | 坐标较差[(2)-(0)]\n";
richTextBox1.Text += "点名 | X(m) Y(m) Z(m) | X(m) Y(m) Z(m) | X(m) Y(m) Z(m) | dx(mm) dy(mm) dz(mm) | dx(mm) dy(mm) dz(mm)\n";
for (i = 0; i < StaKName.Length; i++)
{
richTextBox1.Text +=
string.Format("{0,4}", StaKName[i]) + " | " +
string.Format("{0:0.0000}", StaKXYZ[i, 0]) + " " +
string.Format("{0:0.0000}", StaKXYZ[i, 1]) + " " +
string.Format("{0:0.0000}", StaKXYZ[i, 2]) + " | ";
for (j = 0; j < StaList.Count; j++)
{
StaData D = (StaData)StaList[j];
if (D.StaName == StaKName[i])
{
richTextBox1.Text +=
string.Format("{0:0.0000}", AppXYZ[j, 0]) + " " +
string.Format("{0:0.0000}", AppXYZ[j, 1]) + " " +
string.Format("{0:0.0000}", AppXYZ[j, 2]) + " | ";
richTextBox1.Text +=
string.Format("{0:0.0000}", AdjXYZ[j, 0]) + " " +
string.Format("{0:0.0000}", AdjXYZ[j, 1]) + " " +
string.Format("{0:0.0000}", AdjXYZ[j, 2]) + " |";
richTextBox1.Text +=
string.Format("{0,6:f1}", (AppXYZ[j, 0] - StaKXYZ[i, 0]) * 1000) +
string.Format("{0,7:f1}", (AppXYZ[j, 1] - StaKXYZ[i, 1]) * 1000) +
string.Format("{0,7:f1}", (AppXYZ[j, 2] - StaKXYZ[i, 2]) * 1000) + " |";
richTextBox1.Text +=
string.Format("{0,6:f1}", (AdjXYZ[j, 0] - StaKXYZ[i, 0]) * 1000) +
string.Format("{0,7:f1}", (AdjXYZ[j, 1] - StaKXYZ[i, 1]) * 1000) +
string.Format("{0,7:f1}", (AdjXYZ[j, 2] - StaKXYZ[i, 2]) * 1000) + "\n";
break;
}
}
}
//以下输出方差分量估计信息
richTextBox1.Text += "\n-------------------------------------------------------------\n";
if (IsD0)
{
for (i = 0; i < D0List.Count; i++)
{
StaData D = (StaData)StaList[i];
richTextBox1.Text += "测站 " + D.StaName + " 各次方差分量估计 各类观测值(L V S)的方差估值如下:\n";
richTextBox1.Text += "第*次 水平方向L 天顶距V 斜距S\n";
ArrayList Di = (ArrayList)D0List[i];
for (j = 0; j < Di.Count; j++)
{
DLVS dvs = (DLVS)Di[j];
richTextBox1.Text += string.Format("{0,3}", j + 1);
richTextBox1.Text += string.Format("{0,12:f3}", dvs.DL);
richTextBox1.Text += string.Format("{0,12:f3}", dvs.DV);
richTextBox1.Text += string.Format("{0,12:f3}", dvs.DS) + "\n";
}
richTextBox1.Text += "\n------------------------------------------------------\n";
}
}
else
{
richTextBox1.Text += "未进行方差分量估计!\n";
}
}
public void PrintResultab(
double[,] AppXYZ, double[,] AdjXYZ, double[] StaZ, double[] StaAppZ,double[,] AB,
ArrayList StaList, ArrayList D0List, bool IsD0,
string[] StaKName, double[,] StaKXYZ)
{
richTextBox1.Clear();
int i, j;
//先输出定向角未知数计算结果
richTextBox1.Text += "测站 定向角近似值(°′″) 定向角平差值(°′″) 改正数(″)\n";
for (i = 0; i < StaList.Count; i++)
{
StaData D = (StaData)StaList[i];//获取各测站的数据,需要强制类型转换获取StaList中的数据
richTextBox1.Text +=
string.Format("{0,4}", D.StaName) +
string.Format("{0,17:f4}", StaAppZ[i]) +
string.Format("{0,23:f4}", StaZ[i]) +
string.Format("{0,17:f2}", (StaZ[i] - StaAppZ[i]) * 10000) + "\n";
}
richTextBox1.Text += "\n-------------------------------------------------\n\n";
//输出 固定误差 和 比例误差
richTextBox1.Text += "测站 斜距固定误差(mm) 斜距比例误差(mm/km)\n";
for (i = 0; i < StaList.Count; i++)
{
StaData D = (StaData)StaList[i];//获取各测站的数据,需要强制类型转换获取StaList中的数据
richTextBox1.Text +=
string.Format("{0,4}", D.StaName) +
string.Format("{0,13:f2}", AB[i,0]) +
string.Format("{0,19:f2}", AB[i,1]) +"\n";
}
richTextBox1.Text += "\n-------------------------------------------------\n\n";
richTextBox1.Text += " 既有坐标(0) | 近似坐标(1) | 平差坐标(2) | 坐标较差[(1)-(0)] | 坐标较差[(2)-(0)]\n";
richTextBox1.Text += "点名 | X(m) Y(m) Z(m) | X(m) Y(m) Z(m) | X(m) Y(m) Z(m) | dx(mm) dy(mm) dz(mm) | dx(mm) dy(mm) dz(mm)\n";
for (i = 0; i < StaKName.Length; i++)
{
richTextBox1.Text +=
string.Format("{0,4}", StaKName[i]) + " | " +
string.Format("{0:0.0000}", StaKXYZ[i, 0]) + " " +
string.Format("{0:0.0000}", StaKXYZ[i, 1]) + " " +
string.Format("{0:0.0000}", StaKXYZ[i, 2]) + " | ";
for (j = 0; j < StaList.Count; j++)
{
StaData D = (StaData)StaList[j];
if (D.StaName == StaKName[i])
{
richTextBox1.Text +=
string.Format("{0:0.0000}", AppXYZ[j, 0]) + " " +
string.Format("{0:0.0000}", AppXYZ[j, 1]) + " " +
string.Format("{0:0.0000}", AppXYZ[j, 2]) + " | ";
richTextBox1.Text +=
string.Format("{0:0.0000}", AdjXYZ[j, 0]) + " " +
string.Format("{0:0.0000}", AdjXYZ[j, 1]) + " " +
string.Format("{0:0.0000}", AdjXYZ[j, 2]) + " |";
richTextBox1.Text +=
string.Format("{0,6:f1}", (AppXYZ[j, 0] - StaKXYZ[i, 0]) * 1000) +
string.Format("{0,7:f1}", (AppXYZ[j, 1] - StaKXYZ[i, 1]) * 1000) +
string.Format("{0,7:f1}", (AppXYZ[j, 2] - StaKXYZ[i, 2]) * 1000) + " |";
richTextBox1.Text +=
string.Format("{0,6:f1}", (AdjXYZ[j, 0] - StaKXYZ[i, 0]) * 1000) +
string.Format("{0,7:f1}", (AdjXYZ[j, 1] - StaKXYZ[i, 1]) * 1000) +
string.Format("{0,7:f1}", (AdjXYZ[j, 2] - StaKXYZ[i, 2]) * 1000) + "\n";
break;
}
}
}
//以下输出方差分量估计信息
richTextBox1.Text += "\n-------------------------------------------------------------\n";
if (IsD0)
{
for (i = 0; i < D0List.Count; i++)
{
StaData D = (StaData)StaList[i];
richTextBox1.Text += "测站 " + D.StaName + " 各次方差分量估计 各类观测值(L V S)的方差估值如下:\n";
richTextBox1.Text += "第*次 水平方向L 天顶距V 斜距S\n";
ArrayList Di = (ArrayList)D0List[i];
for (j = 0; j < Di.Count; j++)
{
DLVS dvs = (DLVS)Di[j];
richTextBox1.Text += string.Format("{0,3}", j + 1);
richTextBox1.Text += string.Format("{0,12:f3}", dvs.DL);
richTextBox1.Text += string.Format("{0,12:f3}", dvs.DV);
richTextBox1.Text += string.Format("{0,12:f3}", dvs.DS) + "\n";
}
richTextBox1.Text += "\n------------------------------------------------------\n";
}
}
else
{
richTextBox1.Text += "未进行方差分量估计!\n";
}
}
public void PrintResult(
double[,] AppXYZ, double[,] AdjXYZ, double[] StaZ, double[] StaAppZ,
ArrayList StaList,ArrayList D0List,bool IsD0,
string[] StaKName, double[,] StaKXYZ)
{
richTextBox1.Clear();
int i,j;
//先输出定向角未知数计算结果
richTextBox1.Text += "测站 定向角近似值(°′″) 定向角平差值(°′″) 改正数(″)\n";
for (i = 0; i < StaList.Count; i++)
{
StaData D = (StaData)StaList[i];//获取各测站的数据,需要强制类型转换获取StaList中的数据
richTextBox1.Text +=
string.Format("{0,4}",D.StaName) +
string.Format("{0,17:f4}",StaAppZ[i])+
string.Format("{0,23:f4}",StaZ[i]) +
string.Format("{0,17:f2}",(StaZ[i]-StaAppZ[i])*10000)+"\n";
}
richTextBox1.Text +="\n-------------------------------------------------\n\n";
richTextBox1.Text += " 既有坐标(0) | 近似坐标(1) | 平差坐标(2) | 坐标较差[(1)-(0)] | 坐标较差[(2)-(0)]\n";
richTextBox1.Text += "点名 | X(m) Y(m) Z(m) | X(m) Y(m) Z(m) | X(m) Y(m) Z(m) | dx(mm) dy(mm) dz(mm) | dx(mm) dy(mm) dz(mm)\n";
for (i = 0; i < StaKName.Length; i++)
{
richTextBox1.Text +=
string.Format("{0,4}", StaKName[i]) + " | " +
string.Format("{0:0.0000}", StaKXYZ[i, 0]) + " " +
string.Format("{0:0.0000}", StaKXYZ[i, 1]) + " " +
string.Format("{0:0.0000}", StaKXYZ[i, 2]) + " | ";
for (j = 0; j < StaList.Count;j++)
{
StaData D = (StaData)StaList[j];
if (D.StaName == StaKName[i])
{
richTextBox1.Text +=
string.Format("{0:0.0000}", AppXYZ[j, 0]) + " " +
string.Format("{0:0.0000}", AppXYZ[j, 1]) + " " +
string.Format("{0:0.0000}", AppXYZ[j, 2]) + " | ";
richTextBox1.Text +=
string.Format("{0:0.0000}", AdjXYZ[j, 0]) + " " +
string.Format("{0:0.0000}", AdjXYZ[j, 1]) + " " +
string.Format("{0:0.0000}", AdjXYZ[j, 2]) + " |";
richTextBox1.Text +=
string.Format("{0,6:f1}", (AppXYZ[j, 0] - StaKXYZ[i, 0]) * 1000) +
string.Format("{0,7:f1}", (AppXYZ[j, 1] - StaKXYZ[i, 1]) * 1000) +
string.Format("{0,7:f1}", (AppXYZ[j, 2] - StaKXYZ[i, 2]) * 1000) + " |";
richTextBox1.Text +=
string.Format("{0,6:f1}", (AdjXYZ[j, 0] - StaKXYZ[i, 0]) * 1000) +
string.Format("{0,7:f1}", (AdjXYZ[j, 1] - StaKXYZ[i, 1]) * 1000) +
string.Format("{0,7:f1}", (AdjXYZ[j, 2] - StaKXYZ[i, 2]) * 1000) + "\n";
break;
}
}
}
//以下输出方差分量估计信息
richTextBox1.Text += "\n-------------------------------------------------------------\n";
if (IsD0)
{
for (i = 0; i < D0List.Count;i++ )
{
StaData D = (StaData)StaList[i];
richTextBox1.Text+="测站 "+D.StaName+" 各次方差分量估计 各类观测值(L V S)的方差估值如下:\n";
richTextBox1.Text += "第*次 水平方向L 天顶距V 斜距S\n";
ArrayList Di = (ArrayList)D0List[i];
for (j = 0; j < Di.Count;j++ )
{
DLVS dvs = (DLVS)Di[j];
richTextBox1.Text +=string.Format("{0,3}",j + 1);
richTextBox1.Text += string.Format("{0,12:f3}", dvs.DL);
richTextBox1.Text += string.Format("{0,12:f3}", dvs.DV);
richTextBox1.Text += string.Format("{0,12:f3}", dvs.DS) + "\n";
}
richTextBox1.Text +="\n------------------------------------------------------\n";
}
}
else
{
richTextBox1.Text += "未进行方差分量估计!\n";
}
}
public void FillNSk(double[] NS, double[] WS, int Index, StaData D,//Index表示第几个测站
double[] PS, string[] KPName, double[,] KPXYZ, double[,] StaAppXYZ,
double[,] BS, double[] lS)
{
//用距离观测值填充NS和WS矩阵
double Stax = StaAppXYZ[Index, 0];
double Stay = StaAppXYZ[Index, 1];
double Staz = StaAppXYZ[Index, 2];
double x = 0, y = 0, z = 0;
double a, b, c, S0, ls;//dx,dy,dz误差方程系数,近似距离,距离常数项
for (int i = 0; i < D.Num; i++)
{
GetXYZ(D.CP3Name[i], ref x, ref y, ref z, KPName, KPXYZ);
S0 = Math.Sqrt(Math.Pow(x - Stax, 2) + Math.Pow(y - Stay, 2) + Math.Pow(z - Staz, 2));
ls = D.LVS[i, 2] - S0;
a = (Stax - x) / S0;
b = (Stay - y) / S0;
c = (Staz - z) / S0;
BS[i, 0] = a;
BS[i, 1] = b;
BS[i, 2] = c;
NS[5] += PS[i] * a * a;
NS[8] += PS[i] * a * b;
NS[9] += PS[i] * b * b;
NS[12] += PS[i] * a * c;
NS[13] += PS[i] * b * c;
NS[14] += PS[i] * c * c;
ls *= 1000;
lS[i] = ls;
ls *= PS[i];
WS[2] += a * ls;
WS[3] += b * ls;
WS[4] += c * ls;
}
}
public void FillNSab(double[] NS, double[] WS, int Index, StaData D,//Index表示第几个测站
double[] PS, string[] KPName, double[,] KPXYZ, double[,] StaAppXYZ,
double[,] BS, double[] lS)
{
//用距离观测值填充NS和WS矩阵
double Stax = StaAppXYZ[Index, 0];
double Stay = StaAppXYZ[Index, 1];
double Staz = StaAppXYZ[Index, 2];
double x = 0, y = 0, z = 0;
double a, b, c, S0, S0Km,ls;//dx,dy,dz误差方程系数,近似距离,距离常数项
for (int i = 0; i < D.Num; i++)
{
GetXYZ(D.CP3Name[i], ref x, ref y, ref z, KPName, KPXYZ);
S0 = Math.Sqrt(Math.Pow(x - Stax, 2) + Math.Pow(y - Stay, 2) + Math.Pow(z - Staz, 2));
ls = D.LVS[i, 2] - S0;
a = (Stax - x) / S0;
b = (Stay - y) / S0;
c = (Staz - z) / S0;
BS[i, 0] = a;
BS[i, 1] = b;
BS[i, 2] = c;
S0Km = S0 / 1000;//将斜距单位从m转换为km
BS[i, 3] = S0Km;
NS[2] += PS[i] * a * a;
NS[4] += PS[i] * a * b;
NS[5] += PS[i] * b * b;
NS[7] += PS[i] * a * c;
NS[8] += PS[i] * b * c;
NS[9] += PS[i] * c * c;
NS[11] -= PS[i] * a;
NS[12] -= PS[i] * b;
NS[13] -= PS[i] * c;
NS[14] += PS[i];
NS[16] -= PS[i] * a * S0Km;
NS[17] -= PS[i] * b * S0Km;
NS[18] -= PS[i] * c * S0Km;
NS[19] += PS[i] * S0Km;
NS[20] += PS[i] * S0Km * S0Km;
ls *= 1000;
lS[i] = ls;
ls *= PS[i];
WS[1] += a * ls;
WS[2] += b * ls;
WS[3] += c * ls;
WS[4] -= ls;
WS[5] -= S0Km * ls;
}
}
public void FillNS(double[] NS, double[] WS, int Index,StaData D,//Index表示第几个测站
double[] PS,string[] KPName, double[,] KPXYZ, double[,] StaAppXYZ,
double[,] BS,double[] lS)
{
//用距离观测值填充NS和WS矩阵
double Stax = StaAppXYZ[Index, 0];
double Stay = StaAppXYZ[Index, 1];
double Staz = StaAppXYZ[Index, 2];
double x = 0, y = 0, z = 0;
double a, b, c,S0,ls;//dx,dy,dz误差方程系数,近似距离,距离常数项
for (int i = 0; i < D.Num; i++)
{
GetXYZ(D.CP3Name[i], ref x, ref y, ref z, KPName, KPXYZ);
S0=Math.Sqrt(Math.Pow(x-Stax,2)+Math.Pow(y-Stay,2)+Math.Pow(z-Staz,2));
ls = D.LVS[i, 2] - S0;
a = (Stax - x) / S0;
b = (Stay - y) / S0;
c = (Staz - z) / S0;
BS[i, 0] = a;
BS[i, 1] = b;
BS[i, 2] = c;
NS[2] += PS[i] * a * a;
NS[4] += PS[i] * a * b;
NS[5] += PS[i] * b * b;
NS[7] += PS[i] * a * c;
NS[8] += PS[i] * b * c;
NS[9] += PS[i] * c * c;
ls *= 1000;
lS[i] = ls;
ls *=PS[i];
WS[1] += a * ls;
WS[2] += b * ls;
WS[3] += c * ls;
}
}
public void CalPolXYZ(StaData D,double[] XY,double[] dH)//计算测站坐标系下各点的坐标
{
int i;
double L,V,S,Sij;
for (i=0;i<D.Num;i++)
{
L=D.LVS[i,0]*Pe; //度转换为弧度
V=D.LVS[i,1]*Pe;//度转换为弧度
S=D.LVS[i,2];
Sij=S*Math.Sin(V);
dH[i]=-S*Math.Cos(V);//dH
XY[2*i]=Sij*Math.Cos(L);
XY[2*i+1]=Sij*Math.Sin(L);
}
}
public void CalNewXY(StaData D, double[] NewXY, string[] KPName, double[,] KPXYZ, double[] dH)
{
int i,j;
for (i=0;i<D.Num;i++)
{
for (j=0;j<KPName.Length;j++)
{
if (D.CP3Name[i]==KPName[j])
{
NewXY[2*i]=KPXYZ[j,0];
NewXY[2*i+1]=KPXYZ[j,1];
dH[i] += KPXYZ[j, 2];
break;
}
}
}
}
public void CopyStaAppXY(double[,] XYZ,double[,] XYZCopy)//m-行数;n-列数
{
int i,j;
for (i=0;i<XYZ.GetLength(0);i++)
{
for (j = 0; j < XYZ.GetLength(1); j++)
{
XYZCopy[i, j] = XYZ[i, j];
}
}
}
public void CalStaAppXY(double[,] StaAppXYZ, ArrayList StaList, string[] KPName, double[,] KPXYZ)
{
double[] FourPara=new double[4];
TransXY T = new TransXY();
int i,j,n;
for (i = 0; i < StaList.Count; i++)
{
StaData D = (StaData)StaList[i];
n = D.Num;
double[] NewXY = new double[2 * n];//x,y……x,y存储
double[] OldXY = new double[2 * n];//x,y……x,y存储
double[] dH = new double[n];
CalPolXYZ(D, OldXY, dH);
CalNewXY(D, NewXY, KPName, KPXYZ,dH);
T.Cal_FourPara(NewXY, OldXY, ref FourPara);//没写完
StaAppXYZ[i, 0] = FourPara[0];
StaAppXYZ[i, 1] = FourPara[1];
for (j = 0; j < n; j++)
{
StaAppXYZ[i, 2] +=dH[j];
}
StaAppXYZ[i, 2] /= n;
}
}
public void Cal_SWeight(double[] P2, double ma, double a, double b,StaData D,bool Type)
{
//Type=true:a+b*S; Type=false:sqrt(a*a+(b*S)^2)
int i;
double m2=ma*ma;
b/=1000;
double a2=a*a;
if (Type==false)
{
for (i=0;i<D.Num;i++)
{
P2[i]=m2/(a2+Math.Pow(b*D.LVS[i,2],2));//用仪器标称精度初始化边长权矩阵
}
}
else
{
for (i=0;i<D.Num;i++)
{
P2[i] = Math.Pow(ma / (a + b * D.LVS[i, 2]), 2);//用仪器标称精度初始化边长权矩阵
}
}
}
public string GetStaName(string str)
{
int i,j,k;
string s1 = "站";
string s0="";
i = str.IndexOf(s1)+s1.Length;
j = str.IndexOf("的");
for (k=i;k<j;k++)
{
s0 += str[k];
}
return s0;
}
public string GetData(string str, ref int m)
{
string s0 = "";
for (int i = m; i < str.Length; i++)
{
if (str[i] != ',')
{
s0 += str[i];
}
else
{
m = i + 1;
break;
}
}
return s0;
}
public int GetStrIndex(string str, int m)
{
int i;
for (i = m; i < str.Length; i++)
{
if (str[i]!=' '&&str[i]!='\t')
{
return i;
}
}
return i;
}
public string GetLVSData(string str, ref int m)
{
string s0 = "";
for (int i = m; i < str.Length; i++)
{
if (str[i] != ' ' && str[i] != '\t')
{
s0+=str[i];
}
else
{
m = i;
break;
}
}
return s0;
}
public void GetLVS(string str,ref string s0,ref double L,ref double V,ref double S)
{
str = str.Trim();
int m=0;
s0=GetLVSData(str, ref m);
m = GetStrIndex(str, m);
L =double.Parse(GetLVSData(str, ref m));
m = GetStrIndex(str, m);
V =double.Parse(GetLVSData(str, ref m));
m = GetStrIndex(str, m);
S = double.Parse(GetLVSData(str, ref m));
}
private void 导入CPIII数据ToolStripMenuItem1_Click(object sender, EventArgs e)
{
OpenFileDialog openFileDialog1 = new OpenFileDialog();
openFileDialog1.Multiselect = false;
openFileDialog1.Title = "导入CPIII已知坐标";
openFileDialog1.Filter = "文本文档(*.txt;*.csv)|*.txt;*.csv"
+ "|所有文件(*.*)|*.*";
openFileDialog1.InitialDirectory = @Path;
openFileDialog1.RestoreDirectory = true;
if (openFileDialog1.ShowDialog() == DialogResult.OK)
{
richTextBox1.Clear();
string fname = openFileDialog1.FileName;
StreamReader sr = new StreamReader(fname, Encoding.Default);
richTextBox1.Text = sr.ReadToEnd().ToString();
sr.Close();
String line;
int linecount = 0;
StreamReader A0 = new StreamReader(fname);
while ((line = A0.ReadLine())!= null)
{
if (line.Trim().Length>0)
{
linecount++;
}
}
A0.Close();
KPName = new string[linecount];
KPXYZ = new double[linecount, 3];
int k = 0, j = 0;//j记录逗号下标+1
StreamReader A1 = new StreamReader(fname);
while ((line = A1.ReadLine()) != null)
{
if (line.Trim().Length>0)
{
KPName[k] = GetData(line, ref j);
KPXYZ[k, 1] =double.Parse(GetData(line, ref j));//原文件格式是YXZ
KPXYZ[k, 0] =double.Parse(GetData(line, ref j));//原文件格式是YXZ
KPXYZ[k, 2] = double.Parse(GetData(line, ref j));
k++;
j = 0;
}
}
A1.Close();
}
}
private void 导入观测数据ToolStripMenuItem_Click(object sender, EventArgs e)
{
OpenFileDialog openFileDialog1 = new OpenFileDialog();
openFileDialog1.Multiselect = false;
openFileDialog1.Title = "导入观测文件";
openFileDialog1.Filter = "文本文档(*.txt;*.csv)|*.txt;*.csv"
+ "|所有文件(*.*)|*.*";
openFileDialog1.InitialDirectory = @Path;
openFileDialog1.RestoreDirectory = true;
if (openFileDialog1.ShowDialog() == DialogResult.OK)
{
richTextBox1.Clear();
string fname = openFileDialog1.FileName;
StreamReader sr = new StreamReader(fname, Encoding.UTF8);
richTextBox1.Text = sr.ReadToEnd().ToString();
sr.Close();
String line;
int linecount = 0;
StreamReader A0 = new StreamReader(fname);
while ((line = A0.ReadLine()) != null)
{
if (line.Trim().Length > 0)
{
linecount++;
}
}
A0.Close();
int k = 0;
StaList.Clear();//清空原来的数据
StreamReader A1 = new StreamReader(fname);
line = A1.ReadLine();
while (line!=null)
{
if (line.Trim().Length > 0)
{
if (line[3]=='-')
{
StaData D=new StaData();
D.StaName=GetStaName(line);
line = A1.ReadLine();//跳过"序号 水平方向观测值 天顶距 斜距"
while ((line = A1.ReadLine()) != null)
{
line = line.Trim();
if (line.Length > 0)
{
if (line[3]!='-')
{
GetLVS(line, ref D.CP3Name[k], ref D.LVS[k, 0], ref D.LVS[k, 1], ref D.LVS[k, 2]);
D.LVS[k, 0] *= 0.9; //gon -> ″
D.LVS[k, 1] *= 0.9; //gon -> ″
D.LVS[k,2]-=0.0344;//减去棱镜常数
k++;
}
else
{
D.Num=k;
StaList.Add(D);
k = 0;
break;
}
}
}
if (k>0)//最后一个测站的数据
{
D.Num = k;
StaList.Add(D);
}
}
}
else
{
line = A1.ReadLine();
}
}
A1.Close();
}
}
private void 导入测站ToolStripMenuItem_Click(object sender, EventArgs e)
{
OpenFileDialog openFileDialog1 = new OpenFileDialog();
openFileDialog1.Multiselect = false;
openFileDialog1.Title = "导入既有测站坐标";
openFileDialog1.Filter = "文本文档(*.txt;*.dat)|*.txt;*.dat"
+ "|所有文件(*.*)|*.*";
openFileDialog1.InitialDirectory = @Path;
openFileDialog1.RestoreDirectory = true;
if (openFileDialog1.ShowDialog() == DialogResult.OK)
{
richTextBox1.Clear();
string fname = openFileDialog1.FileName;
StreamReader sr = new StreamReader(fname, Encoding.UTF8);
richTextBox1.Text = sr.ReadToEnd().ToString();
sr.Close();
String line;
int linecount = 0;
StreamReader A0 = new StreamReader(fname);
while ((line = A0.ReadLine()) != null)
{
if (line.Trim().Length > 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)
{