空間索引 - GeoHash演算法及其實現優化

枕邊書發表於2017-05-16
索引演算法優化

前言

上篇部落格中提到了空間索引的用途和多種資料庫對空間索引的支援情況，那麼在應用層以下，好學的小夥伴應該會考慮空間索引的實現原理了。

目前空間索引的實現有 R樹和其變種GIST樹、四叉樹、網格索引等。 網格索引不再多提，使用普通的hash表儲存地點和風格之間的對映來實現。今天要介紹的GeoHash演算法實現的空間索引，雖然是以B樹實現，但我認為它也借用網格索引的一部分思想。

GeoHash

原理

GeoHash 演算法的原理說起來是很簡單的，如下圖：

  1. 從橫向上將整個方形紙分為左右兩份，左側部分為標記為 0， 右側部分標記為 1；
  2. 再將紅點所在的部分劃分為左右兩塊，再對紅點位置做同樣的標識，最後得出紅點在橫向上的標識為 10;
  3. 在縱向上對方形紙做同樣的劃分，左側標識為0，右側標識為 1，得出紅點位置在縱向上的標識為 01;
  4. 將橫向標識和縱向標識合併，規則為 縱向在奇數位，橫向在偶數位 (也可縱橫相反，但要在整個系統內保持一致)，得出紅點在方形紙上的標識為 1001;

只標記一個方格顯得看不出什麼規律，如果我們把這些都空格都標識後會發現 被劃分在角落裡的四個方格會有同樣的字首，如下圖所示。

同樣的字首意味著可以使用 B樹 索引查詢有相同字首的點作為附近的點，GeoHash 演算法便是這些同樣的字首上面做文章。

墨卡託投影

墨卡託投影，是正軸等角圓柱投影。由荷蘭地圖學家墨卡託(G.Mercator)於1569年創立。假想一個與地軸方向一致的圓柱切或割於地球，按等角條件，將經緯網投影到圓柱面上，將圓柱面展為平面後，即得本投影。墨卡託投影在切圓柱投影與割圓柱投影中，最早也是最常用的是切圓柱投影。

墨卡託投影簡單地說，就是可以 把整個地球平面作為一個正方形來處理，當然地球平面不是嚴格的正方形，此投影在兩極附近的點會有誤差，本文專注於原理，糾偏就不多提了（我也不懂，逃）。

實現

按照墨卡託投影的平面，我們可以按照上面劃分方格紙的方式來將整個地球表面劃分為各個小方格。

如(116.276349, 40.040875)這個點的經度劃分：

  1. 經度在 [-180,0) 範圍內的標識為0，經度範圍在 [0, 180) 度的標識為 1;
  2. 繼續劃分，經度範圍在 [0,90) 的標識為 0，經度範圍在 [90,180) 的標識為 1;
  3. 這樣，我們劃分 20 次，方格的精度（見文末對照表）已達到 2m，得到經度的標識二進位制串為11010010101011110111;
  4. 對緯度同樣劃分，得到緯度的標識二進位制串為10111000111100100111;
  5. 我們對它組合，得到40位的二進位制串11011 01110 00010 01110 11100 10111 01001 11111;
  6. 我們將這個二進位制串使用 base32編碼（原理同base64，可以見我的另一篇文章：WEB開發中的字符集和編碼，位編碼對映表見下），得到 GeoHash 編碼為 3OCO4XJ7;

那麼GeoHash編碼字首為 3OCO4XJ7的地理點就是離 (116.276349, 40.040875)兩米內的點。如果我們把地理位置點和其GeoHash編碼存入資料庫的話，我們要查詢 附近兩米點的點，只需要限定條件 geo_code like '3OCO4XJ7%'就行了;

邊界點問題

可是最簡版的 GeoHash 還有一個弱點，如下圖：

如果每個方格的精度為 2km，那麼我們直接按照字首查詢紅點附近 2km 的點是查詢不到離它很近的黑點的。

要解決這個問題，我們就需要所其周邊八個方格也考慮上，將自身方格和周邊八個方格內的點都遍歷一次，再返回符合要求的點。那麼如何知道周邊方格的字首呢？

仔細觀察相鄰方格，我們會發現兩個小方格會在 經度或緯度的二進位制碼上相差1；我們通過 GeoHash 碼反向解析出二進位制碼後，將其經度或緯度（或兩者）的二進位制碼加一，再次組合為 GeoHash 碼。

Redis的GEO函式

問題

我們常見的需求是查詢 n米 範圍內的點，那麼 n米 與 GeoHash 碼位數之間的對映如何實現呢？由於 GeoHash 碼是由5位二進位制碼組成，每少一位，精度就會損失 2e(5/2)。

方法當然有的，我們將二進位制GeoHash碼直接索引就可以，但很長的索引長度會導致 B樹 索引查詢效率會迅速下降。

方案

於是我們接著尋找解決方案，既然使用 base32 轉換為 32進位制碼 會不好控制精度，保持二進位制又導致索引長度過長，那麼進位制位數和索引長度有沒有一個平衡呢？

另外 Redis 的 sorted set 支援 64位 的 double 型別的 score，我們把二進位制的 GeoHash 碼轉為十進位制放入 Redis 的 sorted set 中，不是可以實現 log(n)的查詢效率了麼。

說實話第一次看到 Redis 的 GEO 系列函式的時候我的內心是崩潰的，原來自己感覺極其良好的設計早已被人實現了（雖然這種情況經常出現）。。。

當然不能就這麼算了，於是我使用PHP造了一遍輪子。。。

主要步驟如下：

程式碼實現

實現中我將 GeoHash 的最大精度設定為26位，此時它的距離精度為 0.3m。當然我們也可以充分利用 Redis 的 sorted set 的 score，設定精度為 32 位，剛好使用它的 double 型別。

放上GitHub原始碼地址：空間索引-GeoHash

資料入庫：

將經緯度通過 GeoHash 演算法獲取到二進位制 GeoHash 碼，並將其轉成十進位制作為這個點的 score 存入 Redis 的 sorted set;

// GeoHash核心方法 傳入float型別的度數和其對應的範圍，經度和緯度公用方法
public function getBits($loc, $range, $level = self::LEVEL_MAX) {
    $bits = '';
    for ($i = 0; $i < $level; $i++) {
        $mid = ($range['min'] + $range['max']) / 2;
        if ($loc < $mid) {
            $bits .= '0';
            $range = ['min' => $range['min'], 'max' => $mid];
        } else {
            $bits .= '1';
            $range = ['min' => $mid, 'max' => $range['max']];
        }
    }

    return $bits;
}

另外 php 的 bindec($bin_str) 方法能快速把二進位制字串轉為十進位制數字。

根據查詢範圍半徑獲取精度

上文說過，精度是由地圖的劃分次數決定的，劃分次數多了，範圍就小了，查詢的出的資料就不全；劃分次數少了，範圍就會大了，我們對資料過濾時就會有過多的損耗。

private function getLevel($range_meter){
    $level = 0;
    $global = self::MERCATOR_LENGTH;
    while ($global > $range_meter) {
        $global /= 2;
        $level++;
    }

    return $level;
}

上面程式碼的思想來自redis geo函式原始碼，真的很巧妙。

在墨卡託投影下，地球的表面可以作為一個正方形來看，它的邊是地球周長中最長的一個。而學過初中地理的我們知道：“地球是一個兩極稍扁，赤道略鼓的球體”，那麼它最長的一個周長就是赤道周長了，於是我們得知墨卡託投影的長邊為 2*PI*R=40075452.74M;

於是我們拿正方形的一個邊來不停地進行二次劃分，直到劃分後的結果剛好比範圍半徑長，那麼它構成的一個方塊，便是我們需要的方格。

資料查詢

資料查詢時，我們需要獲取中間方塊的最小 score 值和其範圍，最小 score 值很簡單，直接將二進位制位不足52位的在後面補0。

此外，為了避免邊界點問題，我們還需要把周圍八個方格的 score 值範圍也獲取到。

我們在劃分地圖時，每多劃分一次，會新增經度和緯度兩個二進位制位，在精度最高時，那麼每一個方格的最大值和最小值之間差1。由此，我們通過下面的方法獲取到一個方格的最大和最小 score 值之差。

private function getLevelRange($level) {
    $range = pow(2, 2 * (self::LEVEL_MAX - $level));

    return $range;
}

再由上面提過的邊界點問題的解決方案，獲取到周邊八個方格的最小 score 值。

使用 Redis 的ZRANGEBYSCORE key hashInt hashInt+range命令將這九個方格內的點全部取到，再遍歷九個方格，將距離不符合的資料過濾掉。

小結

花費了十多個小時，總算將 GeoHash 完全整體了一遍，完全理解 GeoHash 並沒有想像中的那麼簡單。除了 GeoHash，四叉樹和R樹據說查詢效率會更高，有時間再研究一下。

如果您覺得本文對您有幫助，可以點選下面的 推薦 支援一下我。部落格一直在更新，歡迎 關注 。

參考：

GeoHash核心原理解析

Redis GEO 原始碼註釋

GeoHash位數精度對照表（wiki百科）：

GeoHash lengthlat bitslng bitslat errorlng errorkm error
1 2 3 ±23 ±23 ±2500
2 5 5 ±2.8 ±5.6 ±630
3 7 8 ±0.70 ±0.70 ±78
4 10 10 ±0.087 ±0.18 ±20
5 12 13 ±0.022 ±0.022 ±2.4
6 15 15 ±0.0027 ±0.0055 ±0.61
7 17 18 ±0.00068 ±0.00068 ±0.076
8 20 20 ±0.000085 ±0.00017 ±0.019

base32 編碼對映表：

ValueSymbolValueSymbolValueSymbolValueSymbol
0 A 9 J 18 S 27 3
1 B 10 K 19 T 28 4
2 C 11 L 20 U 29 5
3 D 12 M 21 V 30 6
4 E 13 N 22 W 31 7
5 F 14 O 23 X    
6 G 15 P 24 Y    
7 H 16 Q 25 Z    
8 I 17 R 26 2    

相關文章