使用 GPS（WGS-84）数学进行位置邻近搜索的方法

• 下载示例类 - 4.6 Kb

引言

这提供了一种快速的邮政编码邻近搜索，并避免了数值方法所需的大量浮点运算 (FLOPS)。我不会称之为真正的距离搜索，因为它是对我们非解析现实的一种解析近似，但它比你在网上找到的大多数假设地球是完美球体的代码要好得多。

对于那些希望大致了解商业实体之间距离，但又不想每年花费 8000 美元购买 MapPoint API，或者根本不需要那么高细节或精度的客户来说，这非常方便。

美国的邮政编码数据每年只需 40 美元即可购买。不要相信联合国在其 LOCODE 系统中使用的经纬度，它们非常不准确。如果你需要国际位置数据，你将需要寻找其他来源。

还要记住，邮政编码的中心点经纬度可能会落在实际邮政编码区域之外，这取决于其几何形状。在实现此功能时，请记住其主要目的是快速进行邻近搜索，并且它基于 GPS 解析器背后的解析模型。如果你的应用程序需要真正的 GPS 级别精度，我建议你查看这篇 Code Project 文章，它使用了现有的 GPS 解析器用于商业目的。

理论

地球是一个扁球体，因为自转导致了扁平化，正如我用 MATLAB 所描绘的那样。图像被挖空，以便你可以清楚地看到它与网上随处可见代码中使用的球体模型之间的区别。

你在网上能找到的许多搜索方法都有一个大问题，那就是这种（球体）假设会引入巨大的误差。你可以想象在上面的图像上叠加一个球体，然后取一个垂直于目标点的平面。首先，用球体模型来表示地球“表面”本身就已经是极大的夸大了。然后，当你取一个法平面时（像许多网上的代码示例那样，用于计算最大和最小经纬度），法平面上的那些点会被进一步夸大。

通过使用 WGS-84 中认可并用于 GPS 的数学方法，我们可以得到更好的近似值。这里没有考虑海拔高度，因为我假设你没有地形数据。这是椭球体解析模型中的一个误差来源，但它也存在于常见的球体模型中，从而进一步加剧了其数学上的缺陷。

下面描述了此 WGS-84 椭球体模型的数学误差。

请记住，在所有情况下，距离都是鸟瞰图距离，因此你需要在搜索结果中明确说明这一点，以免用户将此邻近数值与道路/旅行距离混淆。

现在你每年只需 40 美元就可以获得邮政编码的经纬度数据，无需支付更复杂的服务费用来维护道路距离、地址和海拔数据。这是一个精确的邻近搜索，其性能远超典型的“完美球体”模型，而且不会让你的处理器不堪重负。请下载单个源代码类文件以了解所有详细信息。本文中我将只阐述概念。

Extension

大搜索半径- 如果你需要搜索超过地球 1/4 的范围，我建议将此算法分块处理，而不是抛出异常。你可以分段处理搜索区域，使用这里的相同分析方法，并在完成后合并结果。你需要确定一个最佳频率来创建新的曲率半径（然后使用由这些半径定义的局部椭圆）。通过这种方法，你可以将曲线平面分解为定义明确的小块，而无需进行全面的微分数值近似，后者会消耗更多的浮点运算。我猜这种扩展纯粹是学术性的。如果你的系统要求搜索超过地球 1/4 的范围，那么它可能本身就是一个要求极高精度的系统，并且它的客户可能每年有 8000 美元可以花在 MapPoint API 或其他企业解决方案上。不过……我们都知道，谁能预测客户接下来会要求什么呢？:-) 所以我认为为这个问题给你留下一些线索是明智的。如果你走这条路，请确保准确地对误差进行建模，并向你的客户说明，以防你需要从这个“姜饼屋”中找到回家的路。

访问海拔数据- 如果你有海拔数据，可以修改这些封闭形式的方程以获得更高的精度。你可以参考 WGS-84 手册 [927 Kb] 来推导你的方程。请查看附录中关于坐标系和大地测量学的部分。

Data

我假设你会为自己的数据架构编写脚本，所以我替换了我的存储过程调用和数据层调用。你可能还想在网上找一个为了排序目的而将 DataView 转换为 DataTable 的类。网上有一个非常流行的 C# 函数，很容易找到。

代码

这些方法非常直接。它们有很好的注释，所以我不会在这里深入探讨太多细节。你将需要自己获取数据。联合国维护的“LOCODE”数据对国际业务很有用。你会看到在 location 对象中有一些属性，我是直接从该数据中获取的。我建议在你的架构中使用以下表格：

• 区域
• Countries（国家）
• Subdivisions（行政区划，州、省、郡等的国际术语）
• 货币
• Location Entry Criteria（位置条目标准）
• Counties（县）
• Location（位置，城市/城镇名称和特征，可选的县外键）
• Postal Codes（邮政编码，带有中心点经纬度，LocationID 外键）

我所有感兴趣的代码都在 Location 对象中，它能提供我想知道的关于地球上某个地点或区域的一切信息。

#region Declarations
    private int regionID;
    private string regionName;
    private int countryID;
    private string countryCode;
    private string countryName;
    private int primaryCurrencyID;
    private string primaryCurrencyCode;
    private string primaryCurrencyName;
    private int secondaryCurrencyID;
    private string secondaryCurrencyCode;
    private string secondaryCurrencyName;
    private int subdivisionID;
    private string subdivisionCode;
    private string subdivisionName;
    private string subdivisionType;
    private int countyID;
    private string countyName;
    private int locationID;
    private string cityName;
    private bool hasPort;
    private bool hasRailTerminal;
    private bool hasRoadTerminal;
    private bool hasAirport;
    private bool hasBorderCrossing;
    private bool hasPostalExchangeOffice;
    private string postalCode;
    private double longitude;
    private double latitude;
    private string address;
    private string addressee;
    //private string language;
    //private string productName;
    #endregion

你也可以创建这些对象的集合，因为我使用了 CollectionBase。

#region Collections
    /// <SUMMARY>
    /// Collection of Location objects
    /// </SUMMARY>
    public Location this[int i]
    {
        get    {return (Location)List[i];}
        set    {List[i] = value;}
    }

    /// <SUMMARY>
    /// Adds a Location object to the collection
    /// </SUMMARY>
    /// <PARAM name="L">Completed Location object</PARAM>
    public void AddLocation(Location L)
    {
        List.Add(L);
    }

    /// <SUMMARY>
    /// Removes a Location object from the collection
    /// </SUMMARY>
    /// <PARAM name="L">Completed Location object</PARAM>
    public void RemoveLocation(Location L)
    {
        List.Remove(L);
    }
    #endregion

公式中使用的常量

#region Constants
    //Equatorial radius of the earth from WGS 84 
    //in meters, semi major axis = a
    internal int a = 6378137;
    //flattening = 1/298.257223563 = 0.0033528106647474805
    //first eccentricity squared = e = (2-flattening)*flattening
    internal double e = 0.0066943799901413165;
    //Miles to Meters conversion factor (take inverse for opposite)
    internal double milesToMeters = 1609.347;
    //Degrees to Radians converstion factor (take inverse for opposite)
    internal double degreesToRadians = Math.PI/180;
    #endregion

现在你可以使用搜索邮政编码来返回作为搜索椭圆中心点的经纬度。曲率半径函数决定了你搜索纬度的参考椭球体。如果搜索参数严重影响了精度，则会抛出异常。

函数：public DataTable FindNearbyLocations(string centroidPostalCode, double searchRadius, bool isMetric)。

//lat naught and lon naught are the geodetic parameters in radians
double lat0 = 
   Convert.ToDouble(centroidDT.Rows[0]["Latitude"])*degreesToRadians;
double lon0 = 
   Convert.ToDouble(centroidDT.Rows[0]["Longitude"])*degreesToRadians;

//Find reference ellipsoid radii
double Rm = calcMeridionalRadiusOfCurvature(lat0);
double Rpv = calcRoCinPrimeVertical(lat0);

//Throw exception if search radius is greater than 1/4 of globe and 
//thus breaks accuracy of model (mostly pertinent for russia, alaska, 
//canada, peru, etc.)
if (Rpv*Math.Cos(lat0)*Math.PI/2 < searchRadius)
{
    throw new ApplicationException("Search radius was too great to " + 
                       "achieve an accurate result with this model.");
}

//Find boundaries for curvilinear plane that encloses search ellipse
double latMax = (searchRadius/Rm+lat0)/degreesToRadians;
double lonMax = (searchRadius/(Rpv*Math.Cos(lat0))+lon0)/degreesToRadians;
double latMin = (lat0-searchRadius/Rm)/degreesToRadians;
double lonMin = (lon0-searchRadius/(Rpv*Math.Cos(lat0)))/degreesToRadians;

//Return postal codes in curvilinear plane
SqlParameter one2 = new SqlParameter("@latMin",System.Data.SqlDbType.Decimal, 
                                          9, ParameterDirection.Input, false, 
                                 ((System.Byte)(18)), ((System.Byte)(5)), "", 
                                 System.Data.DataRowVersion.Current, latMin);
SqlParameter two2 = new SqlParameter("@latMax",System.Data.SqlDbType.Decimal, 
                                          9, ParameterDirection.Input, false, 
                                 ((System.Byte)(18)), ((System.Byte)(5)), "", 
                                 System.Data.DataRowVersion.Current, latMax);
SqlParameter three2 = new SqlParameter("@lonMin",System.Data.SqlDbType.Decimal, 
                                            9, ParameterDirection.Input, false, 
                                   ((System.Byte)(18)), ((System.Byte)(5)), "", 
                                   System.Data.DataRowVersion.Current, lonMin);
SqlParameter four2 = new SqlParameter("@lonMax",System.Data.SqlDbType.Decimal, 9, 
                                                 ParameterDirection.Input, false, 
                                     ((System.Byte)(18)), ((System.Byte)(5)), "", 
                                     System.Data.DataRowVersion.Current, lonMax);
object[] Param2 = new object[4]{one2,two2,three2,four2};
DataTable resultPostalCodesDT = new DataTable("ResultPostalCodes");
resultPostalCodesDT = helper.ReadOnlySproc(helper.YOUR_CONNECTION,
                                                      "dbo.[YOUR_SPROC]",Param2);
resultPostalCodesDT.Columns.Add("DistanceToCentroid",
                                           System.Type.GetType("System.Double"));

//Now calc distances from centroid, and remove results that were returned 
//in the curvilinear plane, but are outside of the ellipsoid geodetic
if (resultPostalCodesDT.Rows.Count > 0)
{
  for (int i=0;i<resultPostalCodesDT.Rows.Count;i++)
  {
    resultPostalCodesDT.Rows[i]["DistanceToCentroid"] = calcDistanceLatLons(Rm,
                                                                  Rpv,lat0,lon0,
      Convert.ToDouble(resultPostalCodesDT.Rows[i]["Latitude"])*degreesToRadians,
      Convert.ToDouble(resultPostalCodesDT.Rows[i]["Longitude"])*degreesToRadians);
    if (Convert.ToDouble(resultPostalCodesDT.Rows[i]["DistanceToCentroid"]) > 
                                                                     searchRadius)
    {
          resultPostalCodesDT.Rows.RemoveAt(i);
          i--;
    }
  }
  if (isMetric == false)
  {
      foreach(DataRow r in resultPostalCodesDT.Rows)
      {
          r["DistanceToCentroid"] = 
                       Convert.ToDouble(r["DistanceToCentroid"])/milesToMeters;
      }
  }
  resultPostalCodesDT.DefaultView.Sort = "DistanceToCentroid";
  DataTable finalResults = new DataTable("finalResults");
  finalResults = YOUR_CUSTOM_UTILITY_CLASS.ConvertDataViewToDataTable(
                                               resultPostalCodesDT.DefaultView);
  return finalResults;
}
else
{
  return null;
}

以下是涉及曲率半径和距离计算的函数，供你参考。

/// <summary>
    /// Calculates the Radius of curvature in the 
    /// prime vertical for the reference ellipsoid
    /// </summary>
    /// <remarks>
    /// This is the vector that defines the normal surface 
    /// to any point on the ellipsoid.  It extends from
    /// from the polar axis to that point.  It is used for the 
    /// longitude, in differentiation of east distances, dE
    /// </remarks>
    /// <param name="lat0">Geodetic latitude in radians</param>
    /// <returns>Length of radius of curvature in the 
    /// prime vertical</returns>
    public double calcRoCinPrimeVertical(double lat0)
    {
        double Rn = a/Math.Sqrt(1-e*Math.Pow(Math.Sin(lat0),2));
        return Rn;
    }

    /// <summary>
    /// Calculates the meridional radius of 
    /// curvature for the reference ellipsoid
    /// </summary>
    /// <remarks>
    /// This is the radius of a circle that fits the earth 
    /// curvature in the meridian at the latitude chosen.
    /// It is used for latitude, in differentiation of 
    /// north distances, dN
    /// </remarks>
    /// <param name="lat0">Geodetic latitude in radians</param>
    /// <returns>Length of meridional radius of 
    /// curvature</returns>
    public double calcMeridionalRadiusOfCurvature(double lat0)
    {
        double Rm = 
               a*(1-e)/Math.Pow(1-e*(Math.Pow(Math.Sin(lat0),2)),1.5);
        return Rm;
    }
    
    /// <summary>
    /// Calculates the true birds-eye view length of the arc 
    /// between two positions on the globe using parameters from WGS 84,
    /// used in aviation and GPS.
    /// </summary>
    /// <remarks>
    /// An accurate BIRDS EYE numerical approximation with error 
    /// approaching less than 10 feet at a 50 miles search radius
    /// and only 60 ft at 400 mile search radius.  Error is on the 
    /// order of (searchRadius/equatorialRadius)^2.  
    /// Only accurate for distances
    /// less than 1/4 of the globe (~10,000 km at the equator, 
    /// and approaching 0 km at the poles).
    /// Geoid height above sea level is assumed to be zero, which is 
    /// the only deviation from GPS, and another source of error.
    /// </remarks>
    /// <param name="Rm">Meridional Radius of Curvature 
    /// at the centroid latitude</param>
    /// <param name="Rpv">Radius of Curvature in the Prime 
    /// Vertical at the centroid latitude</param>
    /// <param name="lat0">Centroid latitude</param>
    /// <param name="lon0">Centroid longitude</param>
    /// <param name="lat">Destination latitude</param>
    /// <param name="lon">Destination longitude</param>
    /// <returns>Distance in meters from the arc between 
    /// the two points on the globe</returns>
    public double calcDistanceLatLons(double Rm, double Rpv, double lat0, 
                                         double lon0, double lat, double lon)
    {
        double distance = Math.Sqrt(Math.Pow(Rm,2)*Math.Pow(lat-lat0,2)+
              Math.Pow(Rpv,2)*Math.Pow(Math.Cos(lat0),2)*Math.Pow(lon-lon0,2));
        return distance;
    }