在Java中将ECEF转换为LLA（纬度、经度、高度）

Question

我在网站上查阅了文章，但没有找到解决我的问题的方法......正如标题所说，我正在尝试将坐标从ECEF转换为LLA。

我使用的是这份文档：转换文章中的直接公式，而不是迭代公式，并使用ECEF2LLA网站进行结果比较。

我在Java中编写的代码如下：

public static final double a = 6378137;
public static final double f = 1/298.257223563;
public static final double b = a*(1-f);
public static final double e = Math.sqrt((Math.pow(a, 2)-Math.pow(b, 2))/Math.pow(a, 2));
public static final double e2 = Math.sqrt((Math.pow(a, 2)-Math.pow(b, 2))/Math.pow(b, 2));
public static double[] ecef2lla(double x , double y , double z){
    double[] lla = {0,0,0};
    double lon,lat,height,N;
    double p = Math.sqrt(Math.pow(x, 2)+Math.pow(y, 2));
    double theta = Math.atan((z*a)/(p*b));
    lon = Math.atan(y/x);
    lon = lon*180/Math.PI;

    lat = Math.atan(((z+Math.pow(e2, 2)*b*Math.pow(Math.sin(theta), 3))/((p-Math.pow(e,2)*a*Math.pow(Math.cos(theta), 3)))));
    lat = (lat*180)/Math.PI;

    N= a/(Math.sqrt(1-Math.pow(e*Math.sin(lat), 2)));
    height = (p/Math.cos(theta)) - N;
    lla[0] = lon;
    lla[1] = lat;
    lla[2] = height;
    return lla;

}

我得到了错误的高度数据。 我已经尝试过转换成弧度和角度等方法。

提前感谢您！

- dbkoren

---

请注意，这里展示的是鲍林方法，从技术上讲，它是迭代的。只是通常情况下，执行一次就足够准确了。 - undefined

Answer 1

我发现了这篇文章，准备在我的应用程序中使用“被接受的答案”，但我决定先运行一些测试来验证算法。我使用一个在线转换计算器（http://www.sysense.com/products/ecef_lla_converter/index.html）生成的样本数据。如下所示，输出结果并不太理想。

-----Test 1---------
Inputs:   -576793.17, -5376363.47, 3372298.51
Expected: 32.12345, -96.12345, 500.0
Actuals:  32.12306332822881, 83.87654999786477, 486.5472474489361
-----Test 2---------
Inputs:   2297292.91, 1016894.94, -5843939.62
Expected: -66.87654, 23.87654, 1000.0
Actuals:  -66.876230479461, 23.87653991401422, 959.6879360172898

接着，我使用以下帖子（https://gist.github.com/klucar/1536194）中的代码重新运行了相同的测试，并通过下面的输出结果得到了更好的结果。

-----Test 1---------
Inputs:   -576793.17, -5376363.47, 3372298.51
Expected: 32.12345, -96.12345, 500.0
Actuals:  32.12345004807767, -96.12345000213524, 499.997958839871
-----Test 2---------
Inputs:   2297292.91, 1016894.94, -5843939.62
Expected: -66.87654, 23.87654, 1000.0
Actuals:  -66.87654001741278, 23.87653991401422, 999.9983866894618

我没有花时间去找出“已接受的答案”中提供的解决方案中的错误，但我的建议是：使用这段代码...

/*
*
*  ECEF - Earth Centered Earth Fixed
*   
*  LLA - Lat Lon Alt
*
*  ported from matlab code at
*  https://gist.github.com/1536054
*     and
*  https://gist.github.com/1536056
*/

// WGS84 ellipsoid constants
private final double a = 6378137; // radius
private final double e = 8.1819190842622e-2;  // eccentricity

private final double asq = Math.pow(a,2);
private final double esq = Math.pow(e,2);

private double[] ecef2lla(double[] ecef){
  double x = ecef[0];
  double y = ecef[1];
  double z = ecef[2];

  double b = Math.sqrt( asq * (1-esq) );
  double bsq = Math.pow(b,2);
  double ep = Math.sqrt( (asq - bsq)/bsq);
  double p = Math.sqrt( Math.pow(x,2) + Math.pow(y,2) );
  double th = Math.atan2(a*z, b*p);

  double lon = Math.atan2(y,x);
  double lat = Math.atan2( (z + Math.pow(ep,2)*b*Math.pow(Math.sin(th),3) ), (p - esq*a*Math.pow(Math.cos(th),3)) );
  double N = a/( Math.sqrt(1-esq*Math.pow(Math.sin(lat),2)) );
  double alt = p / Math.cos(lat) - N;

  // mod lat to 0-2pi
  lon = lon % (2*Math.PI);

  // correction for altitude near poles left out.

  double[] ret = {lat, lon, alt};

  return ret;
}

- Kyle Melton

Answer 2

好的，我成功解决了这个问题。

问题是变量放错了位置，所以为了我们未来的方便，这里提供一个可行的JAVA实现：

public static final double a = 6378137;
public static final double f = 0.0034;
public static final double b = 6.3568e6;
public static final double e = Math.sqrt((Math.pow(a, 2) - Math.pow(b, 2)) / Math.pow(a, 2));
public static final double e2 = Math.sqrt((Math.pow(a, 2) - Math.pow(b, 2)) / Math.pow(b, 2));

public static double[] ecef2lla(double x, double y, double z) {

    double[] lla = { 0, 0, 0 };
    double lan, lon, height, N , theta, p;

    p = Math.sqrt(Math.pow(x, 2) + Math.pow(y, 2));

    theta = Math.atan((z * a) / (p * b));

    lon = Math.atan(y / x);

    lat = Math.atan(((z + Math.pow(e2, 2) * b * Math.pow(Math.sin(theta), 3)) / ((p - Math.pow(e, 2) * a * Math.pow(Math.cos(theta), 3)))));
    N = a / (Math.sqrt(1 - (Math.pow(e, 2) * Math.pow(Math.sin(lat), 2))));

    double m = (p / Math.cos(lat));
    height = m - N;


    lon = lon * 180 / Math.PI;
    lat = lat * 180 / Math.PI; 
    lla[0] = lat;
    lla[1] = lon;
    lla[2] = height;
    return lla;
}

注意：ECEF X Y Z的单位为米！

- dbkoren

---

为了获得所有可能经度的正确结果，我认为您需要使用atan2函数。 - Mageek

这里的f是什么？它没有被使用。 - Andrew Torr

在极地附近需要进行修正，其中高度 = dist(x,y,z) - sqrt(N^2 * cos(lat)^2 + (1-e2)^2 * N^2 * sin(lat)^2)。当然，你可以简化这个公式。 - undefined

Answer 3

如果您有兴趣使用已发布的Java库在笛卡尔坐标系和大地坐标系之间进行转换，您可以使用Java轨道传播库orekit。

虽然Orekit提供了一些高级功能，但它不是很容易入门。为了帮助您入门，我编写了一个示例代码来执行一些转换。我的代码是用Scala编写的，但它使用了Java库，所以仍然应该有所帮助。请参见此处的示例代码。

- Graham S