基于MATLAB的批量3度带高斯正算（LB--xy）

---

一、高斯正算原理

（1）高斯-克吕格投影

高斯-克吕格投影是由德国数学家、物理学家、天文学家高斯于19 世纪20 年代拟定，后经德国大地测量学家克吕格于1912 年对投影公式加以补充，故称为高斯-克吕格投影，又名"等角横切椭圆柱投影”，是地球椭球面和平面间正形投影的一种（来自百度百科）。

我们知道，地球是一个椭球体，就拿参考椭球来说，其不考虑真实表面，要想将数学球形表面的地理要素高精度表示在平面地图上，就不得不与数学球体表面无法展开产生矛盾，而这也一直都是当时世界探讨的问题。后面一大批学者加入投影的研究中去，在诸多的成果中，根据我国处于的地理位置，确定我国统一高斯-克吕格投影的方式对中国片区进行投影。其又名“等角横切椭圆柱投影”，故有以下几个概念：

名称	解释	图片解释
等角投影	投影面上两条方向所夹角度与地球面上对应的两条方向线所夹得角度相等 ，但是我们观察到，离最上面越远其面积变形的程度越大。要使面积不发生变化，则采用等面积投影。
横切椭圆柱投影	横轴投影又称“赤道投影”，投影面的轴与地球自转轴相垂直的一类投影，而由于地球并不是规则球形，因此将用于地球的圆柱投影叫做椭圆柱投影

通过上面的描述，可以看出高斯-克吕格投影的在角度上发生的变化较小，而在面积上发生的变化很大。我们把与椭圆柱相切的经线叫做中央经线，由于中央经线投影之后的长度比为1，且离中央经线两侧超过3度的区域，面积变形变形很大，为确保面积变形在可容许范围内，故只选取中央经线两侧3度范围的带（6度带）的经纬数据进行投影。

图1

图2

（2）高斯正算原理以及公式（LB----xy）

解算原理请参考：书本《大地测量学基础》， [链接: https://pan.baidu.com/s/1MF5PcwhNXr2PQHfS5ENIdw 密码: iw17

---

强烈推荐正反算公式参考网站：https://www.likecs.com/show-203333348.html#sc=502.3999938964844

本次程序主要利用的主公式，具体字母含义代码里面写的很清楚。

二、MATLAB高斯正算详解代码

1.源数据

注意：
1、若是待处理数据是txt文档，则采用file= importdata(filename) 函数读取，当数据的组织格式只含有BL数据，则直接使用点的索引data=file()即可读取。若其中含有字符串数据，则需要采用data=file.data()的点索引进行读取。

---

2、若待处理数据是xls或是其他格式文件，可以采用importdata或是xlsread函数进行读取。

2.源码

程序清单：
1、format long g： 临时关闭科学计数法的输出。
2、sind（）： 求取以度为单位的角度的正弦值，与sin（）有不同之处，注意区分。
3、writematrix（）： 将数据以xls格式进行输出，类似与xlswrite，当数据中有较多复杂的字符串，writematrix函数兼容性更高。
4、floor（）： 向下取整函数
5、本文主要采用3度带进行投影，若想采用6度带进行，则需进行以下更改：

N3=[];
ll=[];
Lo=[];
for i1=1:length(L)
    Lo_part=floor((L(i1)+6)/6)*6-3;
    Lo=[Lo;Lo_part]
    N3_part=floor((L(i1)+6)/6);
    ll_part=(L(i1)-Lo(i1))*3600;
    N3=[N3;N3_part];
    ll=[ll;ll_part];
end

6、为了方便初学者理解，我采用循环的形式进行输入和运算，若是要简化代码，可以取消循环采用 ’ . ’ +’ 运算符号 '对矩阵进行运算

%------------------------------------------------------------------------%
    %椭球基础参数
%------------------------------------------------------------------------%
clc;
clear all;
format long g;   
a=6378137;
f=1/298.257222101;
b=double(a*(1-f));
e1=double(((a^2-b^2)^(1/2))/a);
e2=double(((a^2-b^2)^(1/2))/b);
p=180/pi*3600;
%------------------------------------------------------------------------%
    %基本4m常量
%------------------------------------------------------------------------%
m0=a*(1-e1^2);
m2=(3/2)*(e1^2)*m0;
m4=(5/4)*(e1^2)*m2;
m6=(7/6)*(e1^2)*m4;
m8=(9/8)*(e1^2)*m6;
%------------------------------------------------------------------------%
    %基本4a常量
%------------------------------------------------------------------------%   
a0=sum([1,1/2,3/8,5/16,35/128].*([m0,m2,m4,m6,m8]));
a2=sum([1/2,1/2,15/35,7/16].*([m2,m4,m6,m8]));
a4=sum([1/8,3/16,7/32].*([m4,m6,m8]));
a6=sum([1/32,1/16].*([m6,m8]));
a8=(1/128)*m8;
%------------------------------------------------------------------------%
    %数据导入
%------------------------------------------------------------------------%
filename=importdata("WGS_84.txt");
L=filename.data(:,1);
B=filename.data(:,2);
%------------------------------------------------------------------------%

    %公式计算必要参数
    %中央子午线经度(3度带)
    %带号
    %经度差
%------------------------------------------------------------------------%
N3=[];
ll=[];
Lo=[];
for i1=1:length(L)
    Lo_part=floor((L(i1)+1.5)/3)*3;
    Lo=[Lo;Lo_part];
    N3_part=floor((L(i1)+1.5)/3);
    ll_part=(L(i1)-Lo(i1))*3600;
    N3=[N3;N3_part];
    ll=[ll;ll_part];
end
%------------------------------------------------------------------------%
    %子午线弧长
%------------------------------------------------------------------------%
   X=[];
for i2=1:length(B)
    sin2B=2*sind(B(i2))*cosd(B(i2));
    cos2B=2*((cosd(B(i2)))^2)-1;
    sin4B=2*sin2B*cos2B;
    cos4B=2*(cos2B^2)-1;
    sin6B=2*(sind(B(i2))*cos2B+cosd(B(i2))*sin2B)*(cosd(B(i2))*cos2B+sind(B(i2))*sin2B);
    sin8B=2*sin4B*cos4B;
    X_part=a0*(B(i2)*3600/p)-a2/2*sin2B+a4/4*sin4B-a6/6*sin6B+a8/8*sin8B;
    X=[X;X_part];
end
%------------------------------------------------------------------------%
     %卯酉圈曲率半径   
%------------------------------------------------------------------------%
   
N=[];
t=[];
nn=[];
for i3=1:length(B)
    N_part=a*(1-(e1^2)*(sind(B(i3)))^2)^(-1/2);
    t_part=tand(B(i3));
    nn_part=e2*cosd(B(i3));
    N=[N;N_part];
    t=[t;t_part];
    nn=[nn;nn_part];
end
%------------------------------------------------------------------------%
     %正算公式
%------------------------------------------------------------------------%
   x=[];
y_part1=[];
for i4=1:length(B)
    x_part=X(i4)+(N(i4)/(2*p^2))*sind(B(i4))*cosd(B(i4))*ll(i4)^2+(N(i4)/(24*p^4))*sind(B(i4))*(cosd(B(i4))^3)*(5-t(i4)^2+...
        9*nn(i4)^2+4*nn(i4)^4)*ll(i4)^4+(N(i4)/(720*p^6))*sind(B(i4))*(cosd(B(i4))^5)*(61-58*t(i4)^2+t(i4)^4)*ll(i4)^6;
    y_part=(N(i4)/p)*cosd(B(i4))*ll(i4)+(N(i4)/(6*p^3))*(cosd(B(i4))^3)*(1-t(i4)^2+nn(i4)^2)*ll(i4)^3+...
        (N(i4)/(120*p^5))*(cosd(B(i4))^5)*(5-18*t(i4)^2+t(i4)^4+14*nn(i4)^2-58*nn(i4)^2*t(i4)^2)*ll(i4)^5;
    x=[x;x_part];
    y_part1=[y_part1;y_part] ;
end
%------------------------------------------------------------------------%
    %为了方便测绘计算，所以要保证y值为正，即将坐标原点左移500km
    %以及确定带号的位置，因此在前面添加带号N3，此带号不参与后续的运算
    %y有带号
    %yy无带号
%------------------------------------------------------------------------%
y=[];
yy=[];
for i5=1:length(B)    
    y_part2=N3(i5)*1000000+y_part1(i5)+500000;   
    yy_part2=y_part1(i5)+500000;
    y=[y;y_part2];
    yy=[yy;yy_part2];
end
%------------------------------------------------------------------------%
    %数据输出
%------------------------------------------------------------------------%
B1={'B(纬度)'   'L(经度)'   '备注'};
null=zeros(4,1);
C1=[B L null];
E1={'纵坐标（x）'  '有带号横坐标（y）'  '无带号横坐标（yy）'};
F1=[x y yy];
B1_str=string(B1);
C1_str=string(C1);
E1_str=string(E1);
F1_str=string(F1);
endput1=[B1_str;C1_str;E1_str;F1_str];
AA1={'大地源坐标'
     '哈尔滨市'
     '西安市'
     '台北市'
     '乌鲁木齐市'
     '高斯正算之后的坐标'
     '哈尔滨市'
     '西安市'
     '台北市'
     '乌鲁木齐市'};
 AA1_str=string(AA1);
 endput=[AA1_str,endput1];
disp(endput);
%writematrix(endput,"C:\Users\稳魂\Desktop\endput.xls","WriteMode","overwritesheet");
%------------------------------------------------------------------------%
%------------------------------------------------------------------------%

3.运行结果

在测绘上为了方便后续计算而采用投影坐标系，其xy轴的方向刚好与笛卡尔数学坐标系xy的方向相反。

导出的xls表格文件：