1.地理空間距離計算面臨的挑戰

打開美團app，不管是篩選團購還是篩選商家，默認的排序項都是“離我最近”或者“智能排序”（如下圖所示）。

不管是“離我最近”還是“智能排序”，都涉及到計算用戶位置與各個團購單子或者商家的距離（注：在智能排序中距離作為一個重要的參數參與排序打分）。以篩選商家為例，北京地區有5~6w個POI（本文將商家稱之為POI），當用戶進入商家頁，請求北京全城+所有品類+離我最近/智能排序時，我們篩選服務需要計算一遍用戶位置與北京全城所有POI的距離。

這種大量計算距離的場景十分消耗資源，從測試來看目前5w個點僅計算一遍距離就需要7ms，而到100w的時候就需要140多ms，隨著數據的快速增長，篩選服務的性能將會非常堪憂。

如何優化距離的計算，進而提高計算速度、降低cpu使用率已經迫在眉睫。美團移動後台團購組在此方向上進行了些許探索，下文將分為4部分展開：1）地理空間距離計算原理；2）Lucene使用的距離計算公式；3）優化方案；4）實際應用。

2 地理空間距離計算原理

地理空間距離計算方法較多，目前我們使用的可以分為兩類：1）球面模型，這種模型將地球看成一個標准球體，球面上兩點之間的最短距離即大圓弧長，這種方法使用較廣，在我們服務端被廣泛使用；2）橢球模型，該模型最貼近真實地球，精度也最高，但計算較為復雜，目前客戶端有在使用，但實際上我們的應用對精度的要求沒有那麼高。

下面著重介紹我們最常用的基於球面模型的地理空間距離計算公式，推導也只需要高中數學知識即可。

經緯度距離計算示意圖

該模型將地球看成圓球，假設地球上有A(ja,wa)，B(jb,wb)兩點（注：ja和jb分別是A和B的經度，wa和wb分別是A和B的緯度），A和B兩點的球面距離就是AB的弧長，AB弧長=R*角AOB（注：角AOB是A跟B的夾角，O是地球的球心，R是地球半徑，約為6367000米）。如何求出角AOB呢？可以先求AOB的最大邊AB的長度，再根據余弦定律可以求夾角。

如何求出AB長度呢？

1）根據經緯度，以及地球半徑R，將A、B兩點的經緯度坐標轉換成球體三維坐標；

2）根據A、B兩點的三維坐標求AB長度;

3）根據余弦定理求出角AOB;

4）AB弧長=R*角AOB.

3.Lucene使用的地理空間距離算法

目前美團團購app後台使用lucene來篩選團購單子和商家，lucene使用了spatial4j工具包來計算地理空間距離，而spatial4j提供了多種基於球面模型的地理空間距離的公式，其中一種就是上面我們推導的公式，稱之為球面余弦公式。還有一種最常用的是Haversine公式，該公式是spatial4j計算距離的默認公式，本質上是球面余弦函數的一個變換，之前球面余弦公式中有cos(jb-ja)項，當系統的浮點運算精度不高時，在求算較近的兩點間的距離時會有較大誤差，Haversine方法進行了某種變換消除了cos(jb-ja)項，因此不存在短距離求算時對系統計算精度的過多顧慮的問題。

1）Haversine公式代碼

public static double distHaversineRAD(double lat1, double lon1, double lat2, double lon2) {
        double hsinX = Math.sin((lon1 - lon2) * 0.5);
        double hsinY = Math.sin((lat1 - lat2) * 0.5);
        double h = hsinY * hsinY +
                (Math.cos(lat1) * Math.cos(lat2) * hsinX * hsinX);
        return 2 * Math.atan2(Math.sqrt(h), Math.sqrt(1 - h)) * 6367000;
    }

2）Haversine公式性能

目前北京地區在線服務有5w個POI，計算一遍距離需要7ms。現在數據增長特別快，未來北京地區POI數目增大到100w時，我們篩選服務僅計算距離這一項就需要消耗144多ms，性能十分堪憂。（注：本文測試環境的處理器為2.9GHz Intel Core i7，內存為8GB 1600 MHz DDR3，操作系統為OS X10.8.3，實驗在單線程環境下運行）

4.優化方案

通過抓棧我們發現消耗cpu較多的線程很多都在執行距離公式中的三角函數例如atan2等。因此我們的優化方向很直接---消減甚至消除三角函數。基於此方向我們嘗試了多種方法，下文將一一介紹。

4.1 坐標轉換方法

1）基本思路

之前POI保存的是經緯度（lon,lat），我們的計算場景是計算用戶位置與所有篩選出來的poi的距離，這裡會涉及到大量三角函數計算。一種優化思路是POI數據不保存經緯度而保存球面模型下的三維坐標（x,y,z），映射方法如下：

x = Math.cos(lat) Math.cos(lon);

y = Math.cos(lat) Math.sin(lon);

z = Math.sin(lat);

那麼當我們求夾角AOB時，只需要做一次點乘操作。比如求(lon1,lat1）和 (lon2,lat2）的夾角，只需要計算x1x2 + y1y2 + z1*z2， 這樣避免了大量三角函數的計算。

在得到夾角之後，還需要執行arccos函數，將其轉換成角度，AB弧長=角AOB*R（R是地球半徑）。

此方法性能如下：

2）進一步優化

美團是本地生活服務，我們的業務場景是在一個城市范圍內進行距離計算，因此夾角AOB往往會比較小，這個時候sinAOB約等於AOB，因此我們可以將 Math.acos(cosAOB)R 改為Math.sqrt(1 - cosAOBcosAOB)*R，從而完全避免使用三角函數，優化後性能如下。

4.2 簡化距離計算公式方法

1）基本思路

我們的業務場景僅僅是在一個城市范圍內進行距離計算，也就是說兩個點之間的距離一般不會超過200多千米。由於范圍小，可以認為經線和緯線是垂直的，如圖所示，要求A（116.8，39,78）和B（116.9，39.68）兩點的距離，我們可以先求出南北方向距離AM，然後求出東西方向距離BM，最後求矩形對角線距離，即sqrt(AMAM + BMBM)。

簡化距離計算示意圖

南北方向AM = R緯度差Math.PI/180.0；

東西方向BM = R經度差Cos<當地緯度數* 180.0="">

這種方式僅僅需要計算一次cos函數。

public static double distanceSimplify(double lat1, double lng1, double lat2, double lng2, double[] a) {
     double dx = lng1 - lng2; // 經度差值
     double dy = lat1 - lat2; // 緯度差值
     double b = (lat1 + lat2) / 2.0; // 平均緯度
     double Lx = toRadians(dx) * 6367000.0* Math.cos(toRadians(b)); // 東西距離
     double Ly = 6367000.0 * toRadians(dy); // 南北距離
     return Math.sqrt(Lx * Lx + Ly * Ly);  // 用平面的矩形對角距離公式計算總距離
    }
}

我們對這個方法的有效性和性能進行驗證。

1.1）有效性驗證

我們首先檢驗這種簡化是否能滿足我們應用的精度，如果精度較差將不能用於實際生產環境。

我們的方法叫distanceSimplify，lucene的方法叫distHaversineRAD。下表是在不同尺度下兩個方法的相差情況。

可以看到兩者在百米、千米尺度上幾乎沒有差別，在萬米尺度上也僅有分米的差別，此外由於我們的業務是在一個城市范圍內進行篩選排序，所以我們選擇了北京左下角和右上角兩點進行比較，兩點相距有260多千米，兩個方法差別17m。從精度上看該優化方法能滿足我們應用需求。

1.2）性能驗證

2）進一步優化

我們看到這裡計算了一次cos這一三角函數，如果我們能消除此三角函數，那麼將進一步提高計算效率。

如何消除cos三角函數呢？

采用多項式來擬合cos三角函數，這樣不就可以將三角函數轉換為加減乘除了嘛！

首先決定多項式的最高次數，次數為1和2顯然都無法很好擬合cos函數，那麼我們選擇3先嘗試吧，注：最高次數不是越多越好，次數越高會產生過擬合問題。

使用org.apache.commons.math3這一數學工具包來進行擬合。中國的緯度范圍在10~60之間，即我們將此區間離散成Length份作為我們的訓練集。

public static double[] trainPolyFit(int degree, int Length){
    PolynomialCurveFitter polynomialCurveFitter = PolynomialCurveFitter.create(degree);
    double minLat = 10.0; //中國最低緯度
    double maxLat = 60.0; //中國最高緯度
    double interv = (maxLat - minLat) / (double)Length;
    List weightedObservedPoints = new ArrayList();
    for(int i = 0; i < Length; i++) {
        WeightedObservedPoint weightedObservedPoint = new WeightedObservedPoint(1,  minLat + (double)i*interv, Math.cos(toRadians(x[i])));
        weightedObservedPoints.add(weightedObservedPoint); 
    }
    return polynomialCurveFitter.fit(weightedObservedPoints);
}
public static double distanceSimplifyMore(double lat1, double lng1, double lat2, double lng2, double[] a) {
     //1) 計算三個參數
     double dx = lng1 - lng2; // 經度差值
     double dy = lat1 - lat2; // 緯度差值
     double b = (lat1 + lat2) / 2.0; // 平均緯度
     //2) 計算東西方向距離和南北方向距離(單位：米)，東西距離采用三階多項式
     double Lx = (a[3] * b*b*b  + a[2]* b*b  +a[1] * b + a[0] ) * toRadians(dx) * 6367000.0; // 東西距離
     double Ly = 6367000.0 * toRadians(dy); // 南北距離
     //3) 用平面的矩形對角距離公式計算總距離
     return Math.sqrt(Lx * Lx + Ly * Ly);
}

我們對此優化方法進行有效性和性能驗證。

2.1）有效性驗證

我們的優化方法叫distanceSimplifyMore，lucene的方法叫distHaversineRAD，下表是在不同尺度下兩個方法的相差情況。

可以看到在百米尺度上兩者幾乎未有差別，在千米尺度上僅有分米的區別，在更高尺度上如72千米僅有5.6m米別，在264千米也僅有8.1米區別，因此該優化方法的精度能滿足我們的應用需求。

2.2）性能驗證

5.實際應用

坐標轉換方法和簡化距離公式方法性能都非常高，相比lucene使用的Haversine算法大大提高了計算效率，然而坐標轉換方法存在一些缺點：

a）坐標轉換後的數據不能被直接用於空間索引。lucene可以直接對經緯度進行geohash空間索引，而通過空間轉換變成三維數據後不能直接使用。我們的應用有附近范圍篩選功能（例如附近5km的團購單子），通過geohash空間索引可以提高范圍篩選的效率；

b）坐標轉換方法增大內存開銷。我們會將坐標寫入倒排索引中，之前坐標是2列（經度和緯度），現在變成3列（x,y,z），在使用中我們往往會將這數據放入到cache中，因此會增大內存開銷；

c）坐標轉換方法增大建索引開銷。此方法本質上是將計算從查詢階段放至到索引階段，因此提高了建索引的開銷。

基於上述原因我們在實際應用中采用簡化距離公式方法（通過三次多項式來擬合cos三角函數），此方法在團購篩選和商家篩選的距離排序、智能排序中已經開始使用，與之前相比，篩選團購時北京全城美食品類距離排序響應時間從40ms下降為20ms。

6 參考資料

http://blog.csdn.net/liminlu0314/article/details/8553926

http://www.movable-type.co.uk/scripts/gis-faq-5.1.html