基于WGS-84坐标系的GPS坐标下两点间距离计算 VB.net实现

精度更高的一种计算方式.

WGS－84坐标系是一种国际上采用的地心坐标系。坐标原点为地球质心，其地心空间直角坐标系的Z轴指向国际时间局（BIH）1984.0定义的协议地极（CTP）方向，X轴指向BIH1984.0的协议子午面和CTP赤道的交点，Y轴与Z轴、X轴垂直构成右手坐标系，称为1984年世界大地坐标系。这是一个国际协议地球参考系统（ITRS），是目前国际上统一采用的大地坐标系。

] ''' <summary>
''' 通过经纬度计算两点间距离
''' </summary>
''' <param name="lat1">出发点纬度</param>
''' <param name="lon1">出发点经度</param>
''' <param name="lat2">目标点纬度</param>
''' <param name="lon2">目标点经度</param>
''' <returns>距离 单位()</returns>
''' <remarks>采用WGS-84 椭球模型</remarks>
Public Function GPSDistance(ByVal lat1 As Double, ByVal lon1 As Double, ByVal lat2 As Double, ByVal lon2 As Double) As Double
Dim a As Double = 6378137.0, b As Double = 6356752.3142, f As Double = 1 / 298.257223563
Dim rlat1 As Double = GetRad(lat1)
Dim rlon1 As Double = GetRad(lon1)
Dim rlat2 As Double = GetRad(lat2)
Dim rlon2 As Double = GetRad(lon2)
Dim L As Double = GetRad(lon2 - lon1) 'lonDiff
Dim U1 As Double = Math.Atan((1 - f) * Math.Tan(rlat1))
Dim U2 As Double = Math.Atan((1 - f) * Math.Tan(rlat2))
Dim sinU1 As Double = Math.Sin(U1)
Dim cosU1 As Double = Math.Cos(U1)
Dim sinU2 As Double = Math.Sin(U2)
Dim cosU2 As Double = Math.Cos(U2)
Dim lambda As Double = L, lambdap As Double = 0.0
Dim iterLimit As Integer = 100
Dim sinLambda As Double
Dim cosLambda As Double
Dim sinSigma As Double
Dim cosSigma As Double
Dim Sigma As Double
Dim sinAlpha As Double
Dim cosSqAlpha As Double
Dim cos2SigmaM As Double
Dim C As Double
Do
sinLambda = Math.Sin(lambda)
cosLambda = Math.Cos(lambda)
sinSigma = Math.Sqrt((cosU2 * sinLambda) ^ 2 + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) ^ 2)
If sinSigma = 0.0 Then
Return 0.0
End If
cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
Sigma = Math.Atan2(sinSigma, cosSigma)
sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
cosSqAlpha = 1.0 - sinAlpha ^ 2
cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha
If (Double.IsNaN(cos2SigmaM)) Then cos2SigmaM = 0
C = f / 16.0 * cosSqAlpha * (4.0 + f * (4.0 - 3.0 * cosSqAlpha))
lambdap = lambda
lambda = L + (1 - C) * f * sinAlpha * (Sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2.0 * cos2SigmaM ^ 2)))
iterLimit = iterLimit - 1
Loop While ((Math.Abs(lambda - lambdap) > 0.000000000001) And (iterLimit > 0))
If iterLimit = 0 Then Return Double.NaN
Dim uSq As Double = cosSqAlpha * (a ^ 2 - b ^ 2) / (b ^ 2)
Dim aA As Double = 1 + uSq / 16384.0 * (4096.0 + uSq * (-768.0 + uSq * (320.0 - 175.0 * uSq)))
Dim bB As Double = uSq / 1024.0 * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq)))
Dim deltaSigma As Double = bB * sinSigma * (cos2SigmaM + bB / 4.0 * (cosSigma * (-1 + 2 * cos2SigmaM ^ 2) - _
bB / 6.0 * cos2SigmaM * (-3 + 4 * sinSigma ^ 2) * (-3 + 4 * cos2SigmaM ^ 2)))
Dim s As Double = b * aA * (Sigma - deltaSigma)
Return s
End Function

''' <summary>
''' 通过经纬度计算两点间距离
''' </summary>
''' <param name="p0">起始点坐标</param>
''' <param name="p1">目的点坐标</param>
''' <returns>两点间距离,单位(米)</returns>
''' <remarks>采用WGS-84 椭球模型</remarks>
Public Function GPSDistance(ByVal p0 As StationPoint, ByVal p1 As StationPoint) As Double
Return GPSDistance(p0.X, p0.Y, p1.X, p1.Y)
End Function
''' <summary>
''' 获取角度的弧度制表示
''' </summary>
''' <param name="DEG">角度制表示(单位度)</param>
''' <returns></returns>
''' <remarks></remarks>
Function GetRad(ByVal DEG As Double) As Double
Return (Math.PI * DEG / 180)
End Function