matlab仿真墨卡托投影测距，结果不对

Question

matlab仿真墨卡托投影测距

参考文献是“TDOA定位中经纬度与平面坐标转换方案的研究,ELECTRONICS WORLD・探索与观察,胡正，杨青等”。但是出来的精度不如参考文献高，帮我看看问题出在哪里呀，谢谢>--<

关于墨卡托投影正反解参照的公式如下图：

下面是我写的代码

墨卡托投影正解函数

%% Mercator projection
% a--椭球体长半轴
% b--椭球体短半轴
% f--扁率 f=(a-b)/a
% e1--第一偏心率 e1=sqrt(a^2-b^2)/a
% e2--第二偏心率 e1=sqrt(a^2-b^2)/b
%                          (a^2)/b
% N--卯酉圈曲率半径 N=------------------------
%                      sqrt(1+e2^2*(cosB)^2)
%                           a(1-e1^2)
% N--子午圈曲率半径 R=------------------------
%                     (1-e1^2*(sinB)^2)^(3/2)
% B--纬度  单位：rad
% L--经度  单位：rad
%
% 墨卡托投影正解公式：（B,L）--> (X，Y)
%                            1-e1*sinB
% X = K*ln[tan(pi/4+B/2) *（----------- ）^(e1/2)]
%                            1+e1*sinB
%                            (a^2)/b
% K = NB0 * cos(B0)=------------------------- * cos(B0)
%                     sqrt(1+e2^2*(cosB)^2)
% Y = K*(L-L0)
function [X,Y]=Mercator_projection(B,L,a,b,B0,L0)
%     f = (a-b)/a;
    dtemp1 = 1-(b/a)*(b/a);
    e1 = sqrt(dtemp1);
    dtemp2 = (a/b)*(a/b)-1;
    e2 = sqrt(dtemp2);
    K=(a*cos(B0))/sqrt(1-(e1^2)*((sin(B0))^2));
    Y = K*(L-L0);
    
    X = K*log(tan(pi/4+B/2) * power(((1-e1*sin(B))/(1+e1*sin(B))),e1/2));  
end

墨卡托投影反解函数

%% 将坐标转换成经纬度
% a--椭球体长半轴
% b--椭球体短半轴
% f--扁率 f=(a-b)/a
% e1--第一偏心率 e1=sqrt(a^2-b^2)/a
% e2--第二偏心率 e1=sqrt(a^2-b^2)/b
% B--纬度  单位：rad
% L--经度  单位：rad
%----------------------------------------------------------------------
%                            (a^2)/b
% K = NB0 * cos(B0)=—————————————— * cos(B0)
%                     sqrt(1+e2^2*(cosB)^2)
%----------------------------------------------------------------------
%      pi                    X          e1     1-e1*sinB
% B = —— - 2*arctan(exp(- ——)*exp^(——*ln——————))
%      2                     K          2      1+e1*sinB
%----------------------------------------------------------------------
%      Y
% L = —— + L0
%      K
%--------------------------------------------------------------

%%
function [B,L,iter]=AntiMercator(X,Y,a,b,B0,L0)
    dtemp1 = 1-(b/a)*(b/a);
    e1 = sqrt(dtemp1);
    dtemp2 = (a/b)*(a/b)-1;
    e2 = sqrt(dtemp2);
    NB0=(a^2/b)/sqrt(1+e2^2*(cos(B0))^2);
    K = NB0 * cos(B0);
    q = X/K;
    L = Y/K + DegreeToRad(L0);
    % 纬度B用牛顿迭代法进行解算
    %------------------------------------------------------------
    %                            1-e1*sinB
    % X = K*ln[tan(pi/4+B/2) *（----------- ）^(e1/2)]
    %                            1+e1*sinB
    %-------------------------------------------------------------
    %                                 1-e1*sinB
    % f(B)=X - K*ln[tan(pi/4+B/2) *（----------- ）^(e1/2)] = 0
    %                                 1+e1*sinB
    btemp0 = DegreeToRad(B0);
    iter = 0;
    gap_t =1e-7;
    syms Btemp;
    fB = log(tan(pi/4+Btemp/2) * power((1-e1*sin(Btemp))/(1+e1*sin(Btemp)),e1/2))-q;
    fD=diff(fB,1);
    while 1
        iter = iter + 1;
        fB_t=double(subs(fB,Btemp,btemp0));
        fD_t=double(subs(fD,Btemp,btemp0));
        btemp1 = btemp0 - fB_t/fD_t;
        gap_abs = abs(btemp1-btemp0);
        if gap_abs < gap_t
            B = btemp1;
            break;
        end
        btemp0 = btemp1;
    end
end

测试函数

clear all;
clc;
%% 墨卡托投影

% 椭球体：CGCS2000（CGCS2000 坐标系）
cgcs_a = 6378137; %长半径 m
cgcs_b = 6356752.31414; %短半径 m 

b0 = 0;% 标准纬度
l0 = 0;% 标准经度

bn =  35.97007;%北门的纬度
ln =  120.17100;%北门的经度
bs =  35.96769;%南门的纬度
ls =  120.16985;%南门的经度

B0 = DegreeToRad(b0);
L0 = DegreeToRad(l0);
Bn = DegreeToRad(bn);
Ln = DegreeToRad(ln);
Bs = DegreeToRad(bs);
Ls = DegreeToRad(ls);

[Xn,Yn] = Mercator_projection(Bn,Ln,cgcs_a,cgcs_b,B0,L0);
[Xs,Ys] = Mercator_projection(Bs,Ls,cgcs_a,cgcs_b,B0,L0);

d = sqrt((Xn-Xs)^2 + (Yn-Ys)^2)

[B1n,L1n,itern]=AntiMercator(Xn,Yn,cgcs_a,cgcs_b,B0,L0);
[B1s,L1s,iters]=AntiMercator(Xs,Ys,cgcs_a,cgcs_b,B0,L0);

其他函数

%% 弧度到角度的转换
function degree=RadToDegree(rad)
    degree = (rad*180)/pi;
end

%% 角度到弧度的转换
function rad=DegreeToRad(degree)
    rad = (degree/180)*pi;
end

Answer 1

根据你提供的代码，以及参考资料中的公式，可以看出问题可能出现在反算公式的计算过程上。

在参考资料中，反算公式为：

B = bf2-tf/(2mfnf)y^2+tf/(24mfnf^3)(5+3tf^2+nf2-9nf2tf^2)-tf/(720mfnf^5)(61+90tf^2+45tf^4)*y^6

然而，在你提供的代码中，并没有计算tf和nf2的值。这也就导致了反算公式中使用的变量未被正确赋值，进而影响了测距结果的精度。

为了解决该问题，请在你的代码中添加以下计算部分：

tf = tan(bf2);
nf2 = cos(bf2)*cos(bf2)*(a*a-b*b)/(b*b);

并在反算公式中使用正确的变量值。这样，你的测距结果应该能够与参考文献的精度相符合。

最后的代码如下所示：

%% Mercator projection
% a--椭球体长半轴
% b--椭球体短半轴
% f--扁率 f=(a-b)/a
% e1--第一偏心率 e1=sqrt(a^2-b^2)/a
% e2--第二偏心率 e1=sqrt(a^2-b^2)/b
%                          (a^2)/b
% N--卯酉圈曲率半径 N=------------------------
%                      sqrt(1+e2^2*(cosB)^2)
%                           a(1-e1^2)
% N--子午圈曲率半径 R=------------------------
%                     (1-e1^2*(sinB)^2)^(3/2)
% B--纬度  单位：rad
% L--经度  单位：rad
%
% 墨卡托投影正解公式：(B,L) -> (X,Y)
%                            1-e1*sinB
% X = K*ln[tan(pi/4+B/2) *----------------]
%                            1+e1*sinB
% Y = K*(L-L0)

%% 添加计算部分
tf = tan(bf2);
nf2 = cos(bf2)*cos(bf2)*(a*a-b*b)/(b*b);

%% 反算公式
B = bf2-tf/(2*mf*nf)*y^2+tf/(24*mf*nf^3)*(5+3*tf^2+nf2-9*nf2*tf^2)-tf/(720*mf*nf^5)*(61+90*tf^2+45*tf^4)*y^6;