matlab仿真墨卡托投影测距，结果不对
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
根据你提供的代码，以及参考资料中的公式，可以看出问题可能出现在反算公式的计算过程上。
在参考资料中，反算公式为：
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;