七参数尺度因子没除 1e6，20 公里外偏了 24 公里

现象

测区原点附近（1 公里内）转换后的点与已知点相差不到 1 米，看着正常。 到了测区边缘 20 公里处的点，转换结果与已知控制点相差约 24 公里， 而且偏差方向沿原点向外的径向，所有外围点都被「推远」了。

把尺度参数从 -1.2 改成 -1.2e-6 后，全线恢复正常。

原因

布尔莎七参数的尺度因子 m 单位是 ppm（百万分之一）， 在公式里必须先除以 1e6：

X' = ΔX + (1 + m) · (X + ωz·Y - ωy·Z)
Y' = ΔY + (1 + m) · (Y - ωz·X + ωx·Z)
Z' = ΔZ + (1 + m) · (Z + ωy·X - ωx·Y)
其中 m 需代入 m_ppm × 1e-6

直接代入 -1.2 意味着「尺度缩小到 -0.2 倍」，即距离变成原来的 1.2 倍反号， 径向误差 = 距离 × |m|。20 公里 × 1.2 = 24 公里，与观测吻合。 原点附近距离接近 0，误差自然看不出来，这是最容易漏检的地方。

排查过程

1. 取三个不同距离的检核点（1 km、10 km、20 km），分别算转换残差。
2. 若残差与距离近似成正比，就是尺度或旋转参数的单位问题。
3. 把尺度参数按 ppm 与无量纲两种方式各算一次，比较残差量级。
4. 旋转参数同样检查：ω 若以秒给出，需除以 206265 换成弧度。

// 残差随距离线性增长，基本可锁定尺度参数
const pts = [
  { d: 1000,  err: 1180 },
  { d: 10000, err: 12010 },
  { d: 20000, err: 24050 }
];
console.log(pts.map(p => (p.err / p.d).toFixed(4)));
// [1.1800, 1.2010, 1.2025]  -> 比例约 1.2，等于 m 的值

解决方案

封装参数时统一在入口做单位换算，之后内部只使用国际单位。

// 参数原始值（来自测绘成果表）
const P = {
  dx: 12.345, dy: -8.210, dz: 4.670,       // 米
  rx: 0.0000012, ry: -0.0000021, rz: 0.0000009, // 弧度
  m: -1.2                                   // ppm
};

function toCartesian(x, y, z) {
  const m = P.m * 1e-6;                     // 关键：ppm -> 无量纲
  const k = 1 + m;
  return [
    P.dx + k * (x - P.rz * y + P.ry * z),
    P.dy + k * (P.rz * x + y - P.rx * z),
    P.dz + k * (-P.ry * x + P.rx * y + z)
  ];
}

用 proj4js 时不要手写公式，把它写进 +towgs84=， 注意 towgs84 不接受 ppm，只能填三个平移和三个旋转（秒），尺度需另作处理。

proj4.defs('LOCAL_2000',
  '+proj=tmerc +lat_0=0 +lon_0=108 +k=1 +x_0=500000 +y_0=0' +
  ' +ellps=GRS80 +towgs84=12.345,-8.210,4.670,0.247,-0.433,0.186,-1.2' +
  ' +units=m +no_defs');

验证

重新计算 20 公里处的检核点，残差应从 24 公里降到 5 厘米以内； 三个距离档位的残差不应再随距离增长。

const r = [1000, 10000, 20000].map(d => {
  const p = toCartesian(d, 0, 0);
  return Math.hypot(p[0] - d - P.dx, p[1] - P.dy);
});
console.assert(r[2] < 0.05, '远端残差仍过大：' + r[2]);

适用版本