站心直角坐标系转经纬高

序：网上很多算法，基本都是错的，还得去找论文才行。原理就不讲了，我会贴出论文引用。

1.转换过程

(1)确定站心直角坐标系的原点(X0,Y0,Z0)=(0,0,0)的经纬高(B0,L0,H0)作为基准点。

(2)将基准点(B0,L0,H0)转成地心坐标。

(3)将站心直角坐标系的目标点(X1,Y1,Z1)转成地心坐标(eX1,eY1,eZ1)。

(4)将地心坐标(eX1,eY1,eZ1)转成经纬高(B1,L1,H1)。

2.C#代码

    /// <summary>
        /// 使用WGS84椭球体
    /// </summary>
    public class CoorTransform
    {
        #region C#实现方法/三角函数及迭代法
        /// <summary>
        /// 地心坐标转经纬高
        /// </summary>
        /// <param name="X"></param>
        /// <param name="Y"></param>
        /// <param name="Z"></param>
        /// <param name="aAxis"></param>
        /// <param name="bAxis"></param>
        /// <returns></returns>
        public static Vector3D XYZtoBLH(double X, double Y, double Z, double aAxis = 6378137,             double bAxis = 6356752.3142)
		{
            double targetH, targetB, targetL;
            double e1 = (Math.Pow(aAxis, 2) - Math.Pow(bAxis, 2)) / Math.Pow(aAxis, 2);
            double e2 = (Math.Pow(aAxis, 2) - Math.Pow(bAxis, 2)) / Math.Pow(bAxis, 2);

            double S = Math.Sqrt(Math.Pow(X, 2) + Math.Pow(Y, 2));
            double cosL = X / S;
            double B = 0;
            double L = 0;

            L = Math.Acos(cosL);
            L = Math.Abs(L);

            double tanB = Z / S;
            B = Math.Atan(tanB);
            double c = aAxis * aAxis / bAxis;
            double preB0 = 0.0;
            double ll = 0.0;
            double N = 0.0;
            //迭代计算纬度
            do
            {
                preB0 = B;
                ll = Math.Pow(Math.Cos(B), 2) * e2;
                N = c / Math.Sqrt(1 + ll);

                tanB = (Z + N * e1 * Math.Sin(B)) / S;
                B = Math.Atan(tanB);
            } while (Math.Abs(preB0 - B) >= 0.0000000001);

            ll = Math.Pow(Math.Cos(B), 2) * e2;
            N = c / Math.Sqrt(1 + ll);

            targetH = S / Math.Cos(B) - N;// * (1 - e1);
            targetB = B * 180 / Math.PI;
            targetL = L * 180 / Math.PI;
            return new Vector3D(targetB, targetL, targetH);
        }
        /// <summary>
        /// 经纬高转地心坐标
        /// </summary>
        /// <param name="B"></param>
        /// <param name="L"></param>
        /// <param name="H"></param>
        /// <param name="aAxis"></param>
        /// <param name="bAxis"></param>
        /// <returns></returns>
		public static Vector3D BLHtoXYZ(double B, double L, double H, double aAxis = 6378137, double bAxis = 6356752.3142)
		{
			double targetX, targetY, targetZ;
			double dblD2R = Math.PI / 180.0;
            double f = 1 / 298.257223563;
            double e2 = 2 * f - f * f;
            //double e1 = Math.Sqrt(Math.Pow(aAxis, 2) - Math.Pow(bAxis, 2)) / Math.Pow(aAxis, 2);

            double W = Math.Sqrt(1 - e2 * Math.Sin(B * dblD2R) * Math.Sin(B * dblD2R));
            double N = aAxis / W;
			//double N = aAxis / Math.Sqrt(1.0 - Math.Pow(e1, 2) * Math.Pow(Math.Sin(B * dblD2R), 2));
			targetX = (N + H) * Math.Cos(B * dblD2R) * Math.Cos(L * dblD2R);
			targetY = (N + H) * Math.Cos(B * dblD2R) * Math.Sin(L * dblD2R);
			targetZ = (N * (1.0 - e2) + H) * Math.Sin(B * dblD2R);
			return new Vector3D(targetX, targetY, targetZ);
		}
        /// <summary>
        /// 计算站心直角坐标系中的相对坐标在地心坐标系中的坐标
        /// </summary>
        /// <param name="X0">站心原点的地心坐标X</param>
        /// <param name="Y0">站心原点的地心坐标Y</param>
        /// <param name="Z0">站心原点的地心坐标Z</param>
        /// <param name="B0">站心原点的纬度</param>
        /// <param name="L0">站心原点的经度</param>
        /// <param name="H0">站心原点的高度</param>
        /// <param name="rX">相对站心原点的目标坐标X</param>
        /// <param name="rY">相对站心原点的目标坐标Y</param>
        /// <param name="rZ">相对站心原点的目标坐标Z</param>
        /// <returns></returns>
        public static Vector3D CalRelaXYZBasedOnBLH(double X0,double Y0,double Z0,double B0,double L0,double H0,
            double rX,double rY,double rZ)
        {
            double dblD2R = Math.PI / 180.0;
            var AX = -Math.Sin(L0 * dblD2R) * rX + -Math.Sin(B0 * dblD2R) * Math.Cos(L0 * dblD2R) * rY + Math.Cos(B0 * dblD2R) * Math.Cos(L0 * dblD2R) * rZ;
            var AY = -Math.Cos(L0 * dblD2R) * rX + -Math.Sin(B0 * dblD2R) * Math.Sin(L0 * dblD2R) * rY + Math.Cos(B0 * dblD2R) * Math.Sin(L0 * dblD2R) * rZ;
            var AZ = 0 + Math.Cos(B0 * dblD2R) * rY + Math.Sin(B0 * dblD2R) * rZ;
            var relaX = AX + X0;
            var relaY = AY + Y0;
            var relaZ = AZ + Z0;
            return new Vector3D(relaX, relaY, relaZ);
        }
        #endregion
    }

3.转换代码,为了直观变量名使用中文

double 基准点纬度 = 34, 基准点经度 = 112, 基准点高 = 300;
double 目标站心坐标X = 0, 目标站心坐标Y = 100000, 目标站心坐标Z = 8000;
var 站心原点的地心坐标 = CoorTransform.BLHtoXYZ(基准点纬度, 基准点经度, 基准点高);
var 目标的地心坐标 = CoorTransform.CalRelaXYZBasedOnBLH(站心原点的地心坐标.X, 站心原点的地心坐标.Y, 站心原点的地心坐标.Z, 基准点纬度, 基准点经度, 基准点高, 目标站心坐标X, 目标站心坐标Y, 目标站心坐标Z);
var 目标的纬经高 = CoorTransform.XYZtoBLH(目标的地心坐标.X, 目标的地心坐标.Y, 目标的地心坐标.Z);

4.结果

 [注：]

1.直角坐标单位皆为米。

2.本文的站心直角坐标为东、北、天，即X+ => 东，Y+ => 北，Z+ => 天。如果使用北、天、东注意下转换。

引用

[1]凌震莹.大地坐标系与站心地平直角坐标系的坐标转换[J].声学与电子工程,2009(04):31-34.