coordinates - 如何在 WGS84 椭球上插值 ECEF 坐标

Question

是否有直接的方法（不涉及将坐标转换为纬度/经度）在 2 个 ECEF 坐标（xyz）之间进行插值，以便插值点位于 WGS84 椭球上。原始的 2 个点是根据大地坐标计算的。

在球体上插值似乎很明显，但我似乎无法得出椭球的解决方案。

先感谢您。

Answer 1

假设您得到 2 分p0(x,y,z)，p1(x,y,z)并且想要在两者之间插入一些p(t)位置。t=<0.0,1.0>

你可以：

1. 将您的椭球重新缩放为球体

   就像这样：

   const double mz=6378137.00000/6356752.31414; // [m] equatoreal/polar radius of Earth
   p0.z*=mz;
   p1.z*=mz;
   

   现在你得到了指向球形地球模型的笛卡尔坐标。

2. 插

   简单的线性插值就可以了

   p(t) = p0+(p1-p0)*t
   

   但是粗略的你还需要归一化为地球曲率，所以：

   r0 = |p0|
   r1 = |p1|
   
   p(t) = p0+(p1-p0)*t
   r(t) = r0+(r1-r0)*t
   
   p(t)*=r/|p(t)|
   

   其中|p0|表示向量的长度p0。

3. 重新缩放回椭球体

   通过除以相同的值

   p(t).z/=mz
   

这既简单又便宜，但内插路径不会具有线性时间尺度。

这里的 C++ 示例：

void XYZ_interpolate(double *pt,double *p0,double *p1,double t)
    {
    const double  mz=6378137.00000/6356752.31414;
    const double _mz=6356752.31414/6378137.00000;
    double p[3],r,r0,r1;
    // compute spherical radiuses of input points
    r0=sqrt((p0[0]*p0[0])+(p0[1]*p0[1])+(p0[2]*p0[2]*mz*mz));
    r1=sqrt((p1[0]*p1[0])+(p1[1]*p1[1])+(p1[2]*p1[2]*mz*mz));
    // linear interpolation
    r   = r0   +(r1   -r0   )*t;
    p[0]= p0[0]+(p1[0]-p0[0])*t;
    p[1]= p0[1]+(p1[1]-p0[1])*t;
    p[2]=(p0[2]+(p1[2]-p0[2])*t)*mz;
    // correct radius and rescale back
    r/=sqrt((p[0]*p[0])+(p[1]*p[1])+(p[2]*p[2]));
    pt[0]=p[0]*r;
    pt[1]=p[1]*r;
    pt[2]=p[2]*r*_mz;
    }

并预览：

黄色方块是使用的p0,p1笛卡尔坐标，白色曲线是插值路径，其中t=<0.0,1.0>...