利用C#编写一个高斯正反算程序

卤雅少年

已于 2022-10-21 17:10:00 修改

阅读量5.9k

点赞数 14

分类专栏：测绘工程高斯正反算文章标签： c#

于 2022-07-11 18:57:23 首次发布

本文链接：https://blog.csdn.net/zj1638/article/details/125380593

版权

测绘工程同时被 2 个专栏收录

5 篇文章 25 订阅

订阅专栏

高斯正反算

2 篇文章 0 订阅

订阅专栏

一、代码界面展示

在这里插入图片描述
整个界面控件为tabControl，groupBox，label，textbox，comboBox，button，richTextBook。

二、代码运算结果展示

数据结果采用国家统一坐标（横坐标上加上500000m，在坐标前面冠以带号）。
内置了4种椭球参数（CGCS200,WGS84,西安80，克拉索夫斯基椭球体）。
可以计算三度带和六度带。
计算曲率半径功能是我用来做大地测量学作业用的，顺便写在这里的。

2.1 数据相互转换验证

在这里插入图片描述相互转换的结果肯定是一样的（我自己测试了10几次不同的数据，是没有问题的）。

2.2 数据结果正确性验证

2.2.1 高斯正算

数据结果的正确性需要使用其他的软件进行验证。
这里我选取的是https://blog.csdn.net/starmings/article/details/123034189计算的数据进行验证。
在这里插入图片描述
软件结果如下：
采用西安80椭球参数，高斯正算的结果与软件的结果是一致的。

2.2.2 高斯反算

在这里插入图片描述软件结果如下：

采用西安80椭球参数，高斯反算的结果与软件的结果是一致的。（我只保留了0.01，已经可以了。）

三、源代码

3.1 变量设置

        //大地经度 度分秒
        public double L;
        public double Ld;
        public double Lf;
        public double Lm;
        //大地纬度 度分秒
        public double B;
        public double Bd;
        public double Bf;
        public double Bm;
        double L0;//中央子午线度数
        double N1;//带号
        double a;//长半轴
        double b;//短半轴
        double e1;//第一偏心率平方
        double e2; //第二偏心率平方
        double m0,m2,m4,m6,m8;
        double a0, a2, a4, a6, a8;
        int i ;

3.2 高斯正算

   		    //m0-m8的值
            m0 = a * (1 - e1);
            m2 = 1.5 * m0 * e1;
            m4 = 1.25 * m2 * e1;
            m6 = (7.0 / 6) * m4 * e1;
            m8 = (9.0 / 8) * m6 * e1;
            //a0-a8的值
            a0 = m0 + m2 / 2.0 + (3.0 / 8) * m4 + (5.0 / 16) * m6 + (35.0 / 128) * m8;
            a2 = (1.0 / 2) * (m2 + m4) + (15.0 / 32) * m6 + (7.0 / 16) * m8;
            a4 = (1.0 / 8) * m4 + (3.0 / 16) * m6 + (7.0 / 32) * m8;
            a6 = (1.0 / 32) * m6 + (1.0 / 16) * m8;
            a8 = (1.0 / 128) * m8;
            //将大地经度转化为度数
            Ld = double.Parse(textBox1.Text);
            Lf = double.Parse(textBox2.Text);
            Lm = double.Parse(textBox3.Text);
            L = Ld + Lf / 60 + Lm / 3600;
            //将大地纬度转化为度数
            Bd = double.Parse(textBox4.Text);
            Bf = double.Parse(textBox5.Text);
            Bm = double.Parse(textBox6.Text);
            B = Bd + Bf / 60 + Bm / 3600;
            if (comboBox2.Text == "三度带")
            {
                N1 = Math.Truncate((L + 1.5) / 3);//带号
                L0 = 3 * N1;
             
            }
            if (comboBox2.Text == "六度带")
            {
                N1 = Math.Truncate(L / 6) + 1;
                L0 = 6 * N1 - 3;
             
            }
            
            //求椭球面上某点到赤道的子午线弧长，精度至0.001m
            double radB = (Math.PI / 180) * B; //将纬度转换为弧度
            double Δ1, Δ2, Δ3, Δ4, Δ5;
            Δ1 = a0 * radB;
            Δ2 = (a2 / 2) * Math.Sin(2 * radB);
            Δ3 = (a4 / 4) * Math.Sin(4 * radB);
            Δ4 = (a6 / 6) * Math.Sin(6 * radB);
            Δ5 = (a8 / 8) * Math.Sin(8 * radB);
            double X =Math.Round(Δ1 - Δ2 + Δ3 - Δ4 + Δ5, 8);
            //卯酉圈曲率半径N
            double N = a / Math.Sqrt(1 - e1 * Math.Sin(radB) * Math.Sin(radB));
            // tan(B)的值
            double t = Math.Tan(radB);
            // η的值
            double η =Math.Sqrt(e2) * Math.Cos(radB) ;
            //l的值
            double ρm = 180 * 60 * 60 / Math.PI;
            double l = ((L - L0) * 3600) / ρm;
            //高斯坐标中的x值和y值
            double x = Math.Round( X + (N / 2) * t * Math.Cos(radB) * Math.Cos(radB) * l * l + (N / 24) * t * (5 - t * t + 9 * η* η + 4 * Math.Pow(η,4)) * Math.Pow(Math.Cos(radB), 4) * Math.Pow(l, 4) + (N / 720) * t * (61 - 58 * t * t + Math.Pow(t, 4)) * Math.Pow(l, 6) * (Math.Pow(Math.Cos(radB), 6)),4);
            double y = Math.Round( N * Math.Cos(radB) * l + (N / 6.0) * (1 - t * t + η* η) * Math.Pow(Math.Cos(radB), 3) * Math.Pow(l, 3) + (N / 120.0) * (5 - 18 * t * t + Math.Pow(t, 4) + 14 * η * η - 58 * η * η * t * t) * Math.Pow(Math.Cos(radB), 5) * Math.Pow(l, 5),4);
            richTextBox1.AppendText(Environment.NewLine);
            richTextBox1.Text = "大地经度为L:" + Ld + "°" + Lf + "′" + Lm + "″" + "        " + "大地经度为B:" + Bd + "°" + Bf + "′" + Bm+"″" + "\r\n";
            richTextBox1.Text += "大地坐标转高斯平面坐标后：" + "\r\n";
            richTextBox1.Text += "高斯X坐标：" + x + "米" + "        " + "高斯Y坐标：" + +N1 + "" + (y+500000)+  "米" + "\r\n";

3.3 高斯反算

 //输入坐标x,y
            double x = double.Parse(textBox7.Text);
            double y1 = double.Parse(textBox8.Text);
            double y = y1-Math.Truncate(y1 / Math.Pow(10, 6)) * Math.Pow(10, 6)-500000;
            double[] Bf1 = new double[20];
            double[] F = new double[20];
            Bf1[0] = x / a0;//底点纬度
            F[0] = (-a2 / 2.0) * Math.Sin(2 * Bf1[0]) + a4 / 4.0 * Math.Sin(4 * Bf1[0]) - a6 / 6.0 * Math.Sin(6 * Bf1[0]);
            Bf1[1] = (x - F[0]) / a0;
            double Δd =Bf1[1] - Bf1[0];
            for ( i = 0;i<10 ; i++)
            {
                F[i] = (-a2 / 2.0) * Math.Sin(2 * Bf1[i]) + a4 / 4.0 * Math.Sin(4 * Bf1[i]) - a6 / 6.0 * Math.Sin(6 * Bf1[i]);
                Bf1[i + 1] = (x - F[i]) / a0;
                Δd = Bf1[i + 1] - Bf1[i];
                if (Δd < Math.Pow(10, -10))
                {
                    i = i + 1;
                    break;
                }
                
            }
            // tan(Bf)的值
            double tf = Math.Tan(Bf1[i]);
            // ηf的值
            double ηf = e2 * Math.Cos(Bf1[i]) * Math.Cos(Bf1[i]);
            //计算Mf的值
            double Mf = a * (1 - e1) / Math.Pow(1 - e1 * Math.Sin(Bf1[i]) * Math.Sin(Bf1[i]), 1.5);
            //计算Nf的值
            double Nf = a / Math.Sqrt(1 - e1 * Math.Sin(Bf1[i]) * Math.Sin(Bf1[i]));
            //计算y坐标的前2位
            int aa = Convert.ToInt32(Math.Truncate(y1 / Math.Pow(10, 6)));
            if (comboBox2.Text == "三度带")
            {
               
                L0 = 3 * aa;

            }
            if (comboBox2.Text == "六度带")
            {
              
                L0 = 6 * aa - 3;

            }
            //求大地经度L
            L = (Math.PI / 180) * L0 + y / (Nf * Math.Cos(Bf1[i])) - (1 + 2 * tf * tf + ηf ) * Math.Pow(y, 3) / (6 * Math.Pow(Nf, 3) * Math.Cos(Bf1[i])) + (5 + 28 * tf * tf + 24 * Math.Pow(tf, 4) + 6 * ηf + 8 * ηf  * Math.Pow(tf,4)) * Math.Pow(y, 5) / (120 * Math.Pow(Nf, 5) * Math.Cos(Bf1[i]));
         
          
            //转换经度
            double degrees = (180 / Math.PI) * L;//将弧度转化为度
            int d1 = Convert.ToInt32(Math.Truncate(degrees));//度
            int f1 = Convert.ToInt32(Math.Truncate((degrees - d1) * 60));//分
            double m1 = Math.Round(((degrees - d1) * 60 - f1) * 60,2);//秒
            if (m1 == 60)
            {
                m1 = 0;
                f1 += 1;
            }
            //求大地纬度B
            B = Bf1[i] - (tf * y * y) / (2 * Mf * Nf)+ tf * (5 + 3 * tf * tf + ηf  - 9 * ηf  * tf * tf) * Math.Pow(y, 4) / (24 * Mf * Math.Pow(Nf, 3)) + tf * (61 + 90 * tf * tf + 45 * Math.Pow(tf, 4)) * Math.Pow(y, 6) / (720 * Mf * Math.Pow(Nf, 5));
            double degrees1 = (180 / Math.PI) * B;//将弧度转化为度
            int d = Convert.ToInt32(Math.Truncate(degrees1));//度
            int f = Convert.ToInt32(Math.Truncate((degrees1 - d) * 60));//分
            double m = Math.Round(((degrees1 - d) * 60 - f) * 60,2);//秒
            if (m == 60)
            {
                m = 0;
                f += 1;
            }
            richTextBox1.AppendText(Environment.NewLine);
            richTextBox1.Text += "高斯平面坐标转大地坐标后：" + "\r\n";
            richTextBox1.Text +=  "大地坐标L：" + d1+"°   "+f1+"′  "+m1+"″   "+ "大地坐标B：" + d + "°   " + f + "′  " + m + "″   "+"\r\n";

四、结语

最开始写这个高斯正反算程序的时候，公式完全按照大地测量学实习指导书上的进行编写，写完过后进行数据验证，发现相互转换就不能实现，分位就出现了偏差，说明正算公式或者反算公式出现了问题，经过检查，正算公式修改成了与https://blog.csdn.net/du12300/article/details/109307386所展示的公式。反算公式是对的，主要底点纬度的精度，我没在CSDN上找到相关的高斯反算文章，需要你自己发力了。
这个程序没有很困难的地方，就是公式的编写，可能就是需要注意double的精度丢失对数据结果带来的影响。
如果你需要批量处理，其实也简单，把我那个水准或者导线的程序里面的批量导入部分copy过来，设置几个变量，再写几个for语句就完事。
我的其他文章：