如何计算浮点支持较差的处理器上GPS坐标之间的距离?

我需要计算GPS坐标之间的距离来计算行进距离。 我已经尝试了Haversine和Vincenty算法,它们在我的桌面PC上运行良好,但是当我将代码移植到dsPIC时,由于缺少浮点精度,它们返回0表示接近(几米之内)的点和罪和cos的不良实现。

对于我的用例,我的分数不会超过10米,并且都将相距不超过10公里。 我尝试了以下算法,结果似乎没问题:

double dist(double latA, double lonA, double latB, double lonB) { double latD = fabs(latA - latB) * 111.3; double lonD = fabs(lonA - lonB) * 111.3 * cos(latA * 3.14159265 / 180); return sqrt(latD*latD + lonD*lonD) * 1000; } 

假设每1°的距离是111.3km,我用毕达哥拉斯定理来计算距离。 有没有简单的方法来改进我的算法? 或者是否有其他算法不依赖于高度准确的sin / cos?

用于Marine AIS系统(在IEC61193-4中规定)的公认算法是Rhumb线算法。 我已经使用Anthony Williams的定点数学库 (使用CORDIC算法)成功地在目标上实现了这一点 ,我相信通常会比软件浮点数提供5倍的性能提升。

然而,库是C ++而不是C,这使得它易于使用,因为广泛的运算符重载,但可能不是你想要的。 值得考虑为您的C代码使用C ++编译,只是为了这个库的好处。 当然,问题在于Microchip的C31编译器奇怪地不支持C ++。

但需要注意的是,log()函数的查找表太短了一个值,并且最后需要一个额外的元素,其值为零。 安东尼在我找到它后证实了这一点,但我不相信他已经更新了下载。

无论哪种方式,答案可能是使用固定点和CORDIC。

要解析1m的经度或赤道弧,你需要8位精度,所以单精度浮点数不足,使用双精度会大大减慢速度。 检查MikroElectronica的C用户手册显示编译器只支持单精度 – floatdouble精度和long double精度都是32位,因此在任何情况下无法使用内置FP类型实现所需的精度编译器。

如果有任何用处,这是使用Anthony的库的我的Rhumb Line代码:

标题:

 #if !defined cRhumbLine_INCLUDE #define cRhumbLine_INCLUDE #include "fixed.hpp" //! Rhumb Line method for distance and bearing between two geodesic points class cRhumbLine { public: //! @brief Default constructor //! //! Defines a zero length line, bearing zero cRhumbLine() : m_update_bearing(false), m_update_distance(false), m_distance(0), m_bearing(0) {} //! @brief Constructor defining a line //! //! @param startLatDeg Start latitude in degrees, negative values are south of equator //! @param startLonDeg Start longitude in degrees, negative values are west of meridian. //! @param endLatDeg End latitude in degrees, negative values are south of equator //! @param endLonDeg End longitude in degrees, negative values are west of meridian. cRhumbLine( fixed startLatDeg, fixed startLonDeg, fixed endLatDeg, fixed endLonDeg ) { define( startLatDeg, startLonDeg, endLatDeg, endLonDeg ) ; } //! @brief Define a start/ent point //! //! @param startLatDeg Start latitude in degrees, negarive values are south of equator //! @param startLonDeg Start longitude in degrees, negative values are west of meridian. //! @param endLatDeg End latitude in degrees, negarive values are south of equator //! @param endLonDeg End longitude in degrees, negative values are west of meridian. void define( fixed startLatDeg, fixed startLonDeg, fixed endLatDeg, fixed endLonDeg ) ; //! @brief Get start-end distance in nautical miles //! @return Start-end distance in nautical miles. fixed distanceNm() { return distanceMetres() / ONE_NM_IN_METRE ; } //! @brief Get start-end distance in metres. //! @return Start-end distance in metres. fixed distanceMetres() ; //! @brief Get start-end bearing in degreed. //! @return Start-end bearing in degreed (0 <= x < 360). fixed bearingDeg() ; private: static const int ONE_NM_IN_METRE = 1852 ; bool m_update_bearing ; bool m_update_distance ; fixed m_distance ; fixed m_bearing ; fixed m_delta_phi ; fixed m_delta_lon ; fixed m_delta_lat ; fixed m_lat1_radians ; } ; #endif 

身体:

 #include "cRhumbLine.h" void cRhumbLine::define( fixed startLatDeg, fixed startLonDeg, fixed endLatDeg, fixed endLonDeg ) { fixed lat1 = startLatDeg / 180 * fixed_pi ; fixed lon1 = startLonDeg / 180 * fixed_pi ; fixed lat2 = endLatDeg / 180 * fixed_pi ; fixed lon2 = endLonDeg / 180 * fixed_pi ; m_update_bearing = true ; m_update_distance = true ; m_delta_phi = log( tan( lat2 / 2 + fixed_quarter_pi ) / tan( lat1 / 2 + fixed_quarter_pi ) ) ; m_delta_lat = lat1 - lat2 ; m_delta_lon = lon1 - lon2 ; m_lat1_radians = lat1 ; // Deal with delta_lon > 180deg, take shorter route across meridian if( m_delta_lon.abs() > fixed_pi ) { m_delta_lon = m_delta_lon > 0 ? -(fixed_two_pi - m_delta_lon) : (fixed_two_pi + m_delta_lon) ; } } fixed cRhumbLine::distanceMetres() { if( m_update_distance ) { static const fixed mean_radius = 6371009 ; // in metres. Source: International Union of Geodesy and Geophysics (IUGG) fixed q = m_delta_phi != 0 ? m_delta_lat / m_delta_phi : cos( m_lat1_radians ) ; // Deal with lines with constant latitude m_distance = sqrt( ( sqr(m_delta_lat) + sqr(q) * sqr(m_delta_lon) ) ) * mean_radius ; m_update_distance = false ; } return m_distance ; } fixed cRhumbLine::bearingDeg() { if( m_update_bearing ) { static const fixed mean_radius = 6371009 ; // in metres. Source: International Union of Geodesy and Geophysics (IUGG) m_bearing = atan2( m_delta_lon, m_delta_phi ) / fixed_pi * 180 ; if( m_bearing == -180 ) { m_bearing == 0 ; } else if( m_bearing < 0 ) { m_bearing += 360 ; } m_update_bearing = false ; } return m_bearing ; } 

一些评论:

  • 您需要指定计算的范围和精度要求。 范围和精度对于确定用于计算余弦的方法非常重要 。 此外,如果纬度和经度的相对差异与杆的角距离相比较小,那么你发布的毕达哥拉斯近似法效果很好。 如果纬度不靠近,你的伪毕达哥拉斯算法在高纬度地区将无法正常工作。 (例如纬度为43.001和43.002它会很好,但不是89.961和89.962)

  • 需要根据它们的圆度来计算经度 – 你的算法将在国际日期线附近失败,但是通过采用纵向差对称 – 模360可以很容易地解决这个问题,其中smod(x,m) = mod(x+m/2,m)-m/2 。 (例如-179.5 – +179.5 = -359度,但是如果你计算smod(-359,360)你得到+ 1度。)

  • 在设计时,充分利用您的PC。 您有一个非常强大的计算器,您可以评估大量测试点的高精度答案和近似值,并了解它们的不同之处,以评估准确性。 如果您在此信息中推导出模式,则可以使用它来进行二阶近似以提高精度。


更新 :您声明您的范围/精度要求为+/- 60度(在一个半球的范围内没有优势)和1%的准确度。 在该范围内,以度为单位的cos(x)的良好近似是c 2 (x)= 0.9942-1.39 * 10 -4 * x 2 = 0.9942 – (0.01179x) 2 ; 它在此范围内的误差最大值为0.006。

如果你想要更好的精度,使用4次多项式(c 4 (x)= 0.999945-(0.01233015x) 2 +(0.007778x) 4在这个小于6×10 -5的范围内有一个最大误差,但对它更敏感参数错误和算术错误)或拆分成多个范围。

您可能想尝试使用预先计算的表来进行sin和cos。 它使用更多内存,可以将通用处理器(不是你的情况)上的缓存丢弃,但在处理器上具有尽可能高的准确性并且速度非常快。

你是一个固定点DSP(有效),所以你可能想看看定点function; 他们可能表现更好。

事实certificate,Microchip有一个定点库: http : //www.microchip.com/stellent/idcplg ?IdcService = SS_GET_PAGE&nodeId = 2680&dDocName = en552208我不知道它会有多大帮助 – 它可能缺乏你的精确度需要。

这是一个如何自己做的例子: http : //www.coranac.com/2009/07/sines

返回正轨 – Microchip页面建议他们的编译器和库兼容IEEE-754单精度和双精度。 除非他们说的是半真半假,并且正在使用半精度(16位)格式。 如果你还没有得到你需要的结果,我会考虑提交错误报告。