大地坐标与空间直角坐标转换及MATLAB代码

(1)大地坐标与空间直角坐标的转换公式

两者之间的转换公式如下所示
$\left\{ \begin{array}{l} X = (N + H)\cos B\cos L\\ Y = (N + H)\cos B\sin L\\ Z = [N(1 - {e^2}) + H]\sin B \end{array} \right.\\$⎩⎨⎧​X=(N+H)cosBcosLY=(N+H)cosBsinLZ=[N(1−e2)+H]sinB​

$\left\{ \begin{array}{l} L = \arctan (Y/X)\\ B = \arctan \{ Z(N + H)/[\sqrt {{X^2} + {Y^2}} (N(1 - {e^2}) + H)]\} \\ H = Z/\sin B - N(1 - {e^2}) \end{array} \right.$⎩⎨⎧​L=arctan(Y/X)B=arctan{Z(N+H)/[X2+Y2​(N(1−e2)+H)]}H=Z/sinB−N(1−e2)​

式中：$(X,Y,Z)$(X,Y,Z)和$(B,L,H)$(B,L,H)分别为对应点的空间直角坐标和大地坐标，$N = a/\sqrt {(1 - {e^2}{{\sin }^2}B)}$N=a/(1−e2sin2B)​为卯酉圈曲率半径 ，$a$a 、$e$e分别为大地坐标系对应椭球的长半轴和第一偏心率。

(2)大地坐标转换至空间直角坐标

%   intput：
%           (Lon，Lat，H)------------经纬度和大地高，可以为列向量
%
%           ReferenceFrame----------参考框架
%   output：
%           (X,Y,Z)-----------------空间直角坐标
%
function [X,Y,Z]=blh2xyz(Lon,Lat,H,ReferenceFrame)

if nargin<4
    ReferenceFrame='GRS80';
end
switch ReferenceFrame
    %GRS80椭球参数，根据IERS2010得到
    case 'GRS80'
        a=6378137;%长半轴
        f=1/298.257222101;%扁率
    %WGS84椭球参数
    case 'WGS84'
        a=6378137;%长半轴
        f=1/298.257223563;%扁率
end

Lon=Lon*pi/180;
Lat=Lat*pi/180;
b=(1-f)*a;%短半轴
e2=1-b^2/a^2;%偏心率的平方
N=a./sqrt(1-e2*sin(Lat).^2);%卯酉圈曲率半径

%空间直角坐标
X=(N+H).*cos(Lat).*cos(Lon);
Y=(N+H).*cos(Lat).*sin(Lon);
Z=(N*(1-e2)+H).*sin(Lat);

(3)空间直角坐标转换至大地坐标

%   input：
%           (X,Y,Z)-----------------空间直角坐标,可以为列向量
%           ReferenceFrame----------参考框架
%   output：
%           Lon---------------------经度
%           Lat---------------------纬度
%           H-----------------------大地高

function [Lon,Lat,H]=xyz2blh(X,Y,Z,ReferenceFrame)
if nargin<4
    ReferenceFrame='GRS80';
end
switch ReferenceFrame
    %GRS80椭球参数，根据IERS2010得到
    case 'GRS80'
        a=6378137;%长半轴
        f=1/298.257222101;%扁率
        %WGS84椭球参数
    case 'WGS84'
        a=6378137;%长半轴
        f=1/298.257223563;%扁率
        
end

%短半轴、偏心率的平方
b=(1-f)*a;
e2=1-b^2/a^2;

%计算大地经度
Lon=atan(Y./X);%X等于零的情况包括在内
isSecondQuadrant=X<0&Y>0;
isThirdQuadrant=X<0&Y<0;
Lon(isSecondQuadrant)=Lon(isSecondQuadrant)+pi;
Lon(isThirdQuadrant)=Lon(isThirdQuadrant)-pi;

%计算大地纬度和大地高
Lat=atan(Z./(sqrt(X.^2+Y.^2)*(1-e2)));%初值
while(1)%不动点迭代
    N=a./sqrt(1-e2*sin(Lat).^2);%卯酉圈曲率半径
    H=Z./sin(Lat)-N*(1-e2);
    temp=atan(Z.*(N+H)./sqrt(X.^2+Y.^2)./(N*(1-e2)+H));
    notAccurate=abs(temp-Lat)>=1e-10;
    Lat(notAccurate)=temp(notAccurate);
    if(sum(sum(notAccurate))==0)
        break;
    end
end
H=Z./sin(Lat)-N*(1-e2);
Lon=Lon*180/pi;
Lat=Lat*180/pi;