標高タイルで電波の見通しを判定する

地図の上で自宅と中継局を直線で結び、間の標高と比べれば、電波が山にさえぎられるかどうかは分かるはずです。地デジの受信環境を調べるアプリを作り始めたとき、見通し判定はその程度の作業だと見積もっていました。この見積もりは二度崩れます。一度めは標高データの入手と読み方に、二度めは地球の丸みにです。この記事では、国土地理院の標高タイルから地形断面を作り、中継局への見通し(山かげ)を判定した実装を記録します。

PNG画像として届く標高データ

国土地理院は日本全域の標高を標高タイルとして配信しています。地図タイルと同じXYZ形式で、ズームとタイル座標をURLに埋めるとタイルが1枚返ってきます。このアプリでは10mメッシュ由来のDEM10Bというデータを、ズーム12(緯度35度付近で1ピクセル約31m、1タイルで約8km四方)で使いました。なお標高タイルは地理院タイルの一種なので、アプリの画面には出典として地理院タイルを利用している旨を表示しています。

最初のつまずきは、返ってくるものがPNG画像だという点でした。数値の配列を期待して開くと画像で、しかもこの画像は目で見るための地形図ではありません。各ピクセルのRGB値そのものが標高の符号化で、r、g、bを上位から並べた24ビット整数を組み立て、0.01m単位の標高として読みます。

const DEM_NO_DATA = 8388608; // 2^23

export const decodeDemPixel = (r: number, g: number, b: number): number | null => {
  const x = r * 65536 + g * 256 + b;
  if (x === DEM_NO_DATA) return null;
  if (x > DEM_NO_DATA) return (x - 16777216) * 0.01;
  return x * 0.01;
};

たとえばRGBが(0, 247, 168)のピクセルはx = 63400なので標高634.00mです。2の23乗はデータなしを表す無効値としてnullを返し、それを超える値は海面下の標高(xから2の24乗を引いて0.01倍した負値)を表します。ピクセル値の取り出しは、PNGを自前で解析する代わりに、Imageとして読み込んでcanvasに描き、getImageDataでRGB配列を得る形にしました(タイルは別オリジンなのでcrossOriginをanonymousに指定します)。

標高断面のサンプリング

タイルさえ読めれば、任意の緯度経度の標高が引けます。緯度経度をワールドピクセル座標に変換し、タイル座標とタイル内ピクセルに割って、該当ピクセルをデコードするだけです。問題は、測定点から中継局までの断面を何点で作るかでした。

実装では2点間を200m間隔を目安に等分し、点数は最小16、最大400に収めています。中継局までの距離は数kmから100km超まで幅があるので、間隔だけを固定すると近距離で粗すぎ、遠距離では点数が際限なく増えるためです。各点の緯度経度は線形補間で求めています。厳密な大円ではなく航程線の近似ですが、120km以下の径路なら大円との横ズレは200〜300m程度で、サンプル間隔と同じオーダーに収まるため、参考判定としては許容しました。

デコードで無効値になった点(海上や欠測)は0m扱いにしています。海を越える径路では海面が断面に入るだけなので、見通し判定の用途にはこの割り切りで足ります。

地球の丸みと等価地球半径

断面ができれば、あとは両端を結ぶ直線と各点の標高を比べるだけ、のはずでした。実際、10km程度の近距離ならそれでほぼ済みます。ところが100km先の局に対しては、間に山が一つもなくても視線が通らない場合があります。地球が丸いからです。

長い径路では、中間の地面は両端を結ぶ直線より持ち上がって見えます。持ち上がりの量は、両端からの距離をd1、d2、地球半径をRとしてd1 × d2 / (2R)で近似できます。

ここでRに実際の地球半径6371kmをそのまま使わないのが、電波の見通し計算の慣行です。大気の屈折率は高さとともにわずかに下がるため、電波は直進せず、地面に沿う向きにゆるく曲がりながら進みます。曲がる電波と丸い地球をそのまま扱う代わりに、電波を直線と見なし、地球の半径を4/3倍に膨らませて帳尻を合わせる近似が使われます。この膨らませた半径が等価地球半径で、光で見た見通しより電波の見通しが少し遠くまで届くことを表しています。

const EFFECTIVE_EARTH_RADIUS_M = (4 / 3) * 6371000;

// 地球の丸みによる持ち上がり[m](等価地球半径 k=4/3)
export const earthBulgeM = (d1M: number, d2M: number): number =>
  (d1M * d2M) / (2 * EFFECTIVE_EARTH_RADIUS_M);

持ち上がりは直感よりずっと大きい量です。全長120kmの径路の中点(両端から60km)で約212m、全長100kmの中点でも約147mになります。平坦な土地で両端のアンテナ高が10mなら、山がなくても地球の丸みだけで視線はさえぎられます。単体テストにもこのケースを入れて、湾曲の項が効いていることを確かめています。

視線と断面の比較

断面と持ち上がりがそろえば、あとは視線と比べるだけです。受信側の高さは測定点の標高に10mを足します(家屋の屋根上アンテナの想定)。送信側は局の海抜高データがあればその値を使い、なければ送信点の標高に50mを足した概算にして、概算であることをフラグで持ち回ります。視線の高さは受信側から送信側へ距離の比で線形補間し、各サンプル点で「標高 + 持ち上がり − 視線の高さ」を計算します。これが正なら、その地点の地形が視線をさえぎっています。

const t = point.distanceKm / totalDistanceKm;
const rayHeightM = rxHeightM + (txHeightM - rxHeightM) * t;
const bulge = earthBulgeM(point.distanceKm * 1000, (totalDistanceKm - point.distanceKm) * 1000);
const obstruction = point.elevationM + bulge - rayHeightM;

この値の最大値と、それが出た地点の距離を記録しておき、最大値が正なら遮蔽あり、負なら見通しありです。負のときの最大値は、径路上で最も余裕のない地点のマージンとしてそのまま読めます。

一つだけ例外があります。両端から0.5km以内のサンプルは判定から除外しました。測定点の直下や、アンテナ鉄塔が立つ山頂の標高は視線とほぼ同じ高さから始まるため自己交差しやすく、そこまで真面目に判定すると、山頂に立つ送信局が軒並み「遮蔽」になってしまうのです。除外の結果サンプルが一つも残らないほどの近距離は、見通し扱いにしています。

タイルキャッシュと失敗時の再試行

断面1本の判定で、標高の参照は最大400回、タイルにして最大20数枚に及びます。地図をタップし直すたびに同じタイルを取り直すわけにはいかないので、デコード済みのピクセル配列をタイル座標をキーにキャッシュしています。上限は24枚(120km径路1本分)で、参照したタイルをMapの末尾に差し直し、あふれたら先頭の最も古いものから捨てる、Mapの挿入順を利用した小さなLRUです。

キャッシュに入れているのは値ではなくPromiseです。同じタイルへの参照が並行して走っても、実際の取得が1回にまとまります。ただしこの形には罠があって、失敗したPromiseをそのまま残すと、一度の通信エラーが「このタイルは永遠にエラー」として固定されてしまいます。実装では失敗時にキャッシュから取り除き、次のアクセスで再取得させることにしました。判定結果そのもののキャッシュ(局と丸めた座標がキー)も同じ方針で、成功した結果だけを残し、エラーは残しません。

判定結果の表示と判定の限界

画面では、選択中の局について「地形見通し: 良好(地形遮蔽なし)」、または「地形遮蔽の可能性: 測定点から約◯km地点で約◯mの遮蔽」のように、遮蔽の位置と量まで表示しています。どのあたりの山が影を作っているかまで分かると、受信不調の説明として納得しやすいからです。送信点高が概算のときは「※送信点高不明のため概算」を添え、標高データを取得できなかったときは判定不可と表示します。

ただし、この判定が言えるのは地形の話だけです。標高タイルが表すのは地表の高さで、建物や樹木は含まれません。隣のビルや裏山の木立が作る影は、この判定には現れません。また視線1本の幾何判定であり、フレネルゾーンの確保や山岳回折は扱っていないので、「良好」でも受信できる保証にはならず、「遮蔽の可能性」と出た場所でも回折などで受信できていることはあります。アプリ内でもこの判定は、受信可否の答えではなく、受信不調の原因が山かげかどうかのあたりを付ける参考情報として位置づけ、その旨を注記しています。

まとめ

「地図上の2点を直線で結べば見通しが分かる」という最初の見積もりは、それ自体は間違っていませんでした。ただ、直線と比べる相手の断面は、PNGのRGB値から標高を復元し、200m間隔で並べて初めて手に入ります。比べる基準の側にも、地球の丸みと電波の屈折(等価地球半径)の補正が要ります。式にすればどれも数行で、その数行に必要な事実を集めることが実装の大半でした。

この見通し判定は「地デジエリアマップ」で動いています。同じアプリの方位センサーまわりは「スマホのコンパス方位をWebで取得する」に書きました。地形の次にこの判定を狂わせるとしたら、標高タイルに写らない建物の影です。そこまで扱うかどうかは、まだ決めていません。