熱気球のGPSログから風の高度分布を推定

私は現在、[熱気球のパイロット](http://www.jballoon.jp/balloon/pilot.html)を目指してトレーニング中です。1ヶ月に1回ほど北関東の[渡良瀬遊水地](http://www.jballoon.jp/flightarea/kanto.html#watarase)周辺で指導を受けながら飛行しています。

大学では気象の研究をしていたこともあり、「ただ飛ぶだけでは勿体ない」と思い**上空の大気のデータ**を取りたくなりました。今回は手始めに、手元のスマートフォンのGPSで記録した飛行データから風の高度分布を出してみます。

とはいえプログラミング的には難しいことはせず、Rubyで単純にXMLを読み込んで数値を取り出すだけです。

> [ラジオゾンデによる高層気象観測](https://www.jma.go.jp/jma/kishou/know/upper/kaisetsu.html)でも現在はGPSを用いて風を測定しています。

## 熱気球の飛行の仕組み
（GPSログから風を出す原理はとても単純なのですが、一応熱気球の飛行の仕組みが関わるので纏めておきます。）

風船みたいに浮くだけの熱気球が自在に移動できるのかよく質問されます。その答えを知るには、上下方向と水平方向に分けると理解しやすいです。

簡単な紹介に留めますので、専門的に知りたい場合は日本気球連盟の[パイロットハンドブック](http://www.jballoon.jp/safety/handbook/handbook.html)（3.1. 熱気球の飛行原理, 3.2. 熱気球の構造, 3.4. 航法）を参照してください。

### 上下方向
熱気球のメイン部分である、<ruby>球皮<rp>（</rp><rt>きゅうひ</rt><rp>）</rp></ruby>と呼ばれる大きな袋に入れた**空気の浮力[^buoyancy]**で上下に移動します。熱気球という名前通り、**バーナーで空気を加熱する**ことで中の空気を軽くし浮力を制御します。時間が経ち空気が冷めることで自然と下降しますが、ほかに天頂部にある排気機構を使って**熱気を逃がす**ことでより早く下降することもできます。

ちなみにトレーニングを始めて体感していますが、一定高度を維持するのは結構難しいです。バーナーは基本的に火力調節できずON/OFFだけなうえ、焚いてから浮力変化の反応を感じるまで*10秒程度のタイムラグ*があります。また同じ量だけ焚いたとしても、総重量や外気温などによって反応の大きさは変わります。もし自動車の速度調整が同じような性質だったら、車間距離をうまく維持できず追突必至でしょう。

[^buoyancy]: 浮力とは、[アルキメデスの原理](https://ja.wikipedia.org/wiki/アルキメデスの原理)では「物体が押しのけた流体の重さに等しい上向きの力」を指します（物体の重さには依存しません）が、ここでは物体の重さを加味した「外気と内気の密度差による正味の力」という意味で使用しています。

### 水平方向
動力や舵などが無いため**【風まかせ】**です。風に乗る、つまり周囲の空気と共に流れていくのみです。飛行船やグライダーなどのように能動的には動けません。

ただし、風は高度によって方向や速さが異なるため、**上下に移動してうまい風を選ぶ**ことで目的の方向へ飛ぶことができます。そのため、どの高度にどんな風があるか把握することが重要となります。場所や時間によっても変わるのが難しいところです。

風次第なので全方位自由自在とはいきません。そのため熱気球で飛ぶ前には気象データや測風気球（パイバル）で上空の風を調べ、離着陸地点やコースなど飛行計画を立てます。[競技飛行](https://www.honda.co.jp/balloon/library/competition/index.html)の場合もスタッフが直前に調べて、種目（タスク）や目標地点を無理のないものに決めます。また飛行中は他の熱気球の動きがとても参考になります。風の状況が変わってきたら飛行計画も柔軟に修正する必要があります。

---

風に乗るということは、単純には**熱気球の水平移動速度がそのまま水平風の速度**ということになります。そのためGPSから位置と速度を出せばその地点の風がだいたい分かります。

ラジオゾンデ観測と異なる点は主に、熱気球は上下方向に自由に動けることと、観測用ゴム気球と比べて大きく重いことです。これらによって測定に何かしら影響があるかもしれません。

<!--
> もしかしたら熱気球の速さと風速との間には系統誤差があるかもしれません。その場合は風速を正しく測れないことになりますが、熱気球の飛行を助けるデータとしては依然有用です。
-->

## サンプルデータ
トレーニングで解析にちょうど良さそうな飛行ができたため、そのときのGPSログを使います。

### GPXファイル
[GPSLogger](https://play.google.com/store/apps/details?id=com.mendhak.gpslogger) で記録したGPXファイルは以下のようなデータになっています。中身はXMLです。（見やすく整形しました）

```xml
<?xml version="1.0" encoding="UTF-8"?>
<gpx version="1.0" creator="GPSLogger 99 - http://gpslogger.mendhak.com/"
     xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.topografix.com/GPX/1/0"
     xsi:schemaLocation="http://www.topografix.com/GPX/1/0 http://www.topografix.com/GPX/1/0/gpx.xsd">
	<time>2019-02-23T22:10:31.000Z</time>
	<trk>
		<trkseg>
			<trkpt lat="36.25550486" lon="139.65405245">
				<ele>18.0</ele>
				<time>2019-02-23T22:10:31.000Z</time>
				<speed>0.0</speed>
				<geoidheight>39.0</geoidheight>
				<src>gps</src>
				<sat>13</sat>
				<hdop>0.9</hdop>
				<vdop>0.8</vdop>
				<pdop>1.2</pdop>
			</trkpt>
			<trkpt lat="36.25550062" lon="139.65405592">
				<ele>19.0</ele>
				<time>2019-02-23T22:10:33.000Z</time>
				<course>186.5</course>
				<speed>0.28</speed>
				<geoidheight>39.0</geoidheight>
				<src>gps</src>
				<sat>13</sat>
				<hdop>0.9</hdop>
				<vdop>0.8</vdop>
				<pdop>1.2</pdop>
			</trkpt>
			<trkpt lat="36.25549219" lon="139.65406112">
				<ele>19.0</ele>
				<time>2019-02-23T22:10:35.000Z</time>
				<course>192.7</course>
				<speed>0.5</speed>
				<geoidheight>39.0</geoidheight>
				<src>gps</src>
				<sat>13</sat>
				<hdop>0.9</hdop>
				<vdop>0.8</vdop>
				<pdop>1.2</pdop>
			</trkpt>
...
		</trkseg>
	</trk>
</gpx>
```

速度は2点間の差分をとって計算しようかと思っていましたが、**速さと方位**のデータが入っていたためそのまま使えそうです。これが電波のドップラー偏移で測っているなら高精度が期待できる[^accuracy]ものの、スマートフォンのGPSにそこまでの機能があるのかはよく分かりません。

[^accuracy]: 数社のGPSゾンデの製品仕様を見てみたところ、不確かさは風速 0.15m/s ・風向 1～2° となっていました。（受信機の性能のみでなく、観測システムとしての総合的な値と思われます）

### 航跡
以下の図はログを [GPX Viewer](https://play.google.com/store/apps/details?id=com.vecturagames.android.app.gpxviewer) で可視化したものです。基本的に南向きに移動しているのですが、途中で東や西に振れています。これは飛行高度を変えた影響が大きいです（下部のグラフ参照）。

<!-- ![gpx-viewer-0.png](https://qiita-image-store.s3.amazonaws.com/0/127130/0d99f280-c3b1-9fc2-7ff0-066b15c7e12a.png) -->
![gpx-viewer-1-95.png](https://qiita-image-store.s3.amazonaws.com/0/127130/8b872b5f-a7cc-bbc4-65c7-72f9927a5b02.png)

風の解析では、図に描き込んだ**3つの飛行区間**に注目していきます。

0. 高度 250m での高度維持 → 同一高度での速度記録は一定になるか
0. 高度 1000m までの上昇 → 途中の高度での速度分布を得られるか
0. 高度 200m までの急降下 → 上昇時と似た速度分布を得られるか

## データ処理
### XMLの読み込み
ひとまずデータを読み込み、各地点の**経度・緯度・高度・時刻・東西速度・南北速度**の数値を取り出します（速さと方位から、東向き・北向きの速度成分を計算します）。今回のデータはトラックセグメントがひとつだけと分かっているので、全てのトラックポイントをひとつの配列に入れます。

コードは[Rubyリファレンスマニュアルの例](https://docs.ruby-lang.org/ja/latest/library/rexml=2fdocument.html)を真似しています。

```ruby:GPXファイルの読み込み
require "rexml/document"
require "time"

TrackPoint = Struct.new(:lon, :lat, :alt, :time, :u, :v)

doc = REXML::Document.new(File.read(ARGV[0]))

track_points = REXML::XPath.match(doc, "/gpx/trk/trkseg/trkpt").map do |trkpt|
	lon    = trkpt.attribute("lon").value.to_f
	lat    = trkpt.attribute("lat").value.to_f
	alt    = trkpt.elements["ele"].text.to_f
	time   = Time.iso8601(trkpt.elements["time"].text).localtime
	speed  = trkpt.elements["speed"].text.to_f
	course = trkpt.elements["course"]&.text.to_f # speedが0のときは要素が無い → 0と扱う

	v, u = Complex.polar(speed, course * Math::PI / 180).rect
	TrackPoint.new(lon, lat, alt, time, u, v)
end

# この後に計算処理や出力を追加する
```

※経度・緯度は属性読み込みのサンプルなだけで、以降では利用しません。

#### サンプリング間隔の確認
念のため、データが密に取れているか確認します。

```ruby:サンプリング間隔の度数分布
intervals = track_points.each_cons(2).map { |tp1,tp2| tp2.time - tp1.time }
histogram = intervals.each_with_object(Hash.new(0)) { |dt,hash| hash[dt.round] += 1 }

p histogram.sort.to_h
#=> {1=>349, 2=>927, 3=>11, 4=>3, 5=>2, 6=>2, 7=>2, 8=>1, 9=>1, 11=>1, 12=>1, 18=>1, 26=>1, 29=>1, 70=>1}
```

ほとんどのデータは1〜2秒間隔ですが、5秒以上開いている地点も1%程度あります。 GPSLogger の設定（付録参照）で「測位精度が5m以下になるまで測り直す」ようにしていたため、なかなか精度が出なかったときがあったのでしょう。

### データそのまま
熱気球は風に乗るので、**熱気球の移動速度はそのまま風の速度を反映している**と考えられます。なのでまずは高度と速度のデータをそのままプロットしてみます。

```ruby:TSV形式で書き出し（生データ）
track_points.each do |tp|
	puts [tp.time.strftime("%T.%3N"), tp.alt, tp.u, tp.v].join("\t")
end
```

```text:stdout
07:10:31.000	18.0	0.0	0.0
07:10:33.000	19.0	-0.03169689985501386	-0.2782001195894445
07:10:35.000	19.0	-0.10992310217641847	-0.48776727197292835
...
```

![gps-wind-0.png](https://qiita-image-store.s3.amazonaws.com/0/127130/df6ac28e-e88f-3556-8a65-be070b9ae3d6.png)

図の左2つは風の鉛直プロファイル（東西風・南北風）、右は風のホドグラフです。単位は風速が m/s 、高度が m です。サンプリングが5秒以上開いた箇所は線を結んでいません。

* 高度維持時（区間1：赤線）では、**速度のばらつきが ±0.5m/s 程度**と非常に安定しています。ブレの原因はGPSの精度と風の変動の両方が考えられますが、いずれにせよこれが今回の測定で出せる精度と考えていいでしょう。
* 上昇時（区間2：太線）では、線が細かくブレながらも連続的に速度変化を捉えています。ホドグラフからは、南東〜西北西へ行く様々な風や無風に近い層があることがわかります。
* 急降下時（区間3：灰色の細線）では、鉛直プロファイルで見ると上昇時と符号は似ているもののかなり波打っています。ホドグラフではもはや軌跡が分かりません（ので色を薄くしました）。

この結果を踏まえ、上昇・急降下時（区間2～3）の航跡を拡大して再確認してみます。

![gpx-viewer-11.png](https://qiita-image-store.s3.amazonaws.com/0/127130/ec9e2685-b5cc-6415-6522-bdd00ef44cc7.png)

* 上昇時に経路が直角に曲がったり反転したりと、ホドグラフの通り様々な方向へ進んでいます。GPSログの位置と速度に不整合は無さそうです。
* 急降下時は経路がギザギザになっています。熱気球の動きが乱れたのであれば、測定した速度に風以外の影響が含まれているかもしれません。

### 空気抵抗を考慮
今回の実験をするにあたって気になっていたことがあります。「熱気球の速度は*常に*風速と同じ」と考えていいのかということです。

熱気球は風に乗るとはいえ、慣性があるため風の変化に一瞬で追従するとは考えにくく、**何秒か時間をかけて風速に近づいていく**はずです。同じ高度を飛び続けるなら気にしなくていいものの、今回の例のように**上昇中・下降中にすぐ通過する層の風を知る**には補正が必要ではないかと思いました[^feeling]。

空気抵抗の単純なモデルとして、「速度差に比例した抵抗力がかかる」と考えてみます。比例定数 $1/\tau$ を用いると、運動方程式より次の式ができます。

```math
\frac{d\vec{v}}{dt} = -\frac{1}{\tau} \left( \vec{v} - \vec{V} \right)
```

> 例えば風速 $\vec{V}(t)$ が一定値 $\vec{V}_0$ という条件でこの微分方程式を解くと、熱気球の速度 $\vec{v}(t)$ は以下の式になります。
> 
> ```math
> \vec{v}(t) = \vec{V}_0 + \left( \vec{v}(0) - \vec{V}_0 \right) e^{-t/\tau}
> ```
> 
> 右辺第2項は、時間 $\tau$ 経つごとに熱気球と風の速度差が **$1/e$ 倍**（約0.37倍）に減っていくことを表しています（**緩和時間**）。そのため十分に時間が経つと $\vec{v}(t) \to \vec{V}_0$ 、つまり風と同じ速度となります。

今回は $\vec{v}$ が分かっていて $\vec{V}$ を求めたいので、単純に以下の式で計算することになります。熱気球の速度に、その微分（加速度）を考慮した項が加わりました。この項が風の変化に対する熱気球の応答の遅れを補います。

```math
\vec{V} = \vec{v} + \tau \frac{d\vec{v}}{dt}
```

緩和時間 $\tau$ は熱気球の質量などで変わるはずですが、今回は推定方法が思いつかなかったので試しに 2秒 にしておきます。また、データの時間間隔が開いているところは微分が他より平均的な値になってしまうので、欠測ということにします。

```ruby:TSV形式で書き出し（補正あり）
tau = ARGV[1].to_f

track_points.each_cons(2) do |tp1,tp2|
	dt = tp2.time - tp1.time
	dh = tp2.alt  - tp1.alt
	du = tp2.u    - tp1.u
	dv = tp2.v    - tp1.v

	if dt >= 5 # 欠測扱い
		t = tp1.time + dt/2
		h = u = v = -999.0
	else
		t = tp1.time + dt/2
		h = tp1.alt  + dh/2
		u = tp1.u    + du/2 + tau * du/dt
		v = tp1.v    + dv/2 + tau * dv/dt
	end

	puts [t.strftime("%T.%3N"), h, u, v].join("\t")
end
```

![gps-wind-2.png](https://qiita-image-store.s3.amazonaws.com/0/127130/661bac00-c687-eb20-7fed-01884abf16ab.png)

生データを見た段階である程度予想はしていましたが、微分によってブレが増幅されてよりギザギザしたグラフになってしまいました。

高度分布を見てみると、生データと比べて風速はほとんど変わらないものの高度が 数m 変わっています。これで上昇・下降時のずれが改善されれば良かったのですが、元々のグラフがそこまで似ていないため何とも言えません。そもそも場所も時刻も異なる比較なので、重なることが正しいとは限りません。

空気抵抗の考慮は風の測定の補正に使える可能性はありそうですが、効果は高度が僅かに変わるだけで、データのブレを抑えないと却って邪魔かもしれません。ローパスフィルターで平滑化する手もあるものの、データの位相が遅れたら補正の意味が無いため悩み中です。

[^feeling]: 実際熱気球に乗っていると、上昇・下降時は多少の横風を感じます。

## まとめ
熱気球のGPSログから高度と速度の情報を取り出すことで、各高度での風を推定してみました。

* 一定高度を維持したときは速度のばらつきが ±0.5m/s 程度と非常に安定していました。
* 上昇時のデータから、高度によって様々な風があることを捉えられました。
* 急降下時は熱気球の位置自体が複雑に変化していて、風の情報をあまり捉えられなかったかもしれません。
* 熱気球の空気抵抗を考慮して補正するには、もっとデータのブレを抑える必要がありそうです。今回はあまり検証できませんでした。

### 今後したいこと
家に帰ってからログを解析するという方法では、推定した風が実際にあるか確かめることがなかなかできません。そこで、簡単なスマホアプリを作って**リアルタイムで風の高度分布を推定・表示**し、飛行中その場で検証することに挑戦したいです。いずれはセンサーで気温や気圧などのデータも取り、より良く大気の状態を把握して飛行に役立てられればと思っています。

とはいえまだパイロット訓練生の身なので、身につけるべき知識・技術はたくさんありますし、アプリに頼りきりでなく感覚で風をつかめるようにトレーニングしていきたいです。

## 付録
### GPSLogger パフォーマンス設定
GPSLogger は記録方法を細かく設定できる反面、状況に合わない設定をしてしまうとうまく記録されません。次回以降の参考にできるよう、今回の設定を記載しておきます。なるべく記録が途切れないことを重視し、ログ肥大化や電池消費を許容した設定です。

* スマートフォン機種：Xperia Z3 (Android 5.0.2)
* GPSLogger バージョン：99

| 項目 | 値 | 備考 |
|:-:|:-:|:-:|
| 位置情報プロバイダ | GPS のみ |   |
| ロギング・インターバル | 1 秒 | 0秒だと記録し過ぎる |
| GPSを常にON | yes | 測位精度をリセットさせない |
| 距離フィルター | 0 m | 位置が変わらなくても記録させる |
| 精度フィルター | 5 m | 空が開けているので達成しやすい |
| 精度にマッチする<br>までの待ち時間 | 60 秒 |   |
| GPS信号受信までの<br>最大待ち時間 | 30 秒 |   |
| WGS84の代わりに<br>MSLを使用 | yes | 海抜高度で記録させる<br>（参考：[ジオイドとは](http://www.gsi.go.jp/buturisokuchi/geoid.html)） |
| 高度からオフセットを<br>差し引く | 0 m | 〃 |
| 移動していないときは<br>ログファイルを書かない | no | 以前ログが飛んだためnoに変更 |