pythonで電波伝播モデルから発信機の位置を計算してみる[Wi-Fi, Beacon編]

8523 ワード
location Python やってみた beacon WiFi Python テキストリンク

※このBlogは弊社データワイズの技術ブログはじめてみた、のいちエントリになります

概要

目的:
一般的にWi-FiやBeaconを用いた位置計算は裏で何をやっているのかをなんとなくわかった風になって理解してもらう。

方法:
基本的な考え方を説明して、Pythonコードを使って実演してみる。

環境:
MacOS Mojave 10.14.6, Python 3.7.4

基本的な考え方

必要なもの

  1. 電波の受信機(3個)とその位置
  2. 発信機(スマホとか)とその位置
  3. 各受信機に届いた電波の強度
  4. 奥村-秦カーブ

手順

Step1. 発信機から受信した電波強度と受信機の位置情報から奥村-秦カーブで距離を計算
Step2. 計算された距離を半径、各受信機の位置を中心とした円を用意
Step3. 三点測位で円の交点を求める。交点の座標=発信機の位置座標

三点測位の具体的なイメージはこのエントリの一番下にあります。
興味ある方は先にそちらをご覧ください。

コードと結果

Step0. 正解データの準備

実際に計算を行うための正解データを準備します。
今回は弊社がある虎の門ヒルズとその周辺施設を対象としました。
具体的な位置はこんな感じ。

これらのlatitude/longitudeを取得したので、
三平方の定理から距離を、
距離から各受信機における(理論上の)電波強度を求めます。

今回はHata-model:WikipediaのUrban environments, small or middle size cityを採用してます

# -*- coding: utf-8 -*-
import math
import numpy as np
import pandas as pd
import sympy.geometry as sg

# latitude 1 digree = 0.0090133729745762 km
# longitutde 1 degree = 0.010966404715491394

def eucl_dist(pointA, pointB):
    return np.sqrt((pointA[0]-pointB[0])**2+(pointA[1]-pointB[1])**2)
    # euclidean distance

def Oku_Hata(f, h_b, C_H, d):
  Loss = (69.55 + 26.16*math.log10(f)
              - 13.82*math.log10(h_b) - C_H 
              + (44.9 - 6.55*math.log10(h_b))*math.log10(d))
  return Loss
  # Okumura_Hata Curve

def inv_Oku_Hata(f, h_b, C_H, Loss):
    numerator = Loss - 69.55 - 26.16*math.log10(f) + 13.82*math.log10(h_b) + C_H
    denominator = 44.9 - 6.55 * math.log10(h_b)
    d = 10**(numerator/denominator)
    return round(d, 10)
  # Inversed Okumura_Hata Curve

# variables 
f = 150  # frequency of transmissiton: Megahertz
h_b = 1  # height of base station antenna : meter
h_M =  1  # heigth of mobile station antenna ; meter
C_H =0.8 + (1.1 *math.log10(f) - 0.7) * h_M -1.56*math.log10(f)
d = 500

# Test whether or not the functions function?
euDist = eucl_dist([0,0],[3,4])
Loss = Oku_Hata(f, h_b, C_H, d)
d = inv_Oku_Hata(f, h_b, C_H, Loss) 

print('euDist should be 5: ', euDist)
print('Loss from distance :  ', Loss)
print('d should be 500: ', d)  # Should be 500

検算結果

euDist should be 5: 5.0
Loss from distance : 248.5613025107495
d should be 500: 500.0

式はちゃんとできているようです。
次に実際の位置情報から距離を計算し、Lossに直します。

# Set of lattitude/ longitude each GeoData
tmitter = [35.6662344,139.7496597]  # Transmitter from Torano-Mon

sensorLocation = {'S1':['GoodMorningCafe',35.6664959,139.7500406],
                  'S2':['LeParc',35.6666179,139.7479538],
                  'S3':['GrandSuite',35.6675194,139.7489692]}
                  # SensorLocation around Trano-Mon

# Compute distance (= radius) from the transmitter
r1 = eucl_dist(tmitter, sensorLocation['S1'][1:])
r2 = eucl_dist(tmitter, sensorLocation['S2'][1:])
r3 = eucl_dist(tmitter, sensorLocation['S3'][1:])
radii = [r1, r2, r3]  # plural of radius

# Compute Loss
L1 = Oku_Hata(f, h_b, C_H, r1)
L2 = Oku_Hata(f, h_b, C_H, r2)
L3 = Oku_Hata(f, h_b, C_H, r3)

print('List of distance: ', radii)
print('List of Loss: ', L1, L2, L3)

結果

List of distance: [0.0004620249560447818, 0.0017484756389397347, 0.0014587718293119327]
List of Loss: -22.378972679965557 3.572964717331658 0.040582137429893805

Step1. 発信機から受信した電波強度と受信機の位置情報から奥村-秦カーブで距離を計算

Lossが(距離から)計算できたので、
距離を(Lossから)計算します。
正解データでやってるんだから、回りくどくても仕方ない。最初から距離使えよとか言わないで！

d1 = inv_Oku_Hata(f, h_b, C_H, L1)
d2 = inv_Oku_Hata(f, h_b, C_H, L2)
d3 = inv_Oku_Hata(f, h_b, C_H, L3)

print('Computed Value: ', [d1, d2, d3])  # plural of radius
print('Original Value: ', radii)

検算結果

Computed Value: [0.000462025, 0.0017484756, 0.0014587718]
Original Value: [0.0004620249560447818, 0.0017484756389397347, 0.0014587718293119327]

式の途中で値を丸めているので桁数は違いますがちゃんと計算できているようです。

Step2. 計算された距離を半径、各受信機の位置を中心とした円を用意

コード。sympyを使うのが楽そうなので使いました。
Qiita:sympyで円と円の交点

# Step2.
centers = [sg.Point(sensorLocation['S1'][2], sensorLocation['S1'][1]),
           sg.Point(sensorLocation['S2'][2], sensorLocation['S2'][1]),
           sg.Point(sensorLocation['S3'][2], sensorLocation['S3'][1])]
           # x = longitude, y = latitude

circles = [sg.Circle(centers[0], radii[0]), 
           sg.Circle(centers[1], radii[1]),
           sg.Circle(centers[2], radii[2])]
           # define all three circles

print('List of centers (latitude/longitude) :', centers)
print('List of Circles :', circles)

結果

中心の座標

List of centers (latitude/longitude) : [Point2D(698750203/5000000, 356664959/10000000), Point2D(698739769/5000000, 356666179/10000000), Point2D(349372423/2500000, 178337597/5000000)]

円の情報。半径が追加されてます。中心座標と半径のリストですね。

List of Circles : [Circle(Point2D(698750203/5000000, 356664959/10000000), 231012478022391/500000000000000000), Circle(Point2D(698739769/5000000, 356666179/10000000), 174847563893973/100000000000000000), Circle(Point2D(349372423/2500000, 178337597/5000000), 145877182931193/100000000000000000)]

Step3. 三点測位で円の交点を求める。交点の座標=発信機の位置座標

コード。やっぱり sympy (シンパイ？)を使います。
３つの円がありますが、２つずつ交点を計算しているので6つの解が得られます。
なので、6つの解のうち、3回出ている交点を抽出するロジックにしています。
その為にroundして計算の誤差を無視してます。

# Step3. 
crossPoints12 = sg.intersection(circles[0], circles[1])
crossPoints23 = sg.intersection(circles[1], circles[2])
crossPoints31 = sg.intersection(circles[2], circles[0])
# Compute crosspoints for each two circles

crossPointsList = [
    [float(crossPoints12[i].y),float(crossPoints12[i].x)] for i in [0,1]
    ] + [
    [float(crossPoints23[i].y),float(crossPoints23[i].x)] for i in [0,1]
    ] + [
    [float(crossPoints31[i].y),float(crossPoints31[i].x)] for i in [0,1]
    ]
crossPointsList = [[round(i[0],8), round(i[1],8)]for i in crossPointsList]
df = pd.DataFrame(crossPointsList).groupby([0, 1])[0].count()
print('Computed Value: ', df[df.values == max(df.values)].index[0])
print('Original Value: ', tmitter)
# print the most frequent set of lattitude-longitude

結果

Computed Value: (35.6662344, 139.7496597)
Original Value: [35.6662344, 139.7496597]

正しく正解が導けたようです。めでたしめでたし。

補足：三点測位のイメージ、など

"三点測位"でググるといっぱい出てくるんですが、とりあえずここのブログがざっとみた中で一番古いエントリだったのでリスペクト紹介。
参考：iBeaconのRSSIでiPhoneの二次元座標をとれたらいいな

で、本末転倒な話をすると、上記のブログ以外のブログでも話していることですが、
実務としてはこのコード、特にHataカーブのところが使い物になりません。

都会は物が乱立してて電波が乱反射する為、電波伝播が理想の状態で受信機に届かない為です。
その辺の詳しい話は技術ブログの域を超えてしまうので紹介する予定はありませんが、
知りたい方はググるか弊社にジョインをしていただけるといいのかなと思います。

以上です。

どうでもいいけど、デンパデンパって舌噛みそうになる