Webカメラ＋火星でピリオデックモーション測定

シーイングの悪い夜、大気の乱れが収まる兆候はないかと期待しながらWebカメラ越しに火星を眺めていたら、ゆらゆら揺れている火星がスーッと視野の端に近づいた。「極軸合わせをしくじったかな」と思って見ていると、3分位で今度はスーッと戻ってきた。おやおや、これは間違いなくPE（Periodic motion Error、ピリオデックモーション）だ。極軸がずれている場合は行ったきりで、こんなに短期間で戻ってくるはずは無い。

シーイングが改善しそうも無いので今夜はPE測定にしよう、と言うことで、カメラの向きを出来るだけ正確に上=北、左=東に合わせ、ウォーム2回転分、20分間の動画を撮った。先ずは直感的に把握するため、30倍のコマ落とし、今風に言うTime-lapse Videoを作ってみた。

目見当だと直径数個分かな。これは情緒に訴えるには良さそうだが、定量的に測るのには不向きだ。最後のほうで暗くなっているのは雲がかかったせいである。

今回はこの動画を元に出来るだけ正確にPEを測ってみた。また主な原因がわかったので、その対策も併せて紹介しよう。

画像などデータ処理には、動画の変換にffmpegを使ったほかはプログラミング言語Pythonで短いコードを書いた。Pythonは道具として使うには最も適したプログラミング言語の一つで、覚えておいて損はないと思う。これは各種ライブラリーモジュールが充実していることが大きく寄与している。今回使った主要なライブラリーモジュールは以下の通りだ。

• 数値処理：numpy

• 視覚化：matplotlib（SciPy/PyLab）

• 動画表示：PyGame

• 動画・静止画キャプチャーなど：OpenCV、VideoCapture

• 静止画の加工：PIL（Python Image Library）

なおPythonやそのライブラリーモジュールについては、解説やドキュメント類がインターネット上に豊富にあるので、そちらを参照してほしい。

/*******  *******/

伝統的方法のエミュレーション

伝統的PE測定手法は、極軸をわざと東西どちらかにずらして、南中前後の天の赤道付近の恒星を長時間露光するものだ。こうすると、PEが無ければ概ね南北一直線に流れて写るはずの星像が東西方向にくねくねする振れ幅でPEの量を測ることができる。視野の中で角距離のわかっている恒星のペアを角度の基準にする。例えばこんな具合；
http://katabami.blog.so-net.ne.jp/2007-01-08
http://ponkotu-gg.cocolog-nifty.com/blog/2008/08/post_ecde.html

これをエミュレートするため、縦長の真っ黒な背景を用意してコマ落とし動画の各フレームを1ピクセルずつ下にずらしながら比較明合成を行った。従って時間の経過は上が古く、下が新しい。時間の経過が判りやすいよう、60コマ、すなわち1分毎に輝度を強調している。

ぱっと見、火星の直径の4～5倍振れていることがわかる。ウォーム1回転に1回、30秒くらいの間に火星3個分位動いている場所があり、これは見るからに何かがおかしいと思える症状だ。
もう一つ気になるのは全体としてナナメになっていること。これはウォーム1回転によりもっと長期的な追尾誤差の存在を示唆している。火星が画面で西に逃げているので追尾速度が僅かに遅過ぎる。先ず思いつくのがモータードライブの回転数誤差だが、これを制御している水晶発振回路の出力周波数をカウンターで測ると恒星時で追尾するのに必要な周波数に対してプラス0.01%弱の誤差がある。ナナメ加減から予想される数値に対して絶対値で1桁程小さいし、そもそもこれは早過ぎる方向の誤差なので、関係無さそうだ。これ以外の原因もいくつか思い浮かぶが、確かめる方法が無いので原因追究は保留しよう。

PEの量を測る角度の基準に火星の直径を使うつもりだったが、この画像から正確な数値を求めるのはちょっと無理だ。コンポジットした画像を作るべきだろう。伝統に則るとPEも画像の画素数を目で追うことになるが、画像は既にデータなので機械に画素を数えさせてしまえ、と言うことで合成前各フレームの火星の像の横方向の直径と中心位置を測ってプロットするPythonプログラムを書くことにした。

画素を数えるプログラムの概念

FFTを利用することも出来そうだが結果の解釈をしくじると悲惨な結果になるので、もっと素朴な直感的に間違いの無さそうな方法にした。各フレームの画像毎に以下の様に処理する。

1. 画像を「RGBモード」から「Lモード（白黒）」に変換する。これにより画像は[縦の画素数, 横の画素数]の大きさの2次元配列データとして扱える。
2. 画像を「火星の像の部分」と「背景」に2値化する。配列全体の最大値（最も明るい画素の値）の10%よりも大きな画素を「像の部分」とし、それ以外を「背景」とする。ここの閾値が、この処理で唯一任意に調整できるパラメータである。
3. 火星の像がかかっている配列の各行について「像の部分」に該当する列番号の最小値と最大値を求め、横幅＝最大値―最小値＋1、および中心位置＝（最大値＋最小値）÷2を計算する。
4. 横幅の最大値をそのフレームの像の直径とし、最大値が得られた行の中心位置をそのフレームの像の位置とする。

numpyを利用すると各段落の処理はほとんどが1～2行で書ける。ただし段落3の処理は配列の行数分のループ制御および条件分岐があり、段落4の処理と合わせて9行必要だ。また配列行数分のループ制御などをPythonで書いているので全体の処理は遅いが、20分間のデータ処理が3分程度で終わるので許容範囲だと思う。

各フレームについて計算した値から、動画全体について直径の平均値と位置の振れ幅（最大値―最小値）を計算している。

長期的な追尾誤差がありそうだと判っているので、中心位置の時間変化について直線回帰式の係数を求め、それから推定した各時点の長期的追尾誤差の量と、長期的追尾誤差を差し引いた正味PEの量を計算した。正味PE量についても振れ幅（最大値―最小値）を計算している。

結果のプロット（可視化）

計算の結果を時間に対してプロットした。主要な数値などは凡例内に示している。

位置のプロット（赤の実線）は、基本的に比較明合成画像を横倒しにしたものだが、目盛りが入って数値で語れることの違いは大きい。300秒付近で追尾速度がスーッと速くなり、像が数十画素左にずれている。500秒付近で、今度は何かにガツンとぶつかったような唐突な感じで追尾速度が遅くなり、像が数十画素右にずれたことが読み取れる。また同じことが900秒付近および1100秒付近で起こっていることもわかる。

長期的追尾誤差（赤の破線）は20分間で火星の直径1つ分ずれる程度だ。これを差し引いた後の位置の変化（青の実線）は、平均が300画素の位置になるようプロットしているが、概ねウォームの回転に伴う周期的変動になっているように見える。ただし、ウォームの回転より短い周期で変動する成分も相当量あるようだ。

プロットがランダムに数画素程度変動していることが判る。これは大気の乱れによる像の歪みが現れているもので、シーイングが悪いほど変動幅が大きくなる。また最後の30秒間の値が上昇しているのは、火星が雲に隠され像が不鮮明になった影響ではないかと思われる。

凡例に示した数値などは、大気の乱れや雲に隠されたことの影響を考慮していない。これは明らかに手抜きだが、考慮しようとすると結構煩雑だし恣意的な操作も出来てしまうので、このままにしておく。pk-pk値（最大値―最小値）は大気の乱れの分大きく出るが、平均値はほとんど影響を受けていないはずだ。この事を念頭にご覧いただきたい。この夜の火星の見かけの直径10.2秒角が61.2画素に写っているので、

周期的追尾誤差＝±(10.2×195.6÷61.2÷2)＝±16.3秒角
長期的追尾誤差＝10.2×0.0458×1200÷61.2＝9.16秒角/20分
総合的追尾誤差＝±(10.2×233.0÷61.2÷2)＝±19.4秒角/20分

見た目で覚悟していたよりも小さめの数値でホッとしたが、決して誇れたものではない。

PEの原因

初めてPEに気付いたとき、赤道儀のオーバーホールが必要なのだろうと思った。何と言っても30年物だからね。しかし、追尾誤差のプロットが「カックン、カックン」と動いているのを見て思い当たるフシがあった。私の使っているビクセンMD-5モータードライブは、軸もハウジングもリジッドに赤道儀と結合されていたのだ。

回転機屋さんが見たら腰を抜かしそうだが、購入時の取扱説明書に指示されていたか、あるいは取説が無くても付属の部品ではこれ以外の結合方法は考えられなかったかで、当時から私も相当違和感があったが一応動いているので「まぁいいか」で30年経過してしまった。とは言っても、通算してもウォームが何千回転もしていないと思うから、大きなダメージは無いだろう（これは気休め）。

これだと大なり小なり軸の中心がずれた状態なので、回転が渋くなる。偏芯して軸が回転するので、角度によって辛い方向とラクな方向が出てくる。ラクなエリアに入るところで加速して、辛いエリアに入るとき減速していると考えると辻褄が合う。300:1のギアヘッド付きPMステッピングモーターなのでそこそこ駆動トルクがあり、軸の回転数が0.1rpmと超低速なので今まで明るみに出ていなかっただけで、不都合が見逃されていたのだ。

ただし対策を講じて効果を確認する事しか、これが原因と確定する手段は無い。

セオリーはハウジングをリジッドに結合して回転軸はフレキシブルカップリングなどで結合する方法だが、この程度の軽いモーターなら軸をリジッドに結合してハウジングをフローティングにすることも選択できるだろうと考えた。これだとモーターを軸と軸受けで支えることになるが、モーターは軽いし、少なくとも両方ともリジッドに結合していたときよりストレスは少ないだろう、と言う理由。フローティングと言っても完全に自由にしたのでは回転が伝わらないので、回転方向のみ拘束し他の方向には自由に動けるようにする。

モーター側・赤道儀側、両方のステーに幅が6mmよりも少し広い長穴が明いているが、赤道儀側ステーの長穴は何故かテーパーになっていて、赤道儀と反対側では6mmよりも狭くなっている。そのため、赤道儀側ステーの穴にプリント基板をシャーシから浮かして取り付けるときに使う外径6φの金属製スペーサーがちょうど止められ、モーター側のステーに対しては適度なアソビが出る状態になる。写真はナナメから見ているので、アソビが大きく見えるが、実際は回転の接線方向で0.2mm以下だと思われる。

さあ、これでどうだ、と思って試運転すると、意外な伏兵が現れた。モーターの振動で金属同士が軽く触れ合っている箇所がカタカタと耳障りな音をたてるのだ。MD-5は運転中ほとんど音をたてなかった分、余計に耳障りである。対策としてショックコードと言うと聞こえがいいが、要はゴム紐で、回転方向に予圧を与えた。

これで耳障りな音は止む。回転方向以外の抵抗は、ほとんど増えていないはずだ。同じように火星を撮影した。

コマ落とし動画

比較明合成；

機械読み取りのプロット；

この夜の火星の見かけの直径は9.6秒角である。上の図の直径61.4画素が前回の61.2画素より大きくなっているのが気になるが・・・（前回よりシーイングが悪いせい？）；

周期的追尾誤差＝±(9.6×103.7÷61.4÷2)＝±8.1秒角
長期的追尾誤差＝9.6×0.08733×1200÷61.4＝16.4秒角/20分
総合的追尾誤差＝±(9.6×162.5÷61.4÷2)＝±12.7秒角/20分

周期的追尾誤差、すなわちPEが半減したので、対策は効果があったと言える。その一方で長期的追尾誤差が倍近く増加した。今回の対策で増える要素は無いと思うが、元々原因不明なのでなんとも言えない。

短い周期の追尾誤差

対策後の機械読み取りプロットの600秒付近に、10分より短い周期のサイン波のような波形が見える。その振幅は周期的追尾誤差全体の半分を超えているようであり、周期が10分より短い成分が対策後のPEに占める割合がかなり多い様である。この事を確かめるため、対策前後の、長期的追尾誤差を差し引いた後の位置の変動パターンについてフーリエ変換を行ってみた。なおフーリエ周波数ゼロの成分が極端に少ないのは、長期成分を差し引いているためである。またプロットの横軸フーリエ周波数は10分あたりのサイクル数、縦軸は各フーリエ周波数成分のサイン波の片側波高を画像上の画素数で表したものである。

赤の実線が対策前、緑が対策後である。

楽器の音や人間の声に倍音が含まれるのと同様、完全なサイン波ではない周期的波形には高調波が含まれる。対策前の波形はかなり歪んでいるので、多くの高調波が含まれている。その一方対策後では、周波数1の成分が1/3近くまで減っているのに加え、周波数4（周期2分半）の成分が多いのが目に付く。これは周波数1の成分の高調波と考えるには多過ぎる。素朴にMD-5の最終減速比が4:1で、出力軸の4倍の速さで回っている軸についている歯車が偏芯していると考えると辻褄が合う。また周波数10の成分も突出している様なので、次の減速比は2.5:1だろうか。

ギヤトレインの設計には素人なので妥当な推定かどうかわからないが、周波数4の成分がウォーム周り以外で発生していると考えるのは安全な範囲だろう。そう言う目で見ると、対策前の追尾誤差の波形にも周期150秒=2分半のサイン波が現れている、と言うか、コチラの方が明確である。この周波数4のパターンを目印にすると、不適切な軸の結合方法で発生していた追尾誤差成分をうまく分離できそうだ。

なお周波数4の成分は対策後の方が大きくなっているが、対策前では周波数1の成分の高調波とウォーム周り以外で発生した波形の位相が合わず、一部相殺されていたためだと考えられる。

対策により取り除かれた誤差成分の分離

対策前後で追尾状況記録開始時のウォーム回転角を特に合わせるような配慮はしていないので、このまま単純に引き算する訳にはいかない。タイミングを合わせる必要があるが、特に手がかりが無ければ600通りから一つ選ぶ必要がある。しかし、2分半周期の波形が対策前後で同じ位置に出ていると考えると、4通りに選択肢が絞れる。

先ず周期2分半の単位で山の位置を見ると、対策前605秒付近だったのが対策後は565秒付近が最寄りの山である。対策前後で40秒早くなっているが、これは2分半=150秒が単位となっているので、40＋150N（N=0,1,2,3）秒の可能性がある。

対策後の400～600秒の波形に注目すると、対策前の600～800秒とよく似ていることに気付く。これはN=1、すなわち190秒ずれている場合に該当しそうだ。

シーイングの影響を取り除くため、長期的追尾誤差を取り除いた対策後の波形から1分よりも短い周期の成分を取り除き190秒遅らせた波形と、長期的誤差とシーイングの影響を除去した対策前の波形を比較すると、対策により取り除かれた誤差成分が得られる。

一番上から、対策後の生の波形（緑色）とシーイングの影響を取り除いたもの（ピンク色）。生の波形の190秒付近にずらしたことによる継ぎ目が見えるが、赤色の線では取り除かれている。次が対策により取り除かれた成分の波形（黒色）。その下が対策前の生の波形（青色）とそれからシーイングの影響を取り除いたもの（赤色）である。この二つのプロットはほとんど重なっていて見分けにくいので、その差を10倍して一番下に黄色でプロットした。500秒付近と1100秒付近で大きな差が出ているのは、1分よりも短い期間に大きく動いていることが反映されている。

黒色の波形はほぼ矩形波で、想定したシナリオと辻褄が合う。ただし「それらしく見える」だけでN=0, 2, および3の場合の可能性を完全に排除するのは無理だが、特にこれ以上議論する材料も無いので、今回はこれを仮の結論としておく。

ここでシーイングの影響を取り除いた周期的追尾誤差曲線が得られたので、それに基づくと対策前の振れ幅が182.1画素で±15.2秒角、対策後が同90.2画素で±7.1秒角になる。

今回のまとめ

動画の各フレームを機械に読ませることで、比較的正確な数値が得られたと思う。ただし、ここで「正確」と言っているのは画像から数値を得る部分であって、実際に起こっていること全体に関しては見方を変える必要があるだろう。今回特に気になるのが、火星の像の大きさと見かけの大きさの関係だ。恒星ペアの角距離に比べてシーイングなどの影響を受けやすいようである。これを考慮すると、結果には10%程度の誤差があってもおかしくなく、今回一律に小数点1桁まで計算しているのは不適切である。しかしこれは、扱っている桁数があまりに少ないと手を抜いているように見えるかもしれない、と言う見栄もあるので、勘弁して欲しい。

永年見過ごしてきた問題に対策が打てたのは大きな収穫だ。今までも視野の中で星が動いたのに何回か気付いたことがあるが、短時間の内に大きく動いたので「三脚の脚がすべったかな」程度に考え、モーターの取り付けを疑うことは無かった。やはりセオリーを無視してはいけない。

長期的な追尾誤差があるのには困ったが、カメラがデジタルになって銀塩のころに比べるとシャッターを開けている時間が格段に短くなっているので、実用上問題にはならないだろう。今まで星野写真を撮った事が無いが、これならちょっとやってみようか、と言う気になるが、天文道楽が復活してつくづく思うのは「スカッと晴れた夜は稀」と言うこと。天気とシーイングには勝てない。