化学实验的七参数与点坐标，你真的了解吗？

完整代码

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
namespace CoordTrans
{
    public class MyPoint
    {
        public double L =0, B=0, H=0; //大地坐标经度、维度、高度
        public double X=0, Y = 0, Z = 0; //空间坐标系
        public double X2 = 0, Y2 = 0, Z2 = 0;  //七参数转换后的空间坐标
        public double a = 0, f = 0, b = 0, e = 0, e2 = 0;    //椭球参数
        public double rho = 180 / Math.PI, d2r = Math.PI / 180;
        public double xs, ys, hs;
        public double dx = 0, dy = 0, dz = 0, rx = 0, ry = 0, rz = 0, m = 0, k=0;    //七参数 三个线性平移量-单位米 三个旋转平移量-十进制秒为单位（运算时注意转换为度） 比例因子-单位百万分率 (ppm)
                                                                         //测量队给出的七参数单位与计算的单位不同，要进行单位转化 1 秒=0.0000048481373323 弧度
                                                                         //尺度因子有两种单位的表示形式，一种结果约为1，如1.0000045，用k表示；
                                                                         //另一种就是ppm的表示形式，稍微比1大一点，如4.5，用m表示。k=m/1000000
        public void SetLBH_Degree(double L_, double B_, double H_=0)
        {
            // L-经度 B-维度
            //GPS给出的经纬度常用格式有 度.分分秒秒  度.度
            this.L = L_;
            this.B = B_;
            this.H = H_;            
        }
        public void SetLBH_DMS(string L_string, string B_string, string H_string="0")
        {
            //经纬度格式说明 小数点前面是度，小数点后面1、2位是分，3、4位是秒，5位以后是秒的小数部分
            // 解析经度字符串
            int index_L_point = L_string.IndexOf(".");
            string L_deg_str = L_string.Substring(0, index_L_point);
            string L_mmss_str = L_string.Substring(index_L_point + 1);
            double degree_L = double.Parse(L_deg_str);
            string L_min_str = L_mmss_str.Substring(0, 2);
            double minutes_L = double.Parse(L_min_str);
            string L_sec_str1 = L_mmss_str.Substring(2, 2);
            string L_sec_str2 = "0." + L_mmss_str.Substring(4);
            double second_L1 = double.Parse(L_sec_str1);
            double second_L2 = double.Parse(L_sec_str2);
            double second_L = second_L1 + second_L2;
            // 解析纬度字符串
            int index_B_point = B_string.IndexOf(".");
            string B_deg_str = B_string.Substring(0, index_B_point);
            string B_mmss_str = B_string.Substring(index_B_point + 1);
            double degree_B = double.Parse(B_deg_str);
            string B_min_str = B_mmss_str.Substring(0, 2);
            double minutes_B = double.Parse(B_min_str);
            string B_sec_str1 = B_mmss_str.Substring(2, 2);
            string B_sec_str2 = "0." + B_mmss_str.Substring(4);
            double second_B1 = double.Parse(B_sec_str1);
            double second_B2 = double.Parse(B_sec_str2);
            double second_B = second_B1 + second_B2;
            this.L = (degree_L + minutes_L / 60.0 + second_L / 3600.0);
            this.B = (degree_B + minutes_B / 60.0 + second_B / 3600.0);
            this.H = double.Parse(H_string);
        }
        public void BLH2XYZ(string type="WGS-84")
        {
            EllipsoidParam(type);
            double sB = Math.Sin(this.B * d2r);
            double cB = Math.Cos(this.B * d2r);
            double sL = Math.Sin(this.L * d2r);
            double cL = Math.Cos(this.L * d2r);
            double N = this.a / (Math.Sqrt(1 - this.e * this.e * sB * sB));
            this.X = (N + this.H) * cB * cL;
            this.Y = (N + this.H) * cB * sL;
            this.Z = (N * (1 - this.e * this.e) + this.H) * sB;
            this.X2 = this.X;
            this.Y2 = this.Y;
            this.Z2 = this.Z;
        }
        public void XYZ2BLH(string type)
        {
            EllipsoidParam(type);
            // 这里转出来的B L是弧度
            this.L = Math.Atan(this.Y2 / this.X2) + Math.PI;
            //this.L = Math.Atan(this.Y2 / this.X2);
            this.L = this.L * 180 / Math.PI;
            // B需要迭代计算
            double B2 = Math.Atan(Z2 / Math.Sqrt(X2 * X2 + Y2 * Y2));
            double B1;
            double N;
            while (true)
            {
                N = a / Math.Sqrt(1 - f * (2 - f) * Math.Sin(B2) * Math.Sin(B2));
                B1 = Math.Atan((Z2 + N * f * (2 - f) * Math.Sin(B2)) / Math.Sqrt(X2 * X2 + Y2 * Y2));
                if (Math.Abs(B1 - B2) < 0.00000001)
                    break;
                B2 = B1;
            }
            this.B = B2 * 180 / Math.PI ;
            double sB = Math.Sin(this.B * d2r);
            double cB = Math.Cos(this.B * d2r);
            this.H = this.Z2 / sB - N * (1 - this.e * this.e);
        }
        
        public void SevenParamTrans()
        {
            SetSevenParamUnitTrans();
            this.X2 = (1 + k) * (this.X + this.rz * this.Y - this.ry * this.Z) + this.dx;
            this.Y2 = (1 + k) * (-this.rz * this.X + this.Y + this.rx * this.Z) + this.dy;
            this.Z2 = (1 + k) * (this.ry * this.X - this.rx * this.Y + this.Z) + this.dz;
        }
        public void SetSevenParamUnitTrans()
        {
            // 七参数
            this.dx = 123.456;
            this.dy = 78.9999;
            this.dz = 24.5678;
            this.rx = -123.3333333 * 0.0000048481373323;
            this.ry = 338.88888888 * 0.0000048481373323;
            this.rz = 145.5555555 * 0.0000048481373323;
            this.m = 3.6666666;
            this.k = this.m / 1000000;
        }
        public void GaussProjection(string type, double L0)
        {
            this.EllipsoidParam(type);
            double l = (this.L - L0) / this.rho;
            double cB = Math.Cos(this.B * this.d2r);
            double sB = Math.Sin(this.B * this.d2r);
            double s2b = Math.Sin(this.B * this.d2r * 2);
            double s4b = Math.Sin(this.B * this.d2r * 4);
            double s6b = Math.Sin(this.B * this.d2r * 6);
            double s8b = Math.Sin(this.B * this.d2r * 8);
            double N = this.a / Math.Sqrt(1 - this.e * this.e * sB * sB);       // 卯酉圈曲率半径
            double t = Math.Tan(this.B * this.d2r);
            double eta = this.e2 * cB;
            double m0 = this.a * (1 - this.e * this.e);
            double m2 = 3.0 / 2.0 * this.e * this.e * m0;
            double m4 = 5.0 / 4.0 * this.e * this.e * m2;
            double m6 = 7.0 / 6.0 * this.e * this.e * m4;
            double m8 = 9.0 / 8.0 * this.e * this.e * m6;
            double a0 = m0 + 1.0 / 2.0 * m2 + 3.0 / 8.0 * m4 + 5.0 / 16.0 * m6 + 35.0 / 128.0 * m8;
            double a2 = 1.0 / 2.0 * m2 + 1.0 / 2.0 * m4 + 15.0 / 32.0 * m6 + 7.0 / 16.0 * m8;
            double a4 = 1.0 / 8.0 * m4 + 3.0 / 16.0 * m6 + 7.0 / 32.0 * m8;
            double a6 = 1.0 / 32.0 * m6 + 1.0 / 16.0 * m8;
            double a8 = 1.0 / 128.0 * m8;
            // X1为自赤道量起的子午线弧长
            double X1 = a0 * (this.B * this.d2r) - 1.0 / 2.0 * a2 * s2b + 1.0 / 4.0 * a4 * s4b - 1.0 / 6.0 * a6 * s6b + 1.0 / 8.0 * a8 * s8b;
            this.xs = X1 + N / 2 * t * cB * cB * l * l + N / 24 * t * (5 - t * t + 9 * Math.Pow(eta, 2) + 4 * Math.Pow(eta, 4)) * Math.Pow(cB, 4) * Math.Pow(l, 4)
                  + N / 720 * t * (61 - 58 * t * t + Math.Pow(t, 4)) * Math.Pow(cB, 6) * Math.Pow(l, 6);
            //this.ys = N * cB * l + N / 6 * (1 - t * t + eta * eta) * Math.Pow(cB, 3) * Math.Pow(l, 3)
            //    + N / 120 * (5 - 18 * t * t + Math.Pow(t, 4) + 14 * Math.Pow(eta, 2) - 58 * eta * eta * t * t) * Math.Pow(cB, 5) * Math.Pow(l, 5) + 500000;
            this.ys = N * cB * l + N / 6 * (1 - t * t + eta * eta) * Math.Pow(cB, 3) * Math.Pow(l, 3)
                + N / 120 * (5 - 18 * t * t + Math.Pow(t, 4) + 14 * Math.Pow(eta, 2) - 58 * eta * eta * t * t) * Math.Pow(cB, 5) * Math.Pow(l, 5);
            this.hs = this.H;
            // 假东 假北偏移
            double pseudoNorth = -2354048.666;
            double pseudoEast = 234864.444;
            this.xs += pseudoNorth;
            this.ys += pseudoEast;
        }
        // 椭球参数设置
        public void EllipsoidParam(string type)
        {
            // CGCS2000 椭球参数
            if (type == "CGCS2000")
            {
                this.a = 6378137;
                this.f = 1 / 298.257222101;
            }
            // 西安 80
            else if (type == "西安80")
            {
                this.a = 6378140;
                this.f = 1 / 298.257;
            }
            // 北京 54
            else if (type == "北京54")
            {
                this.a = 6378245;
                this.f = 1 / 298.3;
            }
            // WGS-84 
            else
            {
                this.a = 6378137;
                this.f = 1 / 298.257223563;
            }
            this.b = this.a * (1 - this.f);
            this.e = Math.Sqrt(this.a * this.a - this.b * this.b) / this.a;  //第一偏心率
            this.e2 = Math.Sqrt(this.a * this.a - this.b * this.b) / this.b;  //第二偏心率
        }
    }
}

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Windows;
using System.Windows.Controls;
using System.Windows.Data;
using System.Windows.Documents;
using System.Windows.Input;
using System.Windows.Media;
using System.Windows.Media.Imaging;
using System.Windows.Navigation;
using System.Windows.Shapes;
namespace CoordTrans
{
    /// 

    /// MainWindow.xaml 的交互逻辑
    /// 
    public partial class MainWindow : Window
    {
        public MainWindow()
        {
            InitializeComponent();
            Run();
        }
        public void Run()
        {
            MyPoint p1 = new MyPoint();
            MyPoint p2 = new MyPoint();
            MyPoint p3 = new MyPoint();
            List pointList = new List { p1, p2, p3 };
            p1.SetLBH_DMS(113.33333333.ToString(), 22.44444444.ToString(), 119.123.ToString());
            p2.SetLBH_DMS(113.55555555.ToString(), 22.88888888.ToString(), 66.777.ToString());
            p3.SetLBH_DMS(113.66666666.ToString(), 22.55555555.ToString(), 77.888.ToString());
            foreach (MyPoint p in pointList)
            {
                // 经纬度转空间坐标
                p.BLH2XYZ();
                // 七参数转换
                p.SevenParamTrans();
                // 空间坐标转回经纬度
                p.XYZ2BLH("CGCS2000");
                // 高斯投影 经纬度转平面坐标
                p.GaussProjection("CGCS2000", 114);
            }
        }
    }
}

由于保密缘由gps什么时候要转换参数，七参数与点座标都是乱写的。

转换总共有4步，第一步：将GPS经纬度转空间座标，进行第一步转换时，椭球参数选择WGS-84（固定选择这个，不用改）。

第二步：七参数转换。

本文所有公式来自孔祥元的大地测量学基础第二版。

这一步转换注意选择椭球参数，测量队给出的七参数都是对应一个确定的椭球参数，要与七参数的椭球参数相一致。

七参数转换注意单位，测量队给的七参数gps什么时候要转换参数，三个平移量单位是米不用变，3个旋转偏移量单位是秒，要转成弧度能够带入运算。比例因子的单位ppm是百万分率，除1000000才是公式（1+m）中的m。七参数是用于空间坐标系转换，不要对高斯投影后的平面坐标系用。

第三步：转换后的空间坐标系转到经纬度。这一步是为了高斯投影，高斯投影的输入是经纬度座标。

B须要迭代估算

this.L = Math.Atan(this.Y2 / this.X2) + Math.PI;

L的估算中加Pi是为了把纬度转入我国的100多度这个范围。

第四步：经纬度做高斯投影。这一步也要用与七参数对应的椭球参数。

高斯投影的实现可参见博主另一篇文章。这里注意检测队给的七参数是否有假北、假东偏斜，如果没有假北、假东

this.xs = X1 + N / 2 * t * cB * cB * l * l + N / 24 * t * (5 - t * t + 9 * Math.Pow(eta, 2) + 4 * Math.Pow(eta, 4)) * Math.Pow(cB, 4) * Math.Pow(l, 4)
                  + N / 720 * t * (61 - 58 * t * t + Math.Pow(t, 4)) * Math.Pow(cB, 6) * Math.Pow(l, 6);
this.ys = N * cB * l + N / 6 * (1 - t * t + eta * eta) * Math.Pow(cB, 3) * Math.Pow(l, 3) + N / 120 * (5 - 18 * t * t + Math.Pow(t, 4) + 14 * Math.Pow(eta, 2) - 58 * eta * eta * t * t) * Math.Pow(cB, 5) * Math.Pow(l, 5) + 500000;

把y投影加500000是习惯做法，为了让y值为正。

如果给了假北、假东偏斜，这里y不用加500000，给x和y分别加上假北、假东偏斜

// 假东 假北偏移
double pseudoNorth = -2355555.777;
double pseudoEast = 234866.999;
this.xs += pseudoNorth;
this.ys += pseudoEast;

高程写0虽然不影响最后估算结果，但是博主没有严格证明