超腾开源

大地坐标系转换工具 - 北京54西安80转WGS84四参数七参数计算器

大地坐标系转换工具

北京54/西安80 → WGS84/CGCS2000，支持四参数和七参数转换

使用说明
支持四参数和七参数进行坐标系转换，需要输入实测参数。

完整转换（投影坐标→WGS84）

计算原理

完整转换流程

从投影坐标（如北京54）转换为WGS84经纬度需要以下步骤：

  1. 高斯投影反算：投影坐标(X,Y) → 大地坐标(B,L)
  2. 大地转空间：大地坐标 → 空间直角坐标(X,Y,Z)
  3. 七参数转换：源椭球空间坐标 → 目标椭球空间坐标
  4. 空间转大地：空间直角坐标 → WGS84大地坐标

四参数（平面）

        X₂ = ΔX + (1+K×10⁻⁶)(X₁cosθ - Y₁sinθ)
        Y₂ = ΔY + (1+K×10⁻⁶)(X₁sinθ + Y₁cosθ)

七参数（布尔莎模型）

        X₂ = ΔX + (1+K×10⁻⁶)(X₁ - Rz×Y₁ + Ry×Z₁)
        Y₂ = ΔY + (1+K×10⁻⁶)(Rz×X₁ + Y₁ - Rx×Z₁)
        Z₂ = ΔZ + (1+K×10⁻⁶)(-Ry×X₁ + Rx×Y₁ + Z₁)

参数说明

  • ΔX,ΔY,ΔZ：三个平移参数（米）
  • Rx,Ry,Rz：三个旋转参数（秒）
  • K：尺度因子（ppm）

代码实现

JavaScript

// 完整坐标转换：投影坐标 → WGS84经纬度
const ellipsoids = {
  BJ54: { a: 6378245, f: 1 / 298.3 },
  XA80: { a: 6378140, f: 1 / 298.257 },
  WGS84: { a: 6378137, f: 1 / 298.257223563 }
};

// 高斯投影反算
function gaussInverse(x, y, L0, ellipsoid) {
  const { a, f } = ellipsoid;
  const e2 = 2 * f - f * f;
  const e12 = e2 / (1 - e2);
  // 底点纬度 Bf（子午线弧长反算级数展开）
  const m0 = 1 - e2 / 4 - 3 * e2 ** 2 / 64 - 5 * e2 ** 3 / 256;
  const mu = x / (a * m0);
  const e1 = (1 - Math.sqrt(1 - e2)) / (1 + Math.sqrt(1 - e2));
  const Bf = mu
    + (3 * e1 / 2 - 27 * e1 ** 3 / 32) * Math.sin(2 * mu)
    + (21 * e1 ** 2 / 16 - 55 * e1 ** 4 / 32) * Math.sin(4 * mu)
    + (151 * e1 ** 3 / 96) * Math.sin(6 * mu)
    + (1097 * e1 ** 4 / 512) * Math.sin(8 * mu);
  const Wf = Math.sqrt(1 - e2 * Math.sin(Bf) ** 2);
  const Nf = a / Wf;
  const Mf = a * (1 - e2) / Wf ** 3;
  const tf = Math.tan(Bf);
  const eta2 = e12 * Math.cos(Bf) ** 2;
  const B = Bf - (Nf * tf / Mf) *
    (y ** 2 / (2 * Nf ** 2) - (5 + 3 * tf ** 2 + eta2 - 9 * eta2 * tf ** 2) * y ** 4 / (24 * Nf ** 4));
  const L = L0 + y / (Nf * Math.cos(Bf)) - (1 + 2 * tf ** 2 + eta2) * y ** 3 / (6 * Nf ** 3 * Math.cos(Bf));
  return { B: B * 180 / Math.PI, L: L * 180 / Math.PI };
}

// 大地坐标→空间直角坐标
function blhToXyz(B, L, H, ellipsoid) {
  const { a, f } = ellipsoid;
  const e2 = 2 * f - f * f;
  const N = a / Math.sqrt(1 - e2 * Math.sin(B * Math.PI / 180) ** 2);
  const x = (N + H) * Math.cos(B * Math.PI / 180) * Math.cos(L * Math.PI / 180);
  const y = (N + H) * Math.cos(B * Math.PI / 180) * Math.sin(L * Math.PI / 180);
  const z = (N * (1 - e2) + H) * Math.sin(B * Math.PI / 180);
  return { x, y, z };
}

// 七参数转换
function sevenParam(x, y, z, dx, dy, dz, rx, ry, rz, k) {
  const rxr = rx * Math.PI / 648000;
  const ryr = ry * Math.PI / 648000;
  const rzr = rz * Math.PI / 648000;
  const m = 1 + k * 1e-6;
  return {
    x: dx + m * (x - rzr * y + ryr * z),
    y: dy + m * (rzr * x + y - rxr * z),
    z: dz + m * (-ryr * x + rxr * y + z)
  };
}

// 空间直角坐标→大地坐标
function xyzToBlh(x, y, z, ellipsoid) {
  const { a, f } = ellipsoid;
  const e2 = 2 * f - f * f;
  const L = Math.atan2(y, x);
  let B = Math.atan2(z, Math.sqrt(x ** 2 + y ** 2));
  for (let i = 0; i < 5; i++) {
    const N = a / Math.sqrt(1 - e2 * Math.sin(B) ** 2);
    const H = Math.sqrt(x ** 2 + y ** 2) / Math.cos(B) - N;
    B = Math.atan2(z, Math.sqrt(x ** 2 + y ** 2) * (1 - e2 * N / (N + H)));
  }
  return { B: B * 180 / Math.PI, L: L * 180 / Math.PI };
}

// 使用示例：北京54投影坐标转WGS84
const projX = 4426641.2, projY = -47244.5, L0 = 117 * Math.PI / 180;
const bl = gaussInverse(projX, projY, L0, ellipsoids.BJ54);
const xyz1 = blhToXyz(bl.B, bl.L, 0, ellipsoids.BJ54);
const xyz2 = sevenParam(xyz1.x, xyz1.y, xyz1.z, -29.7, -130.6, -95.8, 0.0025, 0.0033, 0, 2.4);
const wgs84 = xyzToBlh(xyz2.x, xyz2.y, xyz2.z, ellipsoids.WGS84);
console.log('WGS84:', wgs84);

Java

public class DatumTransform {
    static class Ellipsoid {
        double a, f;
        Ellipsoid(double a, double f) { this.a = a; this.f = f; }
    }

    static Ellipsoid BJ54 = new Ellipsoid(6378245, 1.0/298.3);
    static Ellipsoid WGS84 = new Ellipsoid(6378137, 1.0/298.257223563);

    // 高斯投影反算
    static double[] gaussInverse(double x, double y, double L0, Ellipsoid e) {
        double e2 = 2 * e.f - e.f * e.f;
        double e12 = e2 / (1 - e2);
        // 底点纬度 Bf（子午线弧长反算级数展开）
        double m0 = 1 - e2 / 4 - 3 * Math.pow(e2, 2) / 64 - 5 * Math.pow(e2, 3) / 256;
        double mu = x / (e.a * m0);
        double e1 = (1 - Math.sqrt(1 - e2)) / (1 + Math.sqrt(1 - e2));
        double Bf = mu
            + (3 * e1 / 2 - 27 * Math.pow(e1, 3) / 32) * Math.sin(2 * mu)
            + (21 * e1 * e1 / 16 - 55 * Math.pow(e1, 4) / 32) * Math.sin(4 * mu)
            + (151 * Math.pow(e1, 3) / 96) * Math.sin(6 * mu)
            + (1097 * Math.pow(e1, 4) / 512) * Math.sin(8 * mu);
        double Wf = Math.sqrt(1 - e2 * Math.sin(Bf) * Math.sin(Bf));
        double Nf = e.a / Wf;
        double Mf = e.a * (1 - e2) / Math.pow(Wf, 3);
        double tf = Math.tan(Bf);
        double eta2 = e12 * Math.cos(Bf) * Math.cos(Bf);
        double B = Bf - (Nf * tf / Mf) *
            (y * y / (2 * Nf * Nf) - (5 + 3 * tf * tf + eta2 - 9 * eta2 * tf * tf) * Math.pow(y, 4) / (24 * Math.pow(Nf, 4)));
        double L = L0 + y / (Nf * Math.cos(Bf)) - (1 + 2 * tf * tf + eta2) * Math.pow(y, 3) / (6 * Math.pow(Nf, 3) * Math.cos(Bf));
        return new double[] {Math.toDegrees(B), Math.toDegrees(L)};
    }

    // 大地坐标→空间直角坐标
    static double[] blhToXyz(double B, double L, double H, Ellipsoid e) {
        double e2 = 2 * e.f - e.f * e.f;
        double Br = Math.toRadians(B), Lr = Math.toRadians(L);
        double N = e.a / Math.sqrt(1 - e2 * Math.sin(Br) * Math.sin(Br));
        return new double[] {
            (N + H) * Math.cos(Br) * Math.cos(Lr),
            (N + H) * Math.cos(Br) * Math.sin(Lr),
            (N * (1 - e2) + H) * Math.sin(Br)
        };
    }

    // 七参数转换
    static double[] sevenParam(double x, double y, double z,
        double dx, double dy, double dz, double rx, double ry, double rz, double k) {
        double rxr = rx * Math.PI / 648000;
        double ryr = ry * Math.PI / 648000;
        double rzr = rz * Math.PI / 648000;
        double m = 1 + k * 1e-6;
        return new double[] {
            dx + m * (x - rzr * y + ryr * z),
            dy + m * (rzr * x + y - rxr * z),
            dz + m * (-ryr * x + rxr * y + z)
        };
    }

    // 空间直角坐标→大地坐标
    static double[] xyzToBlh(double x, double y, double z, Ellipsoid e) {
        double e2 = 2 * e.f - e.f * e.f;
        double L = Math.atan2(y, x);
        double B = Math.atan2(z, Math.sqrt(x * x + y * y));
        for (int i = 0; i < 5; i++) {
            double N = e.a / Math.sqrt(1 - e2 * Math.sin(B) * Math.sin(B));
            double H = Math.sqrt(x * x + y * y) / Math.cos(B) - N;
            B = Math.atan2(z, Math.sqrt(x * x + y * y) * (1 - e2 * N / (N + H)));
        }
        return new double[] {Math.toDegrees(B), Math.toDegrees(L)};
    }

    public static void main(String[] args) {
        double[] bl = gaussInverse(4426641.2, -47244.5, Math.toRadians(117), BJ54);
        double[] xyz1 = blhToXyz(bl[0], bl[1], 0, BJ54);
        double[] xyz2 = sevenParam(xyz1[0], xyz1[1], xyz1[2], -29.7, -130.6, -95.8, 0.0025, 0.0033, 0, 2.4);
        double[] wgs84 = xyzToBlh(xyz2[0], xyz2[1], xyz2[2], WGS84);
        System.out.printf("WGS84: B=%.8f°, L=%.8f°%n", wgs84[0], wgs84[1]);
    }
}

Python

import math

class Ellipsoid:
    def __init__(self, a, f):
        self.a = a
        self.f = f

BJ54 = Ellipsoid(6378245, 1/298.3)
WGS84 = Ellipsoid(6378137, 1/298.257223563)

def gauss_inverse(x, y, L0, ellipsoid):
    """高斯投影反算"""
    e2 = 2 * ellipsoid.f - ellipsoid.f ** 2
    e12 = e2 / (1 - e2)
    # 底点纬度 Bf（子午线弧长反算级数展开）
    m0 = 1 - e2 / 4 - 3 * e2 ** 2 / 64 - 5 * e2 ** 3 / 256
    mu = x / (ellipsoid.a * m0)
    e1 = (1 - math.sqrt(1 - e2)) / (1 + math.sqrt(1 - e2))
    Bf = (mu
        + (3 * e1 / 2 - 27 * e1 ** 3 / 32) * math.sin(2 * mu)
        + (21 * e1 ** 2 / 16 - 55 * e1 ** 4 / 32) * math.sin(4 * mu)
        + (151 * e1 ** 3 / 96) * math.sin(6 * mu)
        + (1097 * e1 ** 4 / 512) * math.sin(8 * mu))
    Wf = math.sqrt(1 - e2 * math.sin(Bf) ** 2)
    Nf = ellipsoid.a / Wf
    Mf = ellipsoid.a * (1 - e2) / Wf ** 3
    tf = math.tan(Bf)
    eta2 = e12 * math.cos(Bf) ** 2
    B = Bf - (Nf * tf / Mf) * (
        y ** 2 / (2 * Nf ** 2) - (5 + 3 * tf ** 2 + eta2 - 9 * eta2 * tf ** 2) * y ** 4 / (24 * Nf ** 4))
    L = L0 + y / (Nf * math.cos(Bf)) - (1 + 2 * tf ** 2 + eta2) * y ** 3 / (6 * Nf ** 3 * math.cos(Bf))
    return math.degrees(B), math.degrees(L)

def blh_to_xyz(B, L, H, ellipsoid):
    """大地坐标→空间直角坐标"""
    e2 = 2 * ellipsoid.f - ellipsoid.f ** 2
    N = ellipsoid.a / math.sqrt(1 - e2 * math.sin(math.radians(B)) ** 2)
    x = (N + H) * math.cos(math.radians(B)) * math.cos(math.radians(L))
    y = (N + H) * math.cos(math.radians(B)) * math.sin(math.radians(L))
    z = (N * (1 - e2) + H) * math.sin(math.radians(B))
    return x, y, z

def seven_param(x, y, z, dx, dy, dz, rx, ry, rz, k):
    """七参数转换"""
    rxr = rx * math.pi / 648000
    ryr = ry * math.pi / 648000
    rzr = rz * math.pi / 648000
    m = 1 + k * 1e-6
    return (
        dx + m * (x - rzr * y + ryr * z),
        dy + m * (rzr * x + y - rxr * z),
        dz + m * (-ryr * x + rxr * y + z)
    )

def xyz_to_blh(x, y, z, ellipsoid):
    """空间直角坐标→大地坐标"""
    e2 = 2 * ellipsoid.f - ellipsoid.f ** 2
    L = math.atan2(y, x)
    B = math.atan2(z, math.sqrt(x ** 2 + y ** 2))
    for _ in range(5):
        N = ellipsoid.a / math.sqrt(1 - e2 * math.sin(B) ** 2)
        H = math.sqrt(x ** 2 + y ** 2) / math.cos(B) - N
        B = math.atan2(z, math.sqrt(x ** 2 + y ** 2) * (1 - e2 * N / (N + H)))
    return math.degrees(B), math.degrees(L)

# 使用示例
B, L = gauss_inverse(4426641.2, -47244.5, math.radians(117), BJ54)
xyz1 = blh_to_xyz(B, L, 0, BJ54)
xyz2 = seven_param(*xyz1, -29.7, -130.6, -95.8, 0.0025, 0.0033, 0, 2.4)
wgs84 = xyz_to_blh(*xyz2, WGS84)
print(f"WGS84: B={wgs84[0]:.8f}°, L={wgs84[1]:.8f}°")

评论 (0)

请 登录 后发表评论

暂无评论，快来发表第一条评论吧！