(57)【特許請求の範囲】
【請求項１】
移動体の運動を解析するための運動解析装置であって、
前記移動体に取り付けられ、角速度センサ、加速度センサ及び地磁気センサの少なくとも１つを含むセンサユニットから出力される時系列の角速度データ、加速度データ及び地磁気データの少なくとも１つを含む時系列のセンサデータを取得する第１取得部と、
第１距離画像センサにより前記移動体の運動を撮影した時系列の深度画像及び二次元画像を含む時系列の画像データを取得する第２取得部と、
前記センサデータ及び前記画像データに基づいて、前記移動体の三次元挙動を導出する導出部と
を備え、
前記導出部は、前記画像データ及び前記センサデータを用いて定義される関数を含む目的関数を最小化又は最大化する最適解として、前記移動体の姿勢を導出する、
運動解析装置。
【請求項２】
前記関数は、前記画像データから導出される前記移動体の向きを用いて定義される、
請求項１に記載の運動解析装置。
【請求項３】
前記関数は、前記地磁気データ及び加速度データの少なくとも一方を用いて定義される、
請求項１又は２に記載の運動解析装置。
【請求項４】
前記導出部は、前記移動体の初期の姿勢及び前記時系列の角速度データを用いて、前記移動体の仮の姿勢を算出し、
前記関数は、前記仮の姿勢を用いて定義される、
請求項１から３のいずれかに記載の運動解析装置。
【請求項５】
前記導出部は、前記画像データから導出される初期の前記移動体の向き及び初期の前記加速度データ、又は、前記角速度データを用いて、前記移動体の初期の姿勢として、前記移動体の静止時の姿勢を導出する、
請求項４に記載の運動解析装置。
【請求項６】
移動体の運動を解析するための運動解析装置であって、
前記移動体に取り付けられ、角速度センサ、加速度センサ及び地磁気センサの少なくとも１つを含むセンサユニットから出力される時系列の角速度データ、加速度データ及び地磁気データの少なくとも１つを含む時系列のセンサデータを取得する第１取得部と、
第１距離画像センサにより前記移動体の運動を撮影した時系列の深度画像及び二次元画像を含む時系列の画像データを取得する第２取得部と、
前記センサデータ及び前記画像データに基づいて、前記移動体の三次元挙動を導出する導出部と
を備え、
前記導出部は、前記移動体の姿勢及び前記移動体の姿勢を微分フィルタによりフィルタリングした微分データに基づいて、前記移動体の角速度を導出する、
運動解析装置。
【請求項７】
移動体の運動を解析するための運動解析装置であって、
前記移動体に取り付けられ、角速度センサ、加速度センサ及び地磁気センサの少なくとも１つを含むセンサユニットから出力される時系列の角速度データ、加速度データ及び地磁気データの少なくとも１つを含む時系列のセンサデータを取得する第１取得部と、
第１距離画像センサにより前記移動体の運動を撮影した時系列の深度画像及び二次元画像を含む時系列の画像データを取得する第２取得部と、
前記センサデータ及び前記画像データに基づいて、前記移動体の三次元挙動を導出する導出部と
を備え、
前記導出部は、前記加速度データを用いて、前記画像データから導出される前記移動体の仮の位置及び当該仮の位置を微分した仮の速度の少なくとも一方をスムージングすることにより、前記移動体の位置及び速度の少なくとも一方を導出する、
運動解析装置。
【請求項８】
前記導出部は、前記移動体の速度を微分フィルタによりフィルタリングすることにより、前記移動体の加速度を導出する、
請求項７に記載の運動解析装置。
【請求項９】
移動体の運動を解析するための運動解析装置であって、
前記移動体に取り付けられ、角速度センサ、加速度センサ及び地磁気センサの少なくとも１つを含むセンサユニットから出力される時系列の角速度データ、加速度データ及び地磁気データの少なくとも１つを含む時系列のセンサデータを取得する第１取得部と、
第１距離画像センサにより前記移動体の運動を撮影した時系列の深度画像及び二次元画像を含む時系列の画像データを取得する第２取得部と、
前記センサデータ及び前記画像データに基づいて、前記移動体の三次元挙動を導出する導出部と
を備え、
前記導出部は、前記画像データから導出される前記移動体の向きの時系列データと、前記センサデータから導出される前記移動体の向きの時系列データとの一致度が最も高くなるように、前記時系列の画像データと前記時系列のセンサデータとの時刻合わせを行う、
運動解析装置。
【請求項１０】
移動体の運動を解析するための運動解析装置であって、
前記移動体に取り付けられ、角速度センサ、加速度センサ及び地磁気センサの少なくとも１つを含むセンサユニットから出力される時系列の角速度データ、加速度データ及び地磁気データの少なくとも１つを含む時系列のセンサデータを取得する第１取得部と、
第１距離画像センサにより前記移動体の運動を撮影した時系列の深度画像及び二次元画像を含む時系列の画像データを取得する第２取得部と、
前記センサデータ及び前記画像データに基づいて、前記移動体の三次元挙動を導出する導出部と
を備え、
前記導出部は、前記画像データに基づいて第１時刻における前記移動体の姿勢である瞬時姿勢を導出し、前記移動体の姿勢が前記第１時刻において前記瞬時姿勢に一致するように、前記第１時刻と異なる第２時刻と前記第１時刻との間の前記移動体の姿勢を補間する、
運動解析装置。
【請求項１１】
第２距離画像センサにより前記第１距離画像センサとは異なる方向から前記移動体の運動を撮影した時系列の深度画像及び二次元画像を含む時系列の画像データを取得する第３取得部
をさらに備え、
前記導出部は、前記センサデータと、前記第２取得部及び前記第３取得部により取得された前記画像データとに基づいて、前記移動体の三次元挙動を導出する、
請求項１から１０のいずれかに記載の運動解析装置。
【請求項１２】
移動体の運動を解析するための運動解析方法であって、
前記移動体に取り付けられ、角速度センサ、加速度センサ及び地磁気センサの少なくとも１つを含むセンサユニットを用いて、時系列の角速度データ、加速度データ及び地磁気データの少なくとも１つを含む時系列のセンサデータを取得するステップと、
距離画像センサを用いて、前記移動体の運動を撮影した時系列の深度画像及び二次元画像を含む時系列の画像データを取得するステップと、
前記センサデータ及び前記画像データに基づいて、前記移動体の三次元挙動を導出するステップと、
を備え、
前記導出するステップでは、前記画像データから導出される前記移動体の向きの時系列データと、前記センサデータから導出される前記移動体の向きの時系列データとの一致度が最も高くなるように、前記時系列の画像データと前記時系列のセンサデータとの時刻合わせを行う、
運動解析方法。
【請求項１３】
移動体の運動を解析するための運動解析プログラムであって、
前記移動体に取り付けられ、角速度センサ、加速度センサ及び地磁気センサの少なくとも１つを含むセンサユニットから出力される時系列の角速度データ、加速度データ及び地磁気データの少なくとも１つを含む時系列のセンサデータを取得するステップと、
距離画像センサにより前記移動体の運動を撮影した時系列の深度画像及び二次元画像を含む時系列の画像データを取得するステップと、
前記センサデータ及び前記画像データに基づいて、前記移動体の三次元挙動を導出するステップと、
をコンピュータに実行させ、
前記導出するステップでは、前記画像データから導出される前記移動体の向きの時系列データと、前記センサデータから導出される前記移動体の向きの時系列データとの一致度が最も高くなるように、前記時系列の画像データと前記時系列のセンサデータとの時刻合わせを行う、
運動解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、移動体の運動を解析するための運動解析装置、方法及びプログラムに関する。
【背景技術】
【０００２】
従来より、角速度センサ、加速度センサ及び地磁気センサを含む慣性センサを移動体に取り付け、当該慣性センサの出力値に基づいて移動体の挙動を推定する技術が知られている（特許文献１等）。また、カメラにより移動体の運動を撮影し、得られた画像データに基づいて移動体の挙動を推定することも行われている（特許文献２，３等）。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２０１１−２２７０１７号公報
【特許文献２】特開２０１３−２１５３４９号公報
【特許文献３】特開２００４−３１３４７９号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかしながら、特許文献１のように慣性センサの出力値に基づいて解析を行う場合、当該出力値には慣性センサのドリフト成分等による誤差が含まれるため、解析の精度が低下することがある。例えば、慣性センサの出力値である加速度のデータを積分して速度や位置を求めるような場合、この積分により誤差が蓄積し得る。また、移動体が高速で運動する場合には、特に角速度センサの誤差が大きくなり、移動体の姿勢の推定が困難となり得る。また、回転運動が大きい場合には特に、加速度センサの出力値に含まれる重力加速度の影響を正しく評価することが困難となり得る。
【０００５】
一方、特許文献２，３のように画像データに基づいて解析を行う場合にも、しばしば精度の面で問題が生じる。例えば、画像処理により移動体の位置を推定し、当該位置のデータを微分することで速度や加速度を求めるような場合、位置の推定における誤差が僅かであっても、微分により誤差が増大し得る。
【０００６】
本発明は、高精度に移動体の運動を解析することができる運動解析装置、方法及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００７】
本発明の第１観点に係る運動解析装置は、移動体の運動を解析するための運動解析装置であって、第１取得部と、第２取得部と、導出部とを備える。前記第１取得部は、センサユニットから出力される時系列の角速度データ、加速度データ及び地磁気データの少なくとも１つを含む時系列のセンサデータを取得する。前記センサユニットは、前記移動体に取り付けられ、角速度センサ、加速度センサ及び地磁気センサの少なくとも１つを含む。前記第２取得部は、第1距離画像センサにより前記移動体の運動を撮影した時系列の深度画像及び二次元画像を含む時系列の画像データを取得する。前記導出部は、前記センサデータ及び前記画像データに基づいて、前記移動体の三次元挙動を導出する。
【０００８】
本発明の第２観点に係る運動解析装置は、第１観点に係る運動解析装置であって、前記導出部は、前記画像データ及び前記センサデータを用いて定義される関数を含む目的関数を最小化又は最大化する最適解として、前記移動体の姿勢を導出する。
【０００９】
本発明の第３観点に係る運動解析装置は、第２観点に係る運動解析装置であって、前記関数は、前記画像データから導出される前記移動体の向きを用いて定義される。
【００１０】
本発明の第４観点に係る運動解析装置は、第２又は第３観点に係る運動解析装置であって、前記関数は、前記地磁気データ及び加速度データの少なくとも一方を用いて定義される。
【００１１】
本発明の第５観点に係る運動解析装置は、第２から第４観点のいずれかに係る運動解析装置であって、前記導出部は、前記移動体の初期の姿勢及び前記時系列の角速度データを用いて、前記移動体の仮の姿勢を算出する。前記関数は、前記仮の姿勢を用いて定義される。
【００１２】
本発明の第６観点に係る運動解析装置は、第５観点に係る運動解析装置であって、前記導出部は、前記画像データから導出される初期の前記移動体の向き及び初期の前記加速度データ、又は、前記角速度データを用いて、前記移動体の初期の姿勢として、前記移動体の静止時の姿勢を導出する。
【００１３】
本発明の第７観点に係る運動解析装置は、第１から第６観点のいずれかに係る運動解析装置であって、前記導出部は、前記移動体の姿勢及び前記移動体の姿勢を微分フィルタによりフィルタリングした微分データに基づいて、前記移動体の角速度を導出する。
【００１４】
本発明の第８観点に係る運動解析装置は、第１から第７観点のいずれかに係る運動解析装置であって、前記導出部は、前記加速度データ用いて、前記画像データから導出される前記移動体の仮の位置及び当該仮の位置を微分した仮の速度の少なくとも一方をスムージングすることにより、前記移動体の位置及び速度の少なくとも一方を導出する。
【００１５】
本発明の第９観点に係る運動解析装置は、第８観点に係る運動解析装置であって、前記導出部は、前記移動体の速度を微分フィルタによりフィルタリングすることにより、前記移動体の加速度を導出する。
【００１６】
本発明の第１０観点に係る運動解析装置は、第１から第９観点のいずれかに係る運動解析装置であって、前記導出部は、前記画像データから導出される前記移動体の向きの時系列データと、前記センサデータから導出される前記移動体の向きの時系列データとの一致度が最も高くなるように、前記時系列の画像データと前記時系列のセンサデータとの時刻合わせを行う。
【００１７】
本発明の第１１観点に係る運動解析装置は、第１観点から第１０観点のいずれかに係る運動解析装置であって、第３取得部をさらに備える。前記第３取得部は、第２距離画像センサにより前記第１距離画像センサとは異なる方向から前記移動体の運動を撮影した時系列の深度画像及び二次元画像を含む時系列の画像データを取得する。前記導出部は、前記センサデータと、前記第２取得部及び前記第３取得部により取得された前記画像データとに基づいて、前記移動体の三次元挙動を導出する。
【００１８】
本発明の第１２観点に係る運動解析装置は、第１観点から第１１観点のいずれかに係る運動解析装置であって、前記導出部は、前記画像データに基づいて第１時刻における前記移動体の姿勢である瞬時姿勢を導出し、前記移動体の姿勢が前記第１時刻において前記瞬時姿勢に一致するように、前記第１時刻と異なる第２時刻と前記第１時刻との間の前記移動体の姿勢を補間する。
【００１９】
本発明の第１３観点に係る運動解析方法は、移動体の運動を解析するための運動解析方法であって、以下の（１）〜（３）のステップを含む。
（１）センサユニットを用いて、時系列の角速度データ、加速度データ及び地磁気データの少なくとも１つを含む時系列のセンサデータを取得するステップ。前記センサユニットは、前記移動体に取り付けられ、角速度センサ、加速度センサ及び地磁気センサの少なくとも１つを含む。
（２）距離画像センサを用いて、前記移動体の運動を撮影した時系列の深度画像及び二次元画像を含む時系列の画像データを取得するステップ。
（３）前記センサデータ及び前記画像データに基づいて、前記移動体の三次元挙動を導出するステップ。
【００２０】
本発明の第１４観点に係る運動解析プログラムは、移動体の運動を解析するための運動解析プログラムであって、以下の（１）〜（３）のステップをコンピュータに実行させる。
（１）センサユニットから出力される時系列の角速度データ、加速度データ及び地磁気データの少なくとも１つを含む時系列のセンサデータを取得するステップ。前記センサユニットは、前記移動体に取り付けられ、角速度センサ、加速度センサ及び地磁気センサの少なくとも１つを含む。
（２）距離画像センサにより前記移動体の運動を撮影した時系列の深度画像及び二次元画像を含む時系列の画像データを取得するステップ。
（３）前記センサデータ及び前記画像データに基づいて、前記移動体の三次元挙動を導出するステップ。
【発明の効果】
【００２１】
本発明によれば、移動体の運動が、（１）時系列の角速度データ、加速度データ及び地磁気データの少なくとも１つを含む時系列のセンサデータ、及び、（２）時系列の深度画像及び二次元画像を含む時系列の画像データに基づいて導出される。その結果、センサデータに含まれる誤差は、画像データに基づいて軽減される。或いは、画像データに含まれる誤差は、センサデータに基づいて軽減される。よって、高精度に移動体の運動を解析することができる。
【図面の簡単な説明】
【００２２】
【図１】本発明の第１実施形態に係る運動解析装置を含む運動解析システムの全体構成を示す図。
【図２】図１の運動解析システムの機能ブロック図。
【図３】ゴルフクラブのグリップ近傍を基準とするｘｙｚ局所座標系を説明する図。
【図４】（Ａ）アドレス状態を示す図。（Ｂ）トップ状態を示す図。（Ｃ）インパクト状態を示す図。（Ｄ）フィニッシュ状態を示す図。
【図５】第１実施形態に係る運動解析方法の流れを示すフローチャート。
【図６】第１実施形態に係る時刻合わせ処理の流れを示すフローチャート。
【図７】二次元のシャフトの向きを説明する図。
【図８】時刻合わせ処理の原理を説明する概念図。
【図９】アドレス期間の姿勢の導出処理の流れを示すフローチャート。
【図１０】アドレスからインパクト付近までの姿勢の導出処理の流れを示すフローチャート。
【図１１】グリップの位置及び速度の導出処理の流れを示すフローチャート。
【図１２】グリップの角速度及び角加速度の導出処理の流れを示すフローチャート。
【図１３】グリップの加速度の導出処理の流れを示すフローチャート。
【図１４】実施例１及び参考例１を示すグラフ。
【図１５】実施例２及び参考例２を示すグラフ。
【図１６Ａ】実施例３及び参考例３を示すグラフ（グリップのＸ軸方向の位置）。
【図１６Ｂ】実施例３及び参考例３を示すグラフ（グリップのＹ軸方向の位置）。
【図１６Ｃ】実施例３及び参考例３を示すグラフ（グリップのＺ軸方向の位置）。
【図１７Ａ】実施例４及び参考例４を示すグラフ（グリップのＸ軸方向の位置）。
【図１７Ｂ】実施例４及び参考例４を示すグラフ（グリップのＹ軸方向の位置）。
【図１７Ｃ】実施例４及び参考例４を示すグラフ（グリップのＺ軸方向の位置）。
【図１８】本発明の第２実施形態に係る運動解析装置を含む運動解析システムの全体構成を示す図。
【図１９】図１８の運動解析システムの機能ブロック図。
【図２０】第２実施形態に係る運動解析方法の流れを示すフローチャート。
【図２１】第２実施形態に係る時刻合わせ処理の流れを示すフローチャート。
【図２２】ダウンスイング期間の姿勢の補正処理の流れを示すフローチャート。
【図２３】実施例５を示すグラフ。
【図２４】変形例に係る初期時刻の姿勢の導出処理の流れを示すフローチャート。
【発明を実施するための形態】
【００２３】
以下、図面を参照しつつ、本発明に係る運動解析装置、方法及びプログラムをゴルフスイングの解析に利用した場合の幾つかの実施形態について説明する。
【００２４】
＜１．第１実施形態＞
＜１−１．運動解析システムの概要＞
図１及び図２に、本発明の第１実施形態に係る運動解析装置１を含む運動解析システム１００の全体構成図を示す。運動解析システム１００は、ゴルファー７によるゴルフクラブ５のスイング動作の解析に適するように構成されている。運動解析装置１は、ゴルフクラブ５に取り付けられた慣性センサユニット４から出力される時系列の角速度データ、加速度データ及び地磁気データ（以下、センサデータということがある）に加えて、スイング動作を撮影した時系列の画像データ（すなわち、動画データ）に基づいて、スイング動作を解析する。以上の撮影は、距離画像センサ２により行われる。運動解析装置１は、慣性センサユニット４及び距離画像センサ２とともに運動解析システム１００を構成する。運動解析装置１による解析の結果は、ゴルファー７に適したゴルフクラブ５のフィッティングや、ゴルファー７のフォームの改善、ゴルフ用品の開発等、様々な用途で利用される。
【００２５】
本実施形態では、主として、スイング動作中に移動するゴルフクラブ５のグリップ５１及びシャフト５２の三次元挙動が解析される。より具体的には、スイング動作中のグリップ５１の位置（三次元座標）、速度、加速度、角速度及び角加速度とともに、シャフト５２の姿勢が導出される。
【００２６】
＜１−２．各部の詳細＞
以下、運動解析システム１００の各部の詳細について説明する。
【００２７】
＜１−２−１．距離画像センサ＞
距離画像センサ２は、ゴルファー７がゴルフクラブ５を試打する様子を二次元画像として撮影するとともに、被写体までの距離を測定する測距機能を有するカメラである。従って、距離画像センサ２は、時系列の二次元画像とともに、時系列の深度画像を出力することができる。なお、ここでいう二次元画像とは、撮影空間の像をカメラの光軸に直交する平面内へ投影した画像である。また、深度画像とは、カメラの光軸方向の被写体の奥行きのデータを、二次元画像と略同じ撮像範囲内の画素に割り当てた画像である。
【００２８】
本実施形態で使用される距離画像センサ２は、二次元画像を赤外線画像（以下、ＩＲ画像という）として撮影する。また、深度画像は、赤外線を用いたタイムオブフライト方式やドットパターン投影方式等の方法により得られる。従って、図１に示すように、距離画像センサ２は、赤外線を前方に向けて発光するＩＲ発光部２１と、ＩＲ発光部２１から照射され、被写体に反射して戻ってきた赤外線を受光するＩＲ受光部２２とを有する。ＩＲ受光部２２は、光学系及び撮像素子等を有するカメラである。ドットパターン投影方式では、ＩＲ発光部２１から照射された赤外線のドットパターンをＩＲ受光部２２で読み取り、距離画像センサ２内部での画像処理によりドットパターンを検出し、これに基づいて奥行きが計算される。本実施形態では、ＩＲ発光部２１及びＩＲ受光部２２は、同じ筐体２０内に収容され、筐体２０の前方に配置されている。また、本実施形態に係る距離画像センサ２では、奥行きを計測可能な有効距離が定められており、深度画像には、奥行きの値が得られない画素が含まれ得る。
【００２９】
距離画像センサ２には、距離画像センサ２の動作全体を制御するＣＰＵ２３の他、撮影された画像データを少なくとも一時的に記憶するメモリ２４が内蔵されている。距離画像センサ２の動作を制御する制御プログラムは、メモリ２４内に格納されている。また、距離画像センサ２には、通信部２５も内蔵されており、当該通信部２５は、撮影された画像データを、有線又は無線の通信線１７を介して、運動解析装置１等の外部のデバイスへと出力することができる。本実施形態では、ＣＰＵ２３及びメモリ２４も、ＩＲ発光部２１及びＩＲ受光部２２とともに、筐体２０内に収納されている。なお、運動解析装置１への画像データの受け渡しは、必ずしも通信部２５を介して行う必要はない。例えば、メモリ２４が着脱式であれば、これを筐体２０内から取り外し、運動解析装置１のリーダー（後述する通信部１５に対応）に挿入する等して、運動解析装置１で画像データを読み出すことができる。
【００３０】
本実施形態では、以上のとおり、距離画像センサ２により赤外線撮影が行われ、撮影されたＩＲ画像に基づいて、グリップ５１の軌道が追跡される。従って、この追跡が容易となるように、グリップ５１（より正確には、グリップ端５１ａ）には、赤外線を効率的に反射する反射シートがマーカーとして貼付される。また、シャフト５２にも、同様の赤外線の反射シートがマーカーとして貼付される。本実施形態では、距離画像センサ２は、ゴルファー７を正面側から撮影すべく、ゴルファー７の前方に設置される。
【００３１】
＜１−２−２．慣性センサユニット＞
慣性センサユニット４は、図１及び図３に示すとおり、ゴルフクラブ５のグリップ端５１ａに取り付けられており、グリップ端５１ａの挙動を計測する。本実施形態に係る慣性センサユニット４は、着脱自在に構成されており、任意のゴルフクラブ５に取り付けることができる。なお、ゴルフクラブ５は、一般的なゴルフクラブであり、シャフト５２と、シャフト５２の一端に設けられたヘッド５３と、シャフト５２の他端に設けられたグリップ５１とから構成される。慣性センサユニット４は、スイング動作の妨げとならないよう、小型且つ軽量に構成されている。図２に示すように、本実施形態に係る慣性センサユニット４には、加速度センサ４１、角速度センサ４２及び地磁気センサ４３が搭載されている。また、慣性センサユニット４には、これらのセンサ４１〜４３から出力されるセンサデータを運動解析装置１等の外部のデバイスに送信するための通信装置４０も搭載されている。なお、本実施形態では、通信装置４０は、スイング動作の妨げにならないように無線式であるが、ケーブルを介して有線式に運動解析装置１に接続するようにしてもよい。
【００３２】
加速度センサ４１、角速度センサ４２及び地磁気センサ４３はそれぞれ、ｘｙｚ局所座標系における加速度、角速度及び地磁気を計測する。より具体的には、加速度センサ４１は、ｘ軸、ｙ軸及びｚ軸方向の加速度ａx，ａy，ａzを計測する。角速度センサ４２は、ｘ軸、ｙ軸及びｚ軸周りの角速度ωx，ωy，ωzを計測する。地磁気センサ４３は、ｘ軸、ｙ軸及びｚ軸方向の地磁気ｍx，ｍy，ｍzを計測する。これらのセンサデータは、所定のサンプリング周期Δｔの時系列データとして取得される。なお、ｘｙｚ局所座標系は、図３に示すとおりに定義される３軸直交座標系である。すなわち、ｚ軸は、シャフト５２の延びる方向に一致し、ヘッド５３からグリップ５１に向かう方向が、ｚ軸正方向である。ｙ軸は、ゴルフクラブ５のアドレス時の飛球方向にできる限り沿うように、すなわち、フェース−バック方向に概ね沿うように配向される。ｘ軸は、ｙ軸及びｚ軸に直交するように、すなわち、トゥ−ヒール方向に概ね沿うように配向され、ヒール側からトゥ側に向かう方向がｘ軸正方向である。
【００３３】
本実施形態では、加速度センサ４１、角速度センサ４２及び地磁気センサ４３からのセンサデータは、通信装置４０を介してリアルタイムに運動解析装置１に送信される。しかしながら、例えば、慣性センサユニット４内の記憶装置にセンサデータを格納しておき、スイング動作の終了後に当該記憶装置からセンサデータを取り出して、運動解析装置１に受け渡すようにしてもよい。
【００３４】
＜１−２−３．運動解析装置＞
運動解析装置１は、ハードウェアとしては汎用のパーソナルコンピュータであり、例えば、タブレットコンピュータ、スマートフォン、ノート型コンピュータ、デスクトップ型コンピュータとして実現される。図２に示すとおり、運動解析装置１は、ＣＤ−ＲＯＭ等のコンピュータで読み取り可能な記録媒体３０から運動解析プログラム３を汎用のコンピュータにインストールすることにより製造される。運動解析プログラム３は、慣性センサユニット４から送られてくるセンサデータ及び距離画像センサ２から送られてくる画像データに基づいて、ゴルフスイングを解析するためのソフトウェアであり、運動解析装置１に後述する動作を実行させる。
【００３５】
運動解析装置１は、表示部１１、入力部１２、記憶部１３、制御部１４及び通信部１５を備える。これらの部１１〜１５は、互いにバス線１６を介して接続されており、相互に通信可能である。表示部１１は、液晶ディスプレイ等で構成することができ、ゴルフスイングの解析の結果等をユーザに対し表示する。なお、ここでいうユーザとは、ゴルファー７自身やそのインストラクター、ゴルフ用品の開発者等、ゴルフスイングの解析の結果を必要とする者の総称である。入力部１２は、マウス、キーボード、タッチパネル等で構成することができ、運動解析装置１に対するユーザからの操作を受け付ける。
【００３６】
記憶部１３は、ハードディスク等で構成することができる。記憶部１３内には、運動解析プログラム３が格納されている他、慣性センサユニット４から送られてくるセンサデータが保存される。また、記憶部１３内には、距離画像センサ２から送られてくる画像データも保存される。制御部１４は、ＣＰＵ、ＲＯＭおよびＲＡＭ等から構成することができる。制御部１４は、記憶部１３内の運動解析プログラム３を読み出して実行することにより、仮想的に第１取得部１４ａ、第２取得部１４ｂ、導出部１４ｃ及び表示制御部１４ｄとして動作する。各部１４ａ〜１４ｄの動作の詳細については、後述する。通信部１５は、慣性センサユニット４や距離画像センサ２等の外部のデバイスから通信線１７を介してデータを受信する通信インターフェースとして機能する。
【００３７】
＜１−３．運動解析方法＞
以下、図５を参照しつつ、ゴルフスイングを対象とする運動解析方法について説明する。この運動解析方法は、第１データ収集処理（Ｓ１）、第２データ収集処理（Ｓ２）及び運動解析処理（Ｓ３〜Ｓ９）により実現される。第１データ収集処理は、慣性センサユニット４から出力されるセンサデータを収集する処理であり、第２データ収集処理は、距離画像センサ２から出力される画像データを収集する処理である。運動解析処理は、当該センサデータ及び当該画像データに基づいて、運動解析装置１によりゴルフクラブ５の三次元挙動を解析する処理である。以下、これらの処理について、順に説明する。
【００３８】
＜１−３−１．第１データ収集処理＞
第１データ収集処理では、ゴルファー７により、上述の慣性センサユニット４付きゴルフクラブ５がスイングされる。このとき、慣性センサユニット４により、ゴルフスイング中の加速度ａx，ａy，ａz、角速度ωx，ωy，ωz及び地磁気ｍx，ｍy，ｍzのセンサデータが検出される。また、これらのセンサデータは、慣性センサユニット４の通信装置４０を介して運動解析装置１に送信される。一方、運動解析装置１側では、第１取得部１４ａが通信部１５を介してこれを受信し、記憶部１３内に格納する。本実施形態では、少なくともアドレスの少し前（以下、初期時刻ｔ1という）からフィニッシュまでの時系列のセンサデータが収集される。
【００３９】
なお、ゴルフクラブのスイング動作は、一般に、アドレス、トップ、インパクト、フィニッシュの順に進む。アドレスとは、図４（Ａ）に示すとおり、ゴルフクラブ５のヘッド５３をボール近くに配置した静止状態を意味し、トップとは、図４（Ｂ）に示すとおり、アドレスからゴルフクラブ５をテイクバックし、最もヘッド５３が振り上げられた状態を意味する。インパクトとは、図４（Ｃ）に示すとおり、トップからゴルフクラブ５が振り下ろされ、ヘッド５３がボールと衝突した瞬間の状態を意味し、フィニッシュとは、図４（Ｄ）に示すとおり、インパクト後、ゴルフクラブ５を前方へ振り抜いた状態を意味する。
【００４０】
＜１−３−２．第２データ収集処理＞
第２データ収集処理では、ゴルフクラブ５のスイング動作が距離画像センサ２により撮影される。すなわち、距離画像センサ２により、ゴルフスイングを捉えたＩＲ画像及び深度画像を含む画像データが検出される。第２データ収集処理は、第１データ収集処理と並行して、第１データ収集処理と同じスイング動作を対象として行われる。検出された画像データは、距離画像センサ２の通信部２５を介して運動解析装置１に送信される。一方、運動解析装置１側では、第２取得部１４ｂが通信部１５を介してこれを受信し、記憶部１３内に格納する。本実施形態では、少なくともアドレスの少し前の初期時刻ｔ1からフィニッシュまでの時系列の画像データが収集される。
【００４１】
＜１−３−３．運動解析処理＞
運動解析処理の大まかな流れは、以下のとおりである。すなわち、まず、第１データ収集処理により取得された時系列のセンサデータと、第２データ収集処理により取得された時系列の画像データとの時刻合わせが行われる（Ｓ３）。言い換えると、導出部１４ｃが、センサデータと画像データとを同期させる。
【００４２】
続くステップＳ４では、導出部１４ｃが、初期時刻ｔ1からアドレスまでの期間（以下、アドレス期間という）のシャフト５２の姿勢を導出する。このとき、本実施形態では、ラグランジュの未定乗数法を用いた最適化が実行される。なお、多くのゴルファー７は、アドレス期間にゴルフクラブ５を軽く前後にスイングさせるワッグル動作を行う。しかしながら、アドレス以降のスイング動作時に比べると、アドレス期間において、ゴルフクラブ５、特にグリップ５１は殆ど運動しない。このため、ステップＳ４では、静止している移動体及びほぼ静止している移動体の運動の解析に適したアルゴリズムが用いられる。
【００４３】
続くステップＳ５では、導出部１４ｃが、アドレスからインパクト付近までのシャフト５２の姿勢を導出する。従って、ステップＳ５では、静止している移動体及びほぼ静止している移動体の運動のみならず、高速で運動する移動体の解析にも適したアルゴリズムが用いられる。本実施形態では、ここでも、ラグランジュの未定乗数法を用いた最適化が実行される。
【００４４】
続くステップＳ６では、導出部１４ｃが、アドレスからインパクト付近までのグリップ５１の位置及び速度を導出する。ここでは、本実施形態では、カルマンフィルタを用いる等したスムージングが実行される。
【００４５】
続くステップＳ７では、導出部１４ｃが、アドレスからインパクト付近までのグリップ５１の角速度及び角加速度を導出する。ここでは、本実施形態では、疑似微分フィルタを用いたフィルタリングが実行される。
【００４６】
続くステップＳ８では、導出部１４ｃが、アドレスからインパクト付近までのグリップ５１の加速度を導出する。本実施形態では、ここでも、疑似微分フィルタを用いたフィルタリングが実行される。
【００４７】
続くステップＳ９では、表示制御部１４ｄが、ステップＳ４〜Ｓ８で導出された三次元挙動を示す情報を表示部１１上に表示する。このとき、例えば、数値情報、文字情報として表示することもできるし、グラフやコンピュータグラフィックス等を用いて図式的に表示することもできる。
【００４８】
＜１−３−３−１．時刻合わせ（Ｓ３）＞
以下、図６を参照しつつ、上記ステップＳ３における、時系列のセンサデータと時系列の画像データとの時刻合わせ処理について説明する。
【００４９】
まず、ステップＳ１１において、導出部１４ｃは、記憶部１３内のセンサデータを参照し、アドレス、トップ及びインパクトの時刻を導出する。これらの時刻の導出方法としては、様々なものが公知であるため、ここでは詳細な説明を省略する。
【００５０】
続くステップＳ１２では、導出部１４ｃは、記憶部１３内のセンサデータに基づいて、アドレスからインパクト付近までの各時刻における姿勢行列Ｎを導出する。今、姿勢行列Ｎを以下の式で表す。姿勢行列Ｎは、ｘｙｚ局所座標系の値をＸＹＺ慣性座標系の値に変換するための行列である。
【数１】
0006776882.tif 000002
【００５１】
ＸＹＺ慣性座標系は、図１に示すとおりに定義される３軸直交座標系である。すなわち、Ｚ軸は、鉛直下方から上方に向かう方向であり、Ｘ軸は、ゴルファー７の背から腹に向かう方向であり、Ｙ軸は、地平面に平行でボールの打球地点から目標地点に向かう方向である。
【００５２】
姿勢行列Ｎの９つの成分の意味は、以下のとおりである。
成分ａ：慣性座標系のＸ軸と、局所座標系のｘ軸とのなす角度の余弦
成分ｂ：慣性座標系のＹ軸と、局所座標系のｘ軸とのなす角度の余弦
成分ｃ：慣性座標系のＺ軸と、局所座標系のｘ軸とのなす角度の余弦
成分ｄ：慣性座標系のＸ軸と、局所座標系のｙ軸とのなす角度の余弦
成分ｅ：慣性座標系のＹ軸と、局所座標系のｙ軸とのなす角度の余弦
成分ｆ：慣性座標系のＺ軸と、局所座標系のｙ軸とのなす角度の余弦
成分ｇ：慣性座標系のＸ軸と、局所座標系のｚ軸とのなす角度の余弦
成分ｈ：慣性座標系のＹ軸と、局所座標系のｚ軸とのなす角度の余弦
成分ｉ：慣性座標系のＺ軸と、局所座標系のｚ軸とのなす角度の余弦
ここで、ベクトル（ａ，ｂ，ｃ）は、ｘ軸方向の単位ベクトルを表し、ベクトル（ｄ，ｅ，ｆ）は、ｙ軸方向の単位ベクトルを表し、ベクトル（ｇ，ｈ，ｉ）は、ｚ軸方向の単位ベクトルを表している。
【００５３】
アドレス時においては、ゴルフクラブ５は静止しているため、加速度センサ４２は、Ｚ軸方向の重力加速度のみを検出することになる。この事実から、アドレス時のｃ，ｆ，ｉは、下式に従って算出することができる。
【数２】
0006776882.tif 000003
【００５４】
また、残りのａ，ｂ，ｄ，ｅ，ｇ，ｈは、ｃ，ｆ，ｉを説明変数とし、ａ，ｂ，ｄ，ｅ，ｇ，ｈのそれぞれを目的変数とする予め定められた重回帰式より算出することができる。この重回帰式は、多数のａ〜ｉの真値のデータを重回帰分析することにより決定することができる。或いは、姿勢行列は、下式のとおり、オイラー変換行列としても表すことができる。従って、姿勢行列の９つの成分のうち３つの値が分かれば、φ，θ，ψを特定できるため、残りの成分を算出することができる。なお、姿勢行列の算出方法としては様々考えられ、いずれの方法を用いてもよいが、ここで説明した方法は、例えば、特開２０１３−５６０７４号公報においても開示されている。
【００５５】
【数３】
0006776882.tif 000004
【００５６】
導出部１４ｃは、アドレス時の姿勢行列Ｎの導出後、アドレス時の姿勢行列Ｎを初期値として、姿勢行列Ｎをサンプリング周期Δｔ間隔で時々刻々更新してゆく。具体的に説明すると、まず、姿勢行列Ｎは、オイラーパラメータｑ＝（ｑ0，ｑ1，ｑ2，ｑ3）を用いて、以下の式で表される。
【数４】
0006776882.tif 000005
【００５７】
従って、ある時刻でのａ〜ｉが既知の場合、当該時刻でのオイラーパラメータｑ＝（ｑ0，ｑ1，ｑ2，ｑ3）は、以下の式に従って算出することができる。
【数５】
0006776882.tif 000006
【００５８】
今、アドレス時の姿勢行列Ｎを規定するａ〜ｉの値は既知である。よって、以上の式に従って、まず、アドレス時のオイラーパラメータｑが算出される。
【００５９】
そして、時刻ｔから微小時刻経過後のオイラーパラメータｑt+1は、時刻ｔにおけるオイラーパラメータｑtを用いて以下の式で表される。
【数６】
0006776882.tif 000007
【００６０】
また、オイラーパラメータｑの時間変化を表す１階微分方程式は、以下の式で表される。
【数７】
0006776882.tif 000008
【００６１】
数６，７の式を用いれば、角速度データに基づいて、時刻ｔのオイラーパラメータｑを順次、次の時刻ｔ＋Δｔのオイラーパラメータｑへと更新することができる。ここでは、アドレスからインパクトまでのオイラーパラメータｑが算出される。そして、アドレスからインパクトまでのオイラーパラメータｑを数４の式に順次代入することにより、アドレスからインパクトまでの姿勢行列Ｎが算出される。
【００６２】
なお、慣性センサユニットから出力されるセンサデータに基づいて、オイラーパラメータｑや姿勢行列Ｎを算出する方法としては、ここで説明した方法に限らず、様々な方法が知られている。例えば、「根田康美，藤橋幸太，尾曲邦之，藤原謙，松永三郎"超小型衛星 Cute-1.7 シリーズにおける姿勢決定制御システムについて"第 16 回スペース・エンジニアリング・コンファレンス，東京，2007 年 1 月 25 日」に開示されている技術を利用することもできる。
【００６３】
続くステップＳ１３では、導出部１４ｃは、記憶部１３内の画像データに基づいて、アドレスからインパクト付近までの各時刻における時系列のシャフト５２の向きθs（図７参照）を導出する。具体的には、本実施形態では、ＩＲ画像を画像処理することにより、直線状のシャフト５２の像を検出する。そして、このシャフト５２の像が延びる方向とＺ方向との為す角度を、シャフト５２の向きθsとして導出する。
【００６４】
続くステップＳ１４では、導出部１４ｃは、ステップＳ１２の時系列の姿勢行列Ｎに基づいて、アドレスからインパクト付近までの各時刻における時系列のシャフト５２の向きθsを導出する。具体的には、θsは、以下の式に従って算出することができる。
【数８】
0006776882.tif 000009
【００６５】
次に、導出部１４ｃは、このセンサデータ由来の時系列のシャフト５２の向きθsと、ステップＳ１３の時系列のシャフト５２の向きθsとの一致度が最も高くなるように、両時系列データを位置合わせする。例えば、導出部１４ｃは、図８に示すように、両時系列データを時間軸に沿って前後に移動させながら、両時系列データの波形が最も一致する点を決定する。そして、そのときのタイミングを基準として、センサデータの時間軸と、画像データの時間軸とを相対的に時刻合わせさせる。
【００６６】
なお、本実施形態では、二次元のシャフト５２の向きθsに基づいて時刻合わせが行われたが、三次元のシャフト５２の向きθsに基づいて時刻合わせを行うこともできる。すなわち、姿勢行列Ｎからは、三次元的なシャフト５２の向き（ｇ，ｈ，ｉ）を導出することが可能であるし、深度画像を用いれば、画像データからも三次元的なシャフト５２の向きを導出することが可能である。従って、センサデータ及び画像データのそれぞれを由来とする三次元的なシャフト５２の向きの時系列データを導出し、これらの時系列データの波形が三次元的に最も一致するタイミングを決定すればよい。
【００６７】
＜１−３−３−２．アドレス期間の姿勢の導出（Ｓ４）＞
以下、図９を参照しつつ、アドレス期間の姿勢の導出処理の流れを説明する。ここで、オイラーパラメータｑは、上記のとおり、姿勢行列Ｎとの間で相互に変換可能である。また、上記のとおり、姿勢行列Ｎに含まれる行ベクトルは、慣性座標系内での局所座標系の三軸の向きを表している。従って、オイラーパラメータｑ及び姿勢行列Ｎの少なくとも一方が分かれば、慣性座標系内での局所座標系の三軸の向き、つまりはｚ軸に平行なシャフト５２の向き（ｇ，ｈ，ｉ）が分かる。そのため、以下では、シャフト５２の姿勢として、オイラーパラメータｑ及び姿勢行列Ｎが導出される。
【００６８】
まず、ステップＳ２１では、導出部１４ｃは、アドレス期間に含まれる各時刻の姿勢ｑ，Ｎを導出する。なお、以下のステップＳ２６では、ラグランジュの未定乗数法により、姿勢ｑが最適化される。この意味で、ステップＳ２１で導出される姿勢ｑ，Ｎは、仮の姿勢と言える。
【００６９】
具体的には、導出部１４ｃは、アドレス期間に含まれる最初の時刻である初期時刻ｔ1のオイラーパラメータｑを算出する。具体的には、導出部１４ｃは、初期時刻ｔ1のＩＲ画像及び深度画像からシャフト５２の向きを三次元的に導出する。すなわち、画像処理によりＩＲ画像からシャフト５２の像を検出し、当該像に含まれる各点のＹ座標及びＺ座標を導出する。そして、深度画像からこれらの点における深度情報を抽出し、Ｘ座標に変換する。その後、これらの三次元座標から三次元のシャフト５２の向き（ｇ，ｈ，ｉ）が導出される。また、導出部１４ｃは、初期時刻ｔ1の加速度ａx，ａy，ａzはほぼ重力のみを表しているという事実に基づいて、初期時刻ｔ1の加速度データ（ａx，ａy，ａz）から三次元の重力の向きを導出する。そして、導出部１４ｃは、シャフト５２の向き（ｇ，ｈ，ｉ）と、重力の向きの外積を導出し、これをｙ軸の方向（ｄ，ｅ，ｆ）とする。さらに、三次元のシャフト５２の向き（ｇ，ｈ，ｉ）と、ｙ軸の方向（ｄ，ｅ，ｆ）との外積を導出し、これをｘ軸の方向（ａ，ｂ，ｃ）とする。これにより、初期時刻ｔ1の姿勢行列Ｎが導出される。また、導出部１４ｃは、数５の式に従って、初期時刻ｔ1の姿勢行列Ｎを初期時刻ｔ1のオイラーパラメータｑに変換する。
【００７０】
また、ステップＳ２１では、導出部１４ｃは、初期時刻ｔ1の姿勢行列Ｎを初期値として、姿勢行列Ｎをサンプリング周期Δｔ間隔で時々刻々更新してゆく。具体的な方法は、上述のステップＳ１２と同様であるため、ここでは詳細な説明を省略する。
【００７１】
続くステップＳ２２では、導出部１４ｃは、初期時刻ｔ1の局所座標系での地磁気ｍB＝（ｍBx，ｍBy，ｍBz）＝（ｍx，ｍy，ｍz）を、慣性座標系での地磁気ｍI＝（ｍIx，ｍIy，ｍIz）に変換する。なお、ｍI，ｍBに限らず、右肩にＢを付けたときは、局所座標系のベクトルであることを表し、右肩にＩを付けた場合には、慣性座標系のベクトルであることを表す。ｍIは、時刻によらず一定とすることができ、オイラーパラメータｑを用いて、初期時刻ｔ1のｍBから下式に従って算出することができる。なお、以下のｑには、ステップＳ２１で導出された初期時刻ｔ1のオイラーパラメータｑを代入する。
【数９】
0006776882.tif 000010
ただし、
【数１０】
0006776882.tif 000011
である。
【００７２】
地磁気ｍIは、別の方法によっても導出可能である。例えば、磁気センサ４３とは異なる磁気センサを、慣性座標系の原点において、当該磁気センサの三軸を慣性座標系の三軸に一致させた状態で配置する。この場合、当該磁気センサの出力値は、慣性座標系での地磁気ｍIを表すことになる。
【００７３】
続くステップＳ２３では、導出部１４ｃは、アドレス期間に含まれる任意の時刻ｔi，ｔjの組み合わせに対し、時刻ｔiでの姿勢ｑiと、時刻ｔjにおける姿勢ｑjとの差である、ｑBjBiを導出する。このとき、ｔj，ｔiの先後は問わない。姿勢の差ｑBjBiは、下式に従って算出することができる。なお、下式に代入されるオイラーパラメータｑi，ｑjは、ステップＳ２１で算出されたオイラーパラメータｑである。
【数１１】
0006776882.tif 000012
【００７４】
姿勢の差ｑBjBiは、別の方法によっても導出可能である。例えば、数１１の式に代入すべきアドレス期間に含まれる各時刻のオイラーパラメータｑを、ステップＳ２１で算出したものではなく、以下の方法で算出することができる。すなわち、上記ステップＳ２１では、初期時刻ｔ1より後の各時刻のオイラーパラメータｑは、角速度データを用いて初期時刻ｔ1のオイラーパラメータｑを順次更新してゆくステップＳ１２と同様の方法により算出された。しかしながら、初期時刻ｔ1より後の各時刻の姿勢行列Ｎについても、初期時刻ｔ1の場合と同様に、画像データ及び加速度データを用いて外積を計算する方法により算出してもよい。
【００７５】
あるいは、アドレス期間における姿勢の差ｑBjBiは、単位クォータニオン（１，０，０，０）と近似することができる。これは、アドレス期間においては、グリップ５１が静止しているか又は殆ど動いておらず、姿勢の差が殆どないためである。
【００７６】
続くステップＳ２４では、導出部１４ｃは、アドレス期間に含まれる任意の時刻ｔi，ｔjの組み合わせに対し、以下に定義される行列Ａkjを導出する。Ａkjは、後のステップＳ２６の最適化処理の中で、目的関数に代入される。
【数１２】
0006776882.tif 000013
【００７７】
以下、Ａkj及び上式中の右辺の各記号ついて説明する。まず、時刻ｔi，ｔjにおける慣性座標系での三次元のシャフト５２の向きを表すベクトルをそれぞれｂIi、ｂIjとし、時刻ｔi，ｔjにおける局所座標系での三次元のシャフト５２の向きを表すベクトルをそれぞれｂBi，ｂBjとする。なお、ｂBi，ｂBjは、時刻によらず一定値［０  ０  −１］Tとすることができ（上付きＴは転置を表す）、従って、ｂBi，ｂBjはいずれも、単にｂBと表記することができる。このとき、下式が成り立つ。
【数１３】
0006776882.tif 000014
【００７８】
また、数１１より、下式が導かれる。
【数１４】
0006776882.tif 000015
【００７９】
ここで、数１３，１４の式を変形すると、以式が成り立つ。
【数１５】
0006776882.tif 000016
ただし、
【数１６】
0006776882.tif 000017
である。数１５の式を変形すると、下式のようになる。
【数１７】
0006776882.tif 000018
【００８０】
数１７の式からは、ｑiに関して線形の関係が成り立っていることが分かる。この式を整理すると、以下の関係式が導かれる。
【数１８】
0006776882.tif 000019
【００８１】
数１８の関係式は、画像データやセンサデータに誤差がなければ、常に成り立つ恒等式である。数１８の関係式は、ステップＳ２６の最適化に用いられる目的関数を定義するのに用いられる。
【００８２】
ここで、数１２の式に戻ると、ｂIjx，ｂIjy，ｂIjzは、時刻ｔjにおける慣性座標系での三次元のシャフト５２の向きを表すｂIjのＸ，Ｙ，Ｚ軸成分を意味する。また、ｂjix，ｂjiy，ｂjizは、数１６により定義されるｂjiのＸ，Ｙ，Ｚ軸成分を意味する。なお、記号のハットは、正規化されていることを意味する。
【００８３】
ステップＳ２４では、導出部１４ｃは、アドレス期間のＩＲ画像及び深度画像に基づいて、アドレス期間の各時刻における三次元のシャフト５２の向きｂIjを導出する。そして、Ａkjを導出すべく、その三軸成分であるｂIjx，ｂIjy，ｂIjzの値を、数１２の式に代入する。また、シャフト５２の向きを表すｂB＝（ｂjx，ｂjy，ｂjz）＝［０  ０  −１］Tと、ステップＳ２３のｑBjBiを数１６の式に代入することにより、ｂjiを導出し、さらにこれを数１２に代入する。
【００８４】
続くステップＳ２５では、導出部１４ｃは、アドレス期間に含まれる任意の時刻ｔi，ｔjの組み合わせに対し、以下に定義される行列Ａmjを導出する。Ａmjは、Ａkjと同様に、後のステップＳ２６の最適化処理の中で目的関数に代入される。
【数１９】
0006776882.tif 000020
【００８５】
ステップＳ２４の中で説明した場合と同様に、下式が導かれる。この関係式は、ステップＳ２６の最適化に用いられる目的関数を定義するのに用いられる。
【数２０】
0006776882.tif 000021
【００８６】
ここで、数１９の式に戻ると、ｍIjx，ｍIjy，ｍIjzは、時刻ｔjにおける慣性座標系での地磁気を表すｍIjのＸ，Ｙ，Ｚ軸成分を意味する。また、ｍjix，ｍjiy，ｍjizは、下式により定義されるｍjiのＸ，Ｙ，Ｚ軸成分を意味する。
【数２１】
0006776882.tif 000022
【００８７】
ステップＳ２５では、導出部１４ｃは、数１９の右辺において、ｍIjとしてステップＳ２２で導出された地磁気のデータを代入する。また、時刻ｔjにおける地磁気のセンサデータであるｍBj＝（ｍjx，ｍjy，ｍjz）と、ステップＳ２３のｑBjBiを数２１の式に代入することにより、ｍjiを導出し、さらにこれを数１９の右辺に代入する。
【００８８】
続くステップＳ２６では、ラグランジュの未定乗数法により、アドレス期間の各時刻ｔiにおける姿勢ｑiが最適化される。具体的には、各時刻ｔiにおいて以下の目的関数Ｊ1を定義する。
【数２２】
0006776882.tif 000023
ただし、
【数２３】
0006776882.tif 000024
である。
【００８９】
ここで、数１８，２０から明らかなとおり、Ｆk，Ｆmは、画像データやセンサデータに誤差がなければ０になるべき関数である。また、（１−ｑiTｑi）は、拘束条件である。そこで、導出部１４ｃは、ステップＳ２４，Ｓ２５で計算したＡkj，Ａmjを用いて、目的関数Ｊ1を最小化するようなオイラーパラメータｑiの最適解を導出する。また、導出部１４ｃは、各時刻ｔiにおける最適化後のオイラーパラメータｑiを数４の式に従って、姿勢行列Ｎに変換する。
【００９０】
なお、ｗk，ｗmは、画像データ、地磁気データのいずれをどれだけ信頼するかに応じて、設定される重み係数である。
【００９１】
また、本実施形態では、目的関数Ｊ1の最小化は、以下のとおり行われる。すなわち、目的関数Ｊ1を最小化するためには、その微分値（偏微分値）がゼロになること、言い換えると、下式が成立することが必要条件となる。ただし、Ｉ4は、４×４の単位行列である。
【数２４】
0006776882.tif 000025
【００９２】
上式は、下式の行列の固有値がλとなり、ｑiが固有ベクトルになることを表している。
【数２５】
0006776882.tif 000026
【００９３】
従って、導出部１４ｃは、ステップＳ２４，Ｓ２５の結果から、数２５の行列を計算し、当該行列の最小の固有値に対応する大きさ１の固有ベクトルを導出して、これをｑiとする。
【００９４】
＜１−３−３−３．アドレスからインパクト付近までの姿勢の導出（Ｓ５）＞
以下、図１０を参照しつつ、アドレスからインパクト付近までの姿勢の導出処理の流れを説明する。本処理では、アドレス期間における姿勢の導出処理と類似のアルゴリズムが用いられる。両処理の主な相違点は、最適化のための目的関数Ｊ2に、上記関数Ｆk，Ｆmに加え、画像データ由来の二次元のシャフト５２の向きに基づく関数Ｆnが含まれる点にある。以下、詳細に説明する。
【００９５】
ステップＳ３１では、導出部１４ｃは、アドレスからインパクト付近までの各時刻の仮の姿勢ｑ，Ｎを導出する。具体的には、導出部１４ｃは、角速度ωx，ωy，ωzのデータを用いて、数６，７の式に従って、最適化後のアドレス時のオイラーパラメータｑを初期値としてオイラーパラメータｑをサンプリング周期Δｔ間隔で時々刻々更新してゆく。その後、数５の式に従って、各時刻の姿勢ｑを、姿勢行列Ｎに換算する。
【００９６】
続くステップＳ３２では、導出部１４ｃは、アドレス時の局所座標系での地磁気ｍB＝（ｍBx，ｍBy，ｍBz）＝（ｍx，ｍy，ｍz）を、慣性座標系での地磁気ｍI＝（ｍIx，ｍIy，ｍIz）に変換する。具体的な変換方法は、ステップＳ２２と同様である。ただし、ここで変換に用いられるオイラーパラメータｑは、ステップＳ３１で算出されたオイラーパラメータｑである。
【００９７】
続くステップＳ３３では、導出部１４ｃは、アドレスからインパクト付近までの期間に含まれる任意の時刻ｔi，ｔjの組み合わせに対し、時刻ｔiでの姿勢ｑiと、時刻ｔjにおける姿勢ｑjとの差である、ｑBjBiを導出する。このとき、ｔj，ｔiの先後は問わない。姿勢の差ｑBjBiは、数１１の式に従って算出することができる。なお、ここで数１１の式に代入されるオイラーパラメータｑi，ｑjは、ステップＳ３１で算出されたオイラーパラメータｑである。
【００９８】
続くステップＳ３４では、導出部１４ｃは、アドレスからインパクト付近までの期間に含まれる任意の時刻ｔi，ｔjの組み合わせに対し、数１２のとおり定義される行列Ａkjを導出する。Ａkjは、後のステップＳ３７の最適化処理の中で、目的関数に代入される。
【００９９】
ステップＳ３４では、導出部１４ｃは、アドレスからインパクト付近までの期間のＩＲ画像及び深度画像に基づいて、当該期間の各時刻における三次元のシャフト５２の向きｂIjを導出する。そして、Ａkjを導出すべく、その三軸成分であるｂIjx，ｂIjy，ｂIjzの値を、数１２の式に代入する。また、シャフト５２の向きを表すｂB＝（ｂjx，ｂjy，ｂjz）＝［０  ０  −１］Tと、ステップＳ３３のｑBjBiを数１６の式に代入することにより、ｂjiを導出し、さらにこれを数１２に代入する。
【０１００】
続くステップＳ３５では、導出部１４ｃは、アドレスからインパクト付近までの期間に含まれる任意の時刻ｔi，ｔjの組み合わせに対し、数１９のとおり定義される行列Ａmjを導出する。Ａmjは、後のステップＳ３７の最適化処理の中で、目的関数に代入される。
【０１０１】
ステップＳ３５では、導出部１４ｃは、数１９の右辺において、ｍIjとしてステップＳ３２で導出された地磁気のデータを代入する。また、時刻ｔjにおける地磁気のセンサデータであるｍBj＝（ｍjx，ｍjy，ｍjz）と、ステップＳ３３のｑBjBiを数２１の式に代入することにより、ｍjiを導出し、さらにこれを数１９の右辺に代入する。
【０１０２】
続くステップＳ３６では、導出部１４ｃは、アドレスからインパクト付近までの期間に含まれる任意の時刻ｔi，ｔjの組み合わせに対し、以下に定義される行列Ａnjを導出する。Ａnjは、後のステップＳ３７の最適化処理の中で、目的関数に代入される。
【数２６】
0006776882.tif 000027
【０１０３】
以下、Ａnj及び上式中の右辺の各記号ついて説明する。まず、ここでは、ＹＺ平面内における二次元のシャフト５２の向きを考える。ここで、ＸＹＺ座標系における三次元座標をＹＺ平面に射影するための射影行列をＰYZとする。このとき、時刻ｔiにおけるＹＺ平面内における二次元のシャフトの向きｂSiは、ＰYZｂIiと表される。ｂSiは、オイラーパラメータｑiを用いて、以下のとおり表すことができる。ただし、ｕ＝[１  ０  ０]Tである。
【数２７】
0006776882.tif 000028
【０１０４】
数２７を変形すると、以下のとおりとなる。
【数２８】
0006776882.tif 000029
さらに、数１４の式から、数２８を変形する。
【数２９】
0006776882.tif 000030
【０１０５】
さらに、数２９の両辺の右からｑBjBiを掛ける。
【数３０】
0006776882.tif 000031
さらに、数１３，１６から、以下のように変形される。
【数３１】
0006776882.tif 000032
【０１０６】
数３１の式からは、ｑiに関して線形の関係が成り立っていることが分かる。この式を整理すると、以下の関係式が導かれる。
【数３２】
0006776882.tif 000033
【０１０７】
数３２の関係式は、画像データやセンサデータに誤差がなければ、常に成り立つ恒等式である。数３２の関係式は、ステップＳ３７の最適化に用いられる目的関数を定義するのに用いられる。
【０１０８】
ここで、数２６の式に戻ると、ｂSjx，ｂSjy，ｂSjzは、時刻ｔjにおけるＹＺ平面内における二次元のシャフトの向きを表すｂSjのＸ，Ｙ，Ｚ軸成分を意味する。従って、Ｘ成分であるｂSjx＝０である。
【０１０９】
ステップＳ３６では、導出部１４ｃは、Ａnjを導出すべく、ｂB＝［０  ０  −１］Tを数１３を用いて慣性座標系の値（ｂIjx，ｂIjy，ｂIjz）に変換し、これを数２６の式に代入する。また、シャフト５２の向きを表すｂB＝（ｂjx，ｂjy，ｂjz）＝［０  ０  −１］Tと、ステップＳ３３のｑBjBiを数１６の式に代入することにより、ｂjiを導出し、さらにこれを数２６に代入する。さらに、導出部１４ｃは、ＩＲ画像から導出されるｂSjを、数２６の式に代入する。
【０１１０】
続くステップＳ３７では、ラグランジュの未定乗数法により、アドレスからインパクト付近までの各時刻ｔiにおける姿勢ｑiが最適化される。具体的には、各時刻ｔiにおいて以下の目的関数Ｊ2を定義する。
【数３３】
0006776882.tif 000034
ただし、
【数３４】
0006776882.tif 000035
である。
【０１１１】
ここで、数３２から明らかなとおり、Ｆnは、Ｆk，Ｆmと同じく、画像データやセンサデータに誤差がなければ０になるべき関数である。また、（１−ｑiTｑi）は、拘束条件である。そこで、導出部１４ｃは、ステップＳ３４〜Ｓ３６で計算したＡkj，Ａmj，Ａnjを用いて、目的関数Ｊ2を最小化するようなオイラーパラメータｑiの最適解を導出する。また、導出部１４ｃは、各時刻ｔiにおける最適化後のオイラーパラメータｑiを数４の式に従って、姿勢行列Ｎに変換する。
【０１１２】
なお、ｗk，ｗm，ｗnは、深度画像のデータ、地磁気データ、ＩＲ画像のデータのいずれをどれだけ信頼するかに応じて、設定される重み係数である。
【０１１３】
また、本実施形態では、目的関数Ｊ2の最小化は、以下のとおり行われる。すなわち、目的関数Ｊ2を最小化するためには、その微分値（偏微分値）がゼロになること、言い換えると、下式が成立することが必要条件となる。
【数３５】
0006776882.tif 000036
【０１１４】
上式は、下式の行列の固有値がλとなり、ｑiが固有ベクトルになることを表している。
【数３６】
0006776882.tif 000037
【０１１５】
従って、導出部１４ｃは、ステップＳ３３〜Ｓ３６の結果から、数３６の行列を計算し、当該行列の最小の固有値に対応する大きさ１の固有ベクトルを導出して、これをｑiとする。
【０１１６】
＜１−３−３−４．グリップの位置及び速度の導出（Ｓ６）＞
以下、図１１を参照しつつ、アドレスからインパクト付近までのグリップ５１の位置及び速度の導出処理の流れを説明する。
【０１１７】
グリップ５１の位置は、画像データにより計測される。しかしながら、距離画像センサ２のサンプリング周期が粗い場合や、グリップ５１が画像上に写らない場合等には、位置の測定が困難になることがある。そこで、ここでは、カルマンフィルタを用いる等して、画像データ由来のグリップ５１の位置が加速度データによりスムージングされる。
【０１１８】
まず、ステップＳ４１では、導出部１４ｃは、下式に従って、アドレスからインパクト付近までの各時刻ｔiにおける局所座標系の加速度ａBi＝（ａx，ａy，ａz）を慣性座標系の加速度に変換し、さらに重力加速度の成分を除去する。これにより、慣性座標系でのグリップ５１の加速度ａIiが導出される。なお、下式のオイラーパラメータｑiには、ステップＳ３７で導出された値を代入することができる。Ｇは、重力加速度の大きさである。
【数３７】
0006776882.tif 000038
【０１１９】
続くステップＳ４２，Ｓ４３では、導出部１４ｃは、ステップＳ４１の加速度ａIiを用いて、ＩＲ画像及び深度画像から導出されるグリップ５１の三次元の位置ｐをスムージングすることにより、誤差の除去された位置ｐ及び速度ｐ'（推定値）を導出する。この意味で、ＩＲ画像及び深度画像から導出されるグリップ５１の三次元の位置ｐは、それぞれ仮の位置ということができる。
【０１２０】
グリップ５１の位置ｐ及び速度ｐ'を推定するためのシステムは、以下のようにモデル化することができる。ただし、ｘは、求めたい位置ｐ及び速度ｐ'、ｙは、ＩＲ画像及び深度画像から導出される位置ｐ、ｕは、加速度データ、ｗは、プロセス雑音、ｖは、観測雑音を表す。添え字のｋは、時系列のデータ番号を表し、ｋ＝１，２，・・・，Ｎ（Ｎは、データ個数）である。なお、ステップＳ４２，Ｓ４３の説明において登場する符号／記号は、特に断らない限り、これまでの説明及びこれよりも後の説明に登場する符号／記号とは独立して定義されているものとする。
【数３８】
0006776882.tif 000039
ただし、
【数３９】
0006776882.tif 000040
とする。
【０１２１】
また、プロセス雑音ｗkの平均値を０、分散をＱkとし、観測雑音ｖkの平均値を０、分散をＲkとする。このＱk及びＲkは、スムージングをかける程度を決定する因子となり、Ｑk／Ｒkが大きければあまりスムージングが行われず、距離画像センサ２の計測結果の比重が大きくなる。一方、Ｑk／Ｒkが小さければ、スムージングが強くなり、加速度センサ４１の計測結果の比重が大きくなる。
【０１２２】
今、時間軸に沿って前向きにフィルタリングを行う場合の推定値をｘ+fk、誤差共分散をＰ+fkとする。導出部１４ｃは、ステップＳ４２において、ｋ＝０のときの初期値を以下のとおり導出する。なお、Ｅは、期待値を表す。
【数４０】
0006776882.tif 000041
【０１２３】
その後、導出部１４ｃは、この初期値を以下の式に従って、更新してゆく。
【数４１】
0006776882.tif 000042
【０１２４】
なお、ｘ+fkは、時刻ｋの状態の時刻ｋ時点での推定値を表し、ｘ-fkは、時刻ｋの状態の時刻（ｋ−１）時点での推定値を表す。誤差共分散についても同様に、Ｐ+fkは、時刻ｋの状態の時刻ｋ時点での推定値を表し、Ｐ-fkは、時刻ｋの状態の時刻（ｋ−１）時点での推定値を表す。また、ｘ-fk、Ｐ-fkは、以下のとおり算出することができる。
【数４２】
0006776882.tif 000043
である。
【０１２５】
以上により、過去及び現在の加速度データ及び画像データから誤差の除去されたｘ＝（ｐ，ｐ'）が推定される。なお、本実施形態では、画像データのサンプリング周期は、慣性センサユニット４の出力値である加速度データのサンプリング周期ΔＴよりも大きい。そのため、ステップＳ４２では、画像データを、線形補完等を適宜行って、加速度データのサンプリング周期ΔＴに合わせたデータセットにしておく前処理が行われる。
【０１２６】
続くステップＳ４３では、導出部１４ｃは、時間軸に沿って後向きにフィルタリングを行うことで、未来の加速度データ及び画像データを用いて、グリップ５１の位置ｐ及び速度ｐ'をさらに更新する。
【０１２７】
具体的には、ｋ＝Ｎ−１，・・・，１，０に対して、以下の式に従ってスムージングを行うことにより、ｘkを導出する。なお、初期値となるｘN＝ｘ+fNである。
【数４３】
0006776882.tif 000044
【０１２８】
以上により、過去及び現在の計測結果のみならず、未来の計測結果を用いて、高精度にグリップ５１の位置ｐ及び速度ｐ'を導出することができる。
【０１２９】
＜１−３−３−５．グリップの角速度及び角加速度の導出（Ｓ７）＞
以下、図１２を参照しつつ、アドレスからインパクト付近までのグリップ５１の角速度及び角加速度の導出処理の流れを説明する。
【０１３０】
まず、ステップＳ５１において、導出部１４ｃは、ステップＳ５で導出された姿勢ｑiに疑似微分フィルタを掛けることにより、アドレスからインパクト付近までの各時刻ｔiの姿勢の微分ｑi'を導出する。なお、疑似微分フィルタは、以下の式で表される。
【０１３１】
【数４４】
0006776882.tif 000045
【０１３２】
ただし、本実施形態では、数４４の式を双１次変換して離散化した以下の式を用いて、時間微分が実行される。
【数４５】
0006776882.tif 000046
【０１３３】
続くステップＳ５２では、導出部１４ｃは、ステップＳ５で導出された姿勢ｑiと、ステップＳ５１で導出された姿勢の微分ｑi'とを用いて、アドレスからインパクト付近までの各時刻ｔiの角速度ωiを導出する。角速度ωiは、下式に従って、導出される。
【数４６】
0006776882.tif 000047
【０１３４】
続くステップＳ５３では、導出部１４ｃは、ステップＳ５２で導出した角速度ωiに疑似微分フィルタを掛けることにより、アドレスからインパクト付近までの各時刻ｔiの角加速度ωi'を算出する。疑似微分フィルタは、数４５と同じである。
【０１３５】
＜１−３−３−６．グリップの加速度の導出（Ｓ８）＞
以下、図１３を参照しつつ、アドレスからインパクト付近までのグリップ５１の加速度の導出処理の流れを説明する。
【０１３６】
まず、ステップＳ６１において、導出部１４ｃは、ステップＳ６で導出されたスムージング後のグリップ５１の速度ｐ'に疑似微分フィルタを掛けることにより、アドレスからインパクト付近までの各時刻のグリップ５１の加速度ｐ''を導出する。疑似微分フィルタは、数４５と同じである。
【０１３７】
続くステップＳ６２において、導出部１４ｃは、ステップＳ５で導出されたオイラーパラメータｑiを用いて、以下の式に従って、アドレスからインパクト付近までの各時刻におけるステップＳ６１の慣性座標系での加速度ｐ''＝ａIi、を、局所座標系の加速度ａBiに変換する。
【数４７】
0006776882.tif 000048
【０１３８】
＜１−４．評価＞
以下、本発明の実施例について説明するが、本発明は、以下の実施例に限定されない。
【０１３９】
＜１−４−１＞
実施例１として、ステップＳ４の最適化の効果を確認すべく、上記実施形態に係るステップＳ１〜Ｓ４を実行して得られた最適化された姿勢行列Ｎを用いて、局所座標系のシャフトの向き[０，０，１]Tを全体座標系でのシャフトの向きに変換した。そして、このシャフトの向きのベクトルと、モーションキャプチャーシステムで計測した全体座標系でのシャフトの向きを表すベクトル（真値）との為す角度を、角度誤差として算出した。この結果を、図１４に実線で示す。なお、ここでいうモーションキャプチャーとは、例えば、特開2011−183090号公報及び特開2011−030761号公報に記載されるような、多数の赤外線カメラを使って高精度で動きを計測するシステムである。また、参考例１として、ステップＳ１〜Ｓ３及びステップＳ２１を実行して姿勢行列Ｎを導出し、この姿勢行列Ｎを用いて、局所座標系のシャフトの向き[０，０，１]Tを全体座標系でのシャフトの向きに変換した。そして、このシャフトの向きのベクトルと、モーションキャプチャーシステムで計測した全体座標系でのシャフトの向きを表すベクトル（真値）との為す角度を、角度誤差として算出した。この結果を、図１４に点線で示す。同図からは、ステップＳ４の最適化処理を経ることにより、真値との誤差が小さくなっていることが分かる。
【０１４０】
＜１−４−２＞
実施例２として、ステップＳ５の最適化の効果を確認すべく、上記実施形態に係るステップＳ１〜Ｓ５を実行して得られた最適化された姿勢行列Ｎを用いて、局所座標系のシャフトの向き[０，０，１]Tを全体座標系でのシャフトの向きに変換した。そして、このシャフトの向きのベクトルと、モーションキャプチャーシステムで計測した全体座標系でのシャフトの向きを表すベクトル（真値）との為す角度を、角度誤差として算出した。この結果を、図１５に実線で示す。また、参考例２として、ステップＳ１〜Ｓ４及びステップＳ３１を実行して姿勢行列Ｎを導出し、この姿勢行列Ｎを用いて、局所座標系のシャフトの向き[０，０，１]Tを全体座標系でのシャフトの向きに変換した。そして、このシャフトの向きのベクトルと、モーションキャプチャーシステムで計測した全体座標系でのシャフトの向きを表すベクトル（真値）との為す角度を、角度誤差として算出した。この結果を、図１５に点線で示す。同図からは、ステップＳ５の最適化処理を経ることにより、真値との誤差が小さくなっていることが分かる。
【０１４１】
＜１−４−３＞
実施例３として、ステップＳ６のスムージングの効果を確認すべく、上記実施形態に係るステップＳ１〜Ｓ６を実行して、グリップの位置ｐを算出した。そして、この位置ｐと真値の位置ｐとの差（位置誤差）を算出した結果を、図１６Ａ〜図１６Ｃに実線で示す。なお、真値は、モーションキャプチャーで計測した値とした。また、参考例３として、ステップＳ１により得られた同じセンサデータから、ステップＳ１１，Ｓ１２と同様の方法で導出された姿勢行列Ｎを用いて絶対座標系の加速度データを算出し、ここから重力成分を除去した後、これを２回積分することによりグリップの位置ｐを算出した。そして、この位置ｐと真値の位置ｐとの差（位置誤差）を算出した結果を、図１６Ａ〜図１６Ｃに点線で示す。同図からは、ステップＳ６のスムージングを経ることにより、真値との誤差が小さくなっていることが分かる。
【０１４２】
＜１−４−４＞
実施例４として、ステップＳ４３のスムージングの効果を確認すべく、上記実施形態に係るステップＳ１〜Ｓ６を実行して、グリップの位置ｐを算出した。そして、この位置ｐと真値の位置ｐとの差（位置誤差）を算出した結果を、図１７Ａ〜図１７Ｃに点線で示す。なお、真値は、モーションキャプチャーで計測した値とした。また、参考例４として、ステップＳ１〜Ｓ５，Ｓ４１，Ｓ４２を実行して、グリップの位置ｐを算出した。そして、この位置ｐと真値の位置ｐとの差（位置誤差）を算出した結果を、図１７Ａ〜図１７Ｃに実線で示す。同図からは、ステップＳ４３のスムージングを経ることにより、真値との誤差が小さくなっていることが分かる。
【０１４３】
＜２．第２実施形態＞
図１８及び図１９に、本発明の第２実施形態に係る運動解析装置１０１を含む運動解析システム２００の全体構成図を示す。これらの図を図１及び図２と比較すれば明らかなように、第２実施形態の運動解析システム２００には、第１実施形態の運動解析システム１００に対し、２台目の距離画像センサ６が追加されている。２台目の距離画像センサ６は、図１９に示すとおり、１台目の距離画像センサ２と同様の構成を有している。距離画像センサ６は、距離画像センサ２とは異なる方向から撮影すべく、具体的には、ゴルファー７を右側から撮影すべく、ゴルファー７の右側に設置される。距離画像センサ６も、通信線１７を介して運動解析装置１０１に接続される。なお、運動解析装置１０１は、複数台のコンピュータから構成されていてもよく、例えば、距離画像センサ２，６が異なるコンピュータに接続されていてもよい。
【０１４４】
第２実施形態の運動解析装置１０１と、第１実施形態の運動解析装置１とは、ハードウェア構成は同じである。しかしながら、後述するとおり、一部ソフトウェア処理が異なるため、異なる運動解析プログラム１０３がインストールされている。それにより、制御部１４は、第１及び第２取得部１４ａ，１４ｂとして動作する他、第３取得部１４ｅとしても動作し、導出部１４ｃに代えて、導出部１１４ｃとして動作する。以下では、第１及び第２実施形態間で共通の要素には共通の符号を付し、適宜説明を省略しつつ、主として両実施形態の相違点について説明する。
【０１４５】
なお、本明細書では様々な記号が登場し、同じ記号が異なる意味で、或いは異なる記号が同じ意味で用いられる場合があるが、最新の定義が優先されるものとする。
【０１４６】
＜２−１．運動解析方法＞
図２０を参照しつつ、第２実施形態に係る運動解析方法について説明する。図５及び図２０を比較すれば明らかなとおり、第１及び第２実施形態に係る運動解析方法は類似しているが、ステップＳ２，Ｓ３，Ｓ６が、ステップＳ１０２，Ｓ１０３，Ｓ１０６に改良されている点で異なる。また、第２実施形態では、ステップＳ５の後に、ステップＳ５で導出されたシャフト５２の姿勢を補正するステップＳ１０５が挿入されている。以下、相違点に係るステップＳ１０２，Ｓ１０３，Ｓ１０５，Ｓ１０６について順に説明する。
【０１４７】
＜２−１−１．第２データ収集処理（Ｓ１０２）＞
第２実施形態では、第２データ収集処理として、ゴルフクラブ５のスイング動作が２台の距離画像センサ２，６により撮影される。すなわち、距離画像センサ２だけでなく、同時に距離画像センサ６も、ゴルフスイングを捉えた時系列のＩＲ画像及び深度画像を含む画像データを検出する。距離画像センサ６により検出された画像データは、距離画像センサ２により検出された画像データと同様に、距離画像センサ６の通信部２５を介して運動解析装置１０１に送信される。一方、運動解析装置１０１側では、第３取得部１４ｅが通信部１５を介してこれを受信し、記憶部１３内に格納する。本実施形態では、少なくともアドレスの少し前の初期時刻ｔ1からフィニッシュまでの時系列の画像データが収集される。
【０１４８】
＜２−１−２．時刻合わせ（Ｓ１０３）＞
本実施形態では、２台の距離画像センサ２，６にそれぞれ由来する２系列の画像データが取得される。そのため、第１実施形態では、画像データとセンサデータとの時刻合わせのみが行われたが、本実施形態では、ステップＳ１０３として、これら２系列の画像データの時刻合わせも行われる。
【０１４９】
具体的には、ステップＳ１０３は、図２１に示すフローチャートに従って進行する。同図に示されるとおり、ステップＳ１０３では、導出部１１４ｃにより、ステップＳ３と同じくステップＳ１１〜Ｓ１４が実行される。これにより、センサデータと、距離画像センサ２により正面から撮影された画像データ（以下、正面画像データという）との同期が取られる。続いて、ステップＳ１１１,Ｓ１１２において、正面画像データと、距離画像センサ６により右側から撮影された画像データ（以下、右側画像データという）との同期が取られる。まず、ステップＳ１１１では、導出部１１４ｃが、正面画像データを画像処理することにより、アドレスからインパクト付近までの各時刻におけるグリップ端５１ａの三次元座標を導出する。同様に、導出部１１４ｃが、右側画像データを画像処理することにより、アドレスからインパクト付近までの各時刻におけるグリップ端５１ａの三次元座標を導出する。
【０１５０】
続くステップＳ１１２では、導出部１１４ｃは、正面画像データ由来の時系列のグリップ端５１ａの三次元座標と、右側画像データ由来の時系列のグリップ端５１ａの三次元座標との一致度が最も高くなるように、両時系列データを位置合わせする。つまり、ステップＳ１４と同様に、両時系列データを時間軸に沿って前後に移動させながら、両時系列データの波形が最も一致する点を決定する。そして、そのときのタイミングを基準として、正面画像データの時間軸と、右側画像データの時間軸とを相対的に時刻合わせさせる。
【０１５１】
以上により、時刻合わせは終了する。なお、以上の処理により、正面画像データを介して、右側画像データとセンサデータも時刻合わせされる。
【０１５２】
＜２−１−３．ダウンスイング期間の姿勢の補正（Ｓ１０５）＞
図１５を参照されたい。同図は、上述したとおり、ステップＳ１〜Ｓ５により導出された姿勢行列Ｎに基づくシャフト５２の向きの角度誤差を表している。同図から分かるように、この誤差は、インパクト（図中、０秒に対応する）の少し前までは効果的に低減されているが、インパクト直前においてはやや大きい。これは、ダウンスイング期間に概ね相当するインパクト直前においてはゴルフクラブ５が高速に運動し、角速度センサ４２（ジャイロセンサ）のバイアス成分が大きくなるためと考えられる。
【０１５３】
ステップＳ１０５は、このような誤差が比較的大きくなるインパクト直前の期間、言い換えると、（インパクト−Ｔ0秒）からインパクトまでの期間（以下、ダウンスイング期間という）におけるシャフト５２の姿勢を補正するステップである。Ｔ0は、例えば、０．１秒〜０．５秒と設定することができ、０．２秒〜０．３秒がより好ましい。
【０１５４】
ステップＳ１０５は、図２２に示すとおり進行する。まず、ステップＳ１２１では、導出部１１４ｃは、画像データに基づいてインパクト時の姿勢ｑ（オイラーパラメータ）を導出する。より具体的には、導出部１１４ｃは、右側画像データに含まれるインパクト時の深度画像を画像処理することにより、シャフト５２の向きを表す三次元ベクトルを導出するとともに、右側画像データに含まれるインパクト時のＩＲ画像を画像処理することにより、シャフト５２の向きを表す二次元ベクトルを導出する。同様に、導出部１１４ｃは、正面画像データに含まれるインパクト時の深度画像を画像処理することにより、シャフト５２の向きを表す三次元ベクトルを導出するとともに、正面画像データに含まれるインパクト時のＩＲ画像を画像処理することにより、シャフト５２の向きを表す二次元ベクトルを導出する。さらに、導出部１１４ｃは、右側画像データ及び正面画像データに含まれるインパクト時の対のＩＲ画像を画像処理することにより、シャフト５２の向きを表す三次元ベクトルを導出する。このとき、ＤＬＴ法等を用いて、２枚のＩＲ画像に由来する２方向からのシャフト５２の像の２次元座標に基づいて、シャフト５２の像の３次元座標が導出される。
【０１５５】
導出部１１４ｃは、以上の５つのシャフト５２の向きと、未知の姿勢ｑとの誤差の大きさを表わす目的関数を定義する。そして、当該目的関数を最適化することにより、当該誤差を最小にする姿勢ｑを導出する。なお、本実施形態では、シャフト５２の向きを表す上記５つのベクトルを全て用いて目的関数が定義されるが、例えば、これらのうちの２つ以上のベクトルを任意に組み合わせて、目的関数を定義することもできる。
【０１５６】
続くステップＳ１２２では、導出部１１４ｃは、ステップＳ１２１で導出されたインパクト時の姿勢ｑ（以下、ｑIBsと表記する）に基づいて、ステップＳ５で導出されたインパクト時の姿勢ｑ（以下、ｑIBbと表記する）を補正する。ｑIBsは、Ｂs系の慣性系（ＸＹＺ慣性座標系に一致する。以下同じ）に対する姿勢を表すオイラーパラメータであり、ｑIBbは、Ｂb系の慣性系に対する姿勢を表すオイラーパラメータである。Ｂb系とは、ステップＳ５で得られる局所座標系であり、Ｂs系とは、ステップＳ１２１で得られる局所座標系である。
【０１５７】
ここで、Ｂs系のＢb系に対する姿勢を表すオイラーパラメータをｑBbBsとおく。このとき、ｑBbBsは、以下のとおり表される。
【数４８】
0006776882.tif 000049
【０１５８】
姿勢誤差は、Ｂs系又はＢb系におけるシャフト軸と、Ｂ系におけるシャフト軸とをそれぞれ慣性系に変換したときの軸どうしの為す角度で評価される。なお、Ｂ系とは、局所座標系（真値）である。姿勢誤差をシャフト軸の向きの違いで評価するために、ｑBbBsを次のようにシャフト軸を一致させる回転(シャフト軸の向きの違い)と、シャフト軸周りの回転とに分解して考える。ｑsがシャフト軸を一致させる回転に対応し、ｑψがシャフト軸周りの回転に対応する。
【数４９】
0006776882.tif 000050
【０１５９】
ｓとは、−ｚ軸の向きを表す単位ベクトルであり、ｓ＝（０，０,−１）である。ｑsは、ｘｙ平面内の単位ベクトルα＝（αx，αy,０）の周りにＢb系のｓ軸を角度φだけ回転させる操作を表す。ｑψは、ｓ軸の周りに角度ψだけ回転させる操作を表す。従って、姿勢誤差は、角度φを用いて評価される。また、αとｓは直交するため、α・ｓ＝０となる。数４９の式を展開すると、以下のとおりとなる。
【数５０】
0006776882.tif 000051
【０１６０】
上式の右辺第１項〜第３項は、以下のとおり変形される。
【数５１】
0006776882.tif 000052
以上より、オイラーパラメータｑのベクトル部をＶ（ｑ）と表し、スカラー部をＳ（ｑ）と表すとき、以下のとおりとなる。
【数５２】
0006776882.tif 000053
【０１６１】
数５２の第１式を第２式で割ることにより、ψを表すことができる。
【数５３】
0006776882.tif 000054
【０１６２】
ステップＳ１２２では、導出部１１４ｃは、数４８の式にステップＳ１２１で導出されたインパクト時の姿勢ｑIBs、及びステップＳ５で導出されたインパクト時の姿勢ｑIBbを代入することにより、インパクト時のｑBbBsを導出する。続いて、導出部１１４ｃは、このｑBbBsを数５３の式に代入することにより、インパクト時のψを導出する。さらに、導出部１１４ｃは、数４９の第１式を変形した以下の式に、このｑBbBsとψとを代入することにより、インパクト時のｑsを導出する。
【数５４】
0006776882.tif 000055
【０１６３】
さらに、導出部１１４ｃは、ステップＳ５で導出されたインパクト時の姿勢ｑIBbと、インパクト時のｑsとを下式に代入することにより、インパクト時の姿勢ｑIBbを補正したｑIBb’を導出する。数５５の式による変換により、Ｂb系の姿勢誤差とＢs系の姿勢誤差とが等しくなる。
【数５５】
0006776882.tif 000056
【０１６４】
続くステップＳ１２３，Ｓ１２４では、導出部１１４ｃは、ダウンスイング期間の始端の時刻（すなわち、インパクトのＴ0前）と終端の時刻（すなわち、インパクトの時刻）との間の姿勢ｑIBbを補間する。このとき、インパクトの時刻における姿勢ｑIBbが、ステップＳ１２２で導出された補正後の姿勢ｑIBb’に一致するようにする。具体的には、まず、ｑIBbの微分を考える。
【数５６】
0006776882.tif 000057
ここで、ωBmは、角速度センサ４２の検出値である。ωbiasは、ωBmのバイアス成分である。ステップＳ１２３では、このバイアス成分ωbiasが導出される。なお、数５６では、バイアス成分が考慮されているという意味で、ｑIBbにハットの記号が付されている。
【０１６５】
具体的には、導出部１１４ｃは、ダウンスイング期間に含まれる各時刻において、ｑIBb及びωBmを数５６の右辺に代入する。このとき代入されるｑIBbは、ダウンスイング期間の始端の時刻においては、ステップＳ５で算出された値である。それ以降の時刻においては、直前の時刻において数５６の式から算出されたｑIBbの微分と、数６の式とに基づいて更新されるクォータニオンが代入される。そして、これらの全ての時刻の代入値を加算する。言い換えると、ダウンスイング期間を積分期間として、数５６の式が積分される。この加算値すなわち積分値は、ωbiasを未知数として含む式であり、インパクト時の姿勢ｑIBbを表している。従って、導出部１１４ｃは、この積分値（式）を、インパクト時の姿勢ｑIBbとして数４８の式に代入して、ｑBbBsの値（式）を導出する。さらに、このｑBbBsの値（式）を用いて、数５３,５４の式からｑsの値（式）を導出する。そして、このｑsの値（式）が、ステップＳ１２２で先に導出したｑsの値と等しくなるように、ωbiasを決定する。ここで、ｑsのｚ成分は０となるため、満たすべき終端条件は２成分である。ωbiasは３つの成分を有するので、３つの成分を２つの終端条件を満たすように計算する。このとき、ニュートン法等を用いて効率的に、ωbiasの３つの成分を探索することが好ましい。
【０１６６】
次に、Ｓ１２４では、ステップＳ１２３で導出されたバイアス成分ωbiasに基づいて、ダウンスイング期間における姿勢ｑIBbを導出し直す。具体的には、ステップＳ５で導出されたダウンスイング期間の始端の時刻における姿勢ｑIBbを、数６，７の式を用いて時々刻々更新してゆく。このとき、更新に用いられる角速度は、バイアス成分を加味した角速度（ωBm＋ωbias）とされる。
【０１６７】
なお、本実施形態では、ステップＳ１２４において、バイアス成分を加味した角速度（ωBm＋ωbias）を算出するに当たり、各時刻におけるバイアス成分ωbiasが以下の式で補正される。Ｔは、ダウンスイング期間の長さであり、ｔは、ダウンスイング期間の始端からの経過時間である。ωbiasfは、インパクト時のバイアス成分である。
【数５７】
0006776882.tif 000058
これにより、３段落前の説明で導出されたバイアス成分ωbiasは、インパクトの時刻のバイアス成分ωbiasfとされ、ダウンスイング期間の始端の時刻のバイアス成分は、０とされる。また、ダウンスイング期間中にバイアス成分が徐々に大きくなるように設定される。これにより、ダウンスイング期間とそれ以外の期間との境界（始期の時刻とその直前の時刻との境界）で角速度が急激に変化することを避け、角速度を連続的に変化させることができる。
【０１６８】
また、バイアス成分の補正方法としては、以下の方法も考えられ、以下の補正方法は、以上の補正方法と組み合わせて使用することもできるし、単独で使用することもできる。すなわち、バイアス成分は、高速域では大きくなり、低速域では小さくなる。従って、ωBmは下式のとおり補正することができる。
ωBm’＝ＴＨ＋（ωBm−ＴＨ）×ｋ      （ωBm＞ＴＨの場合）
ωBm’＝−ＴＨ＋（ωBm＋ＴＨ）×ｋ    （ωBm＜ＴＨの場合）
【０１６９】
ここで、ωBm’は、補正後の角速度（バイアス成分を加味した角速度）であり、ＴＨは閾値であり、ｋは、補正係数である。以上により、ステップＳ１０５の補正処理が終了する。
【０１７０】
＜２−１−４．グリップの位置及び速度の導出（Ｓ１０６）＞
グリップの位置及び速度の導出処理（ステップＳ１０６）は、ステップＳ６の処理と以下の点を除き同様である。すなわち、ステップＳ１０６では、正面画像データのみならず、側面画像データに基づいて、グリップ５１の位置及び速度が導出される。具体的には、数３９の式が以下のとおり変更される。ｂは、加速度センサのバイアス成分を意味している。
【数５８】
0006776882.tif 000059
【０１７１】
＜２−２．評価＞
以下、本発明の実施例について説明するが、本発明は、以下の実施例に限定されない。
【０１７２】
実施例５として、ステップＳ１０５の最適化の効果を確認すべく、第２実施形態に係るステップＳ１〜Ｓ１０５を実行して得られた最適化された姿勢ｑを用いて、局所座標系のシャフトの向き[０，０，１]Tを全体座標系でのシャフトの向きに変換した。そして、このシャフトの向きのベクトルと、モーションキャプチャーシステムで計測した全体座標系でのシャフトの向きを表すベクトル（真値）との為す角度を、角度誤差として算出した。この結果を、図２３に一点鎖線で示す。なお、図２３には、図１５の実施例２及び参考例２も示している。同図からは、ステップＳ１０５の補正処理を経ることにより、インパクト直前の期間においても、真値との誤差が小さくなっていることが分かる。
【０１７３】
＜３．変形例＞
以上、本発明の一実施形態について説明したが、本発明は上記実施形態に限定されるものではなく、その趣旨を逸脱しない限りにおいて、種々の変更が可能である。例えば、以下の変更が可能である。
【０１７４】
＜３−１＞
上記ステップＳ２６の最適化処理では、地磁気データ由来の関数Ｆmと、画像データ由来の関数Ｆkを含む目的関数Ｊ1が最小化された。しかしながら、目的関数Ｊ1は、関数Ｆm,Ｆkの一方を含まないように定義することもできる。同様に、上記ステップＳ３７の最適化処理では、Ｆn，Ｆm，Ｆkを含む目的関数Ｊ2が最小化されたが、目的関数Ｊ2は、関数Ｆn，Ｆm，Ｆkの少なくとも１つを含まないように定義することもできる。また、地磁気データ由来の関数Ｆmを、加速度データ由来の関数として定義することもできる。この場合、例えば、加速度センサ４１により計測されたａx，ａy，ａzから重力成分を抽出する。そして、地磁気ｍx，ｍy，ｍzのデータを、それぞれ抽出された重力の向きのデータに置き換えることにより、関数Ｆmを同様に算出することができる。
【０１７５】
＜３−２＞
上記実施形態では、距離画像センサ２によりＩＲ画像が撮影されたが、ＩＲ画像に代えてカラー画像を撮影し、これを運動の解析に用いることもできる。この場合、距離画像センサ２には、可視光を受光する可視光受光部（例えば、ＲＧＢカメラ）を搭載すればよい。
【０１７６】
＜３−３＞
上記実施形態に係る運動解析装置、方法及びプログラムは、ゴルフスイングの運動を解析するのに適するように構成されていたが、同様のアルゴリズムは、任意の運動を解析するのに用いることができる。
【０１７７】
＜３−４＞
上記ステップＳ２１では、初期時刻ｔ1の姿勢（オイラーパラメータｑ及び姿勢行列Ｎ）が、初期時刻ｔ1の画像データ、及び、初期時刻ｔ1の加速度データから特定される重力の向きに基づいて導出された。しかしながら、初期時刻ｔ1の姿勢は、角速度データから導出することもできる。図２４のフローチャートは、その一例の処理の流れを示している。
【０１７８】
まず、ステップＳ７１において、導出部１４ｃは、初期時刻ｔ1のオイラーパラメータｑを［０，０，０，１］とし、当該オイラーパラメータｑを、初期時刻ｔ1からインパクト付近までの角速度ωx，ωy，ωzのデータを用いて、数６，７の式に従って時々刻々更新してゆく。その後、数５の式に従って、各時刻の姿勢ｑを、姿勢行列Ｎに換算する。続いて、導出部１４ｃは、ステップＳ７２において、ステップＳ２３と同様に、時刻ｔi，ｔjにおける姿勢の差ｑBjBiを導出する。また、ステップＳ７３において、ステップＳ２４と同様に、行列Ａkjを導出する。その後、ステップＳ７４において、ステップＳ２６と同様に、ラグランジュの未定乗数法により、ステップＳ７１で導出された姿勢ｑiが最適化される。具体的には、以下の評価関数Ｊ3を最小化するように最適化が行われる。
【数５９】
0006776882.tif 000060
【０１７９】
この最適化の方法も、ステップＳ２６と同様であり、評価関数Ｊ3が最小化されるためには、目的関数Ｊ3の微分に関し、数６０が成立する必要がある。そして、この式は、数６１の行列の固有値がλとなり、ｑiが固有ベクトルになることを表しているから、導出部１４ｃは、数６１の行列を計算し、当該行列の最小の固有値に対応する大きさ１の固有ベクトルを導出して、これをｑiとする。導出部１４ｃは、このｑiを少なくとも初期時刻ｔ1について導出し、さらに姿勢行列Ｎに変換する。これにより、初期時刻ｔ1における姿勢が導出される。
【数６０】
0006776882.tif 000061
【数６１】
0006776882.tif 000062
【０１８０】
＜３−５＞
第１及び第２実施形態におけるアドレスからインパクト付近までの姿勢の導出処理（ステップＳ５）において、ステップＳ３２〜Ｓ３７を省略することもできる。この場合、ステップＳ６以降の処理においては、ステップＳ３７の結果に代えて、ステップＳ３１で導出された仮の姿勢を用いることができる。
【符号の説明】
【０１８１】
１，１０１  運動解析装置
２  距離画像センサ
６  距離画像センサ
３  運動解析プログラム
４  慣性センサユニット
４１  加速度センサ
４２  角速度センサ
４３  地磁気センサ
５  ゴルフクラブ
７  ゴルファー
１４ａ  第１取得部
１４ｂ  第２取得部
１４ｅ  第３取得部
１４ｃ，１１４ｃ  導出部
５１  グリップ
５２  シャフト
１００，２００  運動解析システム
【図１】
0006776882.tif 000063
【図２】
0006776882.tif 000064
【図３】
0006776882.tif 000065
【図４】
0006776882.tif 000066
【図５】
0006776882.tif 000067
【図６】
0006776882.tif 000068
【図７】
0006776882.tif 000069
【図８】
0006776882.tif 000070
【図９】
0006776882.tif 000071
【図１０】
0006776882.tif 000072
【図１１】
0006776882.tif 000073
【図１２】
0006776882.tif 000074
【図１３】
0006776882.tif 000075
【図１４】
0006776882.tif 000076
【図１５】
0006776882.tif 000077
【図１６Ａ】
0006776882.tif 000078
【図１６Ｂ】
0006776882.tif 000079
【図１６Ｃ】
0006776882.tif 000080
【図１７Ａ】
0006776882.tif 000081
【図１７Ｂ】
0006776882.tif 000082
【図１７Ｃ】
0006776882.tif 000083
【図１８】
0006776882.tif 000084
【図１９】
0006776882.tif 000085
【図２０】
0006776882.tif 000086
【図２１】
0006776882.tif 000087
【図２２】
0006776882.tif 000088
【図２３】
0006776882.tif 000089
【図２４】
0006776882.tif 000090