JP6954039B2 - 波形推定装置、波形推定方法および波形推定プログラム - Google Patents

Description

本発明は、波形推定装置、波形推定方法および波形推定プログラムに関する。

従来から、運転への集中力低下で発生する、走行中の危険なふらつき挙動を検出してフィードバックすることで、安全運転を支援する装置やシステムが知られている。例えば、操舵角や横変位の時間変化による波形を高周波成分と低周波成分に分類し、低周波成分を高周波成分で除算したパワー比が高いときに低覚醒状態と判定する技術が知られている。

特開平７−１２５５６０号公報 特開平６−１９７８９１号公報 特開平３−２０３４１２号公報

しかしながら、上記の技術では、操舵行動に関連した波形成分をリアルタイムに推定することが難しく、結果として、リアルタイムに安全運転を支援することが難しい。例えば、上記技術では、覚醒度を精度良く推定するために、５０秒から８０秒程度の長い時間、変位量データの蓄積が必要となり、波形分類のための信号処理に時間がかかってしまう。

１つの側面では、操舵行動に関連した波形成分をリアルタイムに推定することができる波形推定装置、波形推定方法および波形推定プログラムを提供することを目的とする。

第１の案では、波形推定装置は、車道を走行する車両の走行軌跡に応じた入力波形と、前記入力波形から算出される算出済みの第１周波数未満の第１振動成分とを用いて、現時刻における第１振動成分を推定する第１推定部を有する。波形推定装置は、前記入力波形と、算出済みの前記第１振動成分と、前記入力波形から算出される算出済みの第１周波数以上の第２振動成分から予測される前記現時刻における第２振動成分とを用いて、前記現時刻における第１振動成分を推定する第２推定部を有する。波形推定装置は、推定された各前記第１振動成分に基づいて、前記入力波形から前記現時刻における第２振動成分を算出する算出部を有する。

一実施形態によれば、操舵行動に関連した波形成分をリアルタイムに推定することができる。

図１は、実施例１にかかる走行制御システムの機能構成を示す機能ブロック図である。 図２は、車両の走行軌跡を説明する図である。 図３は、安全運転のリアルタイム判定を説明する図である。 図４は、実施例１にかかる波形推定を説明する図である。 図５は、実施例１にかかる車載装置の機能構成を示す機能ブロック図である。 図６は、走行軌跡の各成分波形を説明する図である。 図７は、振動成分の抽出時刻を説明する図である。 図８は、第１波形推定部の機能構成を示す機能ブロック図である。 図９は、方式Ａによる第１平均波形の推定を説明する図である。 図１０は、第１平均波形の推定を説明する図である。 図１１は、方式Ｂによる第１平均波形の推定を説明する図である。 図１２は、第１揺動波形の算出を説明する図である。 図１３は、方式Ａの入力波形適合の特性を説明する図である。 図１４は、方式Ｂの予測適合の特性を説明する図である。 図１５は、振幅拡大時のハイブリッド波形推定の特性を説明する図である。 図１６は、振幅縮小時のハイブリッド波形推定の特性を説明する図である。 図１７は、方式Ａのシミュレーション結果（振幅拡大時）を説明する図である。 図１８は、方式Ｂのシミュレーション結果（振幅拡大時）を説明する図である。 図１９は、ハイブリッド波形推定のシミュレーション結果（振幅拡大時）を説明する図である。 図２０は、方式Ａのシミュレーション結果（振幅縮小時）を説明する図である。 図２１は、方式Ｂのシミュレーション結果（振幅縮小時）を説明する図である。 図２２は、ハイブリッド波形推定のシミュレーション結果（振幅縮小時）を説明する図である。 図２３は、ハイブリッド波形推定の波形パワー推定のシミュレーション結果を説明する図である。 図２４は、各方式の波形パワー推定のシミュレーション結果を説明する図である。 図２５は、危険判定部の機能構成を示す機能ブロック図である。 図２６は、処理の流れを示すフローチャートである。 図２７は、波形推定処理の流れを示すフローチャートである。 図２８は、危険判定処理の流れを示すフローチャートである。 図２９は、操舵評価部の処理を説明する図である。 図３０は、急操舵の評価に使用するフィルタを説明する図である。 図３１は、急操舵の判定を説明する図である。 図３２は、揺動評価部の処理を説明する図である。 図３３は、全体揺動の評価を説明する図である。 図３４は、ハードウェア構成例を説明する図である。

以下に、本願の開示する波形推定装置、波形推定方法および波形推定プログラムの実施例を図面に基づいて詳細に説明する。なお、この実施例によりこの発明が限定されるものではない。各実施例は、矛盾のない範囲内で適宜組み合わせることができる。

［全体構成］
図１は、実施例１にかかる走行制御システムの機能構成を示す機能ブロック図である。図１に示すように、走行制御システムは、ＧＰＳセンサ１と、車載装置１０と、自動運転制御装置２と、アクチュエータ３とを有する。走行制御システムは例えば車両に搭載される。

ＧＰＳ（Global Positioning System）センサ１は、車両に搭載されるセンサであり、車両の位置情報を取得する。ＧＰＳセンサ１は、取得した位置情報を車載装置１０および自動運転制御装置２に出力する。

車載装置１０は、車両に搭載されるコンピュータであり、ＧＰＳセンサ１から位置情報を取得し、車両の走行軌跡などを生成する。また、車載装置１０は、例えば白線検出処理などの公知の手法により、車両の走行軌跡を生成することもできる。また、車載装置１０は、走行軌跡を生成するために使用される車両の運行情報などを様々な公知の手法やＧＰＳセンサ１等により取得することもできる。

また、車載装置１０は、車両の走行軌跡や運行情報に基づき、運転者の運転を評価する。例えば、車載装置１０は、走行軌跡の波形を、無意識の舵取りの大きさが表れる長周期、意識的な舵取りの大きさが表れる中周期、突発的な舵取りの大きさが表れる短周期の３種の成分波形に分離する。このとき、車載装置１０は、各周期の半波長前から現時刻までの各成分波形を推定することで、走行軌跡の波形を現時刻までの各成分波形に分離することができる。そのため、車載装置１０は、波形を各成分波形に分離する処理（信号処理）にかかる時間を短縮することができる。

続いて、車載装置１０は、各波形の特徴量から、全体揺動および急操舵を評価する。車載装置１０は、走行軌跡の波形を現時刻までの各成分波形に分離することができるため、各評価をリアルタイムで行うことができる。そして、車載装置１０は、リアルタイムでの評価結果からふらつきの危険度を算出し、閾値を超えた場合に危険ふらつきと判定するので、危険な車両のふらつきを検出するまでの時間を短縮することができる。

ここで、走行軌跡と検出ポイントとの関係について説明する。図２は、車両の走行軌跡を説明する図である。図２に示すように、車両は、白線の破線と白線の実線との間の斜線を走行する。車両が破線または実線を外れて走行すると、他車線を走行する他車両などとの衝突などの危険が発生する。

例えば、図２に示す走行軌跡ａにおいて、危険ふらつきとは、全体揺動（図２のｂ）が高まっており、かつ、急操舵軌跡（図２のｃ）が道路端（図２のｄ）かつ揺動端で発生している状態である。ここで、全体揺動の高まりは車線維持集中力の低下、急操舵軌跡は反応遅れの突発的対応、道路端での急操舵は逸脱を回避する操舵タイミングを示し、車載装置１０はこれを評価することで危険度を判定する。車載装置１０は、判定結果を自動運転制御装置２に出力する。

自動運転制御装置２は、車両に搭載されるコンピュータであり、車両の自動運転を制御する。自動運転制御装置２は、例えば車載装置１０が危険ふらつきと判定した場合に、車両の運転モードを運転者が手動運転を行う手動運転モードから、車両を自動的に走行させる自動運転モードに切り替える。また、自動運転制御装置２は、例えば運転者からの指示にしたがって運転モードを切り替えることもできる。自動運転制御装置２は、自動運転モードに切り替えると、ＧＰＳセンサ１からの位置情報等に基づき、アクチュエータ３を制御して舵角や車速を制御することで、車両が自動的に走行するように自動運転制御を行う。また、自動運転制御装置２は、例えば白線検出処理などの公知の手法により、車両の自動運転制御を行うこともできる。

アクチュエータ３は、例えばアクセルアクチュエータ、ブレーキアクチュエータおよび操舵アクチュエータ等を含む。アクセルアクチュエータは、エンジンのスロットルバルブの開度を調整する。アクセルアクチュエータは、自動運転制御装置２からの指示に従い、スロットルバルブの開度を調整する。ブレーキアクチュエータは、各車輪のホイールシリンダそれぞれに油圧を加えることで、各車輪と共に回転するディスクやドラム等の回転を抑制する。ブレーキアクチュエータは、自動運転制御装置２からの指示に従い、油圧を調整する。操舵アクチュエータは、ステアリング機構に操舵トルクを付与する。操舵アクチュエータは、自動運転制御装置２からの指示に従い、操舵トルクを調整する。

［リアルタイム判定］
続いて、車載装置１０が実行する安全運転のリアルタイム判定の手法について説明する。図３は、安全運転のリアルタイム判定を説明する図である。安全運転の判定としては、一般的に、走行軌跡を蓄積した上で安全運転を判定する蓄積手法が知られている。具体的には、蓄積手法は、走行軌跡から操舵波形（揺動成分）を分解する際に、高周波（微舵振動）から低周波（中心揺動）に順に濾波する。そして、入力波形から波長推定、平均波長推定、波形分解、パワー算出順に、全フレームを走査しながら順次算出する。

例えば、入力波形から平均波形Ａと揺動波形Ａを抽出し、揺動波形Ａを用いて、細かなハンドル修正やノイズなどを示す微舵操舵成分を抽出する。さらに、入力波形から抽出された平均波形Ａから、平均波形Ｂと揺動波形Ｂを抽出し、揺動波形Ｂを用いて、舵取りのための主たる操作などを示す修正操舵成分を抽出する。さらに、平均波形Ｂから、平均波形Ｃと揺動波形Ｃを抽出し、揺動波形Ｃを用いて、車線全体をふらつきながら走行などを示す路内揺動成分を抽出する。さらに、平均波形Ｃを用いて、路内振動の中心の変化を示す中心振動成分を抽出する。このようにして、入力波形から安全運転が判定される。なお、波形分解とは、平均波形と揺動波形へ分解することをいい、平均波形は、より低周波の波形分解の入力として使用され、揺動波形は、着目周波数帯でのパワー算出に使用される。

上述したように、揺動波形の推定を行うためには、着目時間の前後１／２波長が観測されている必要があり、着目時間より未来の１／２波長のデータが必要となるため、蓄積データのように、取りためられたデータでのみ実施することができる。したがって、走行中のリアルタイム分析ができない。具体的には、図３に示すように、現在（最新時刻）で観測された入力波形からでは領域（ｚ）に該当する波形の情報が不足する。すなわち、波長推定と平均波長推定と波形分解を行うために用いる１／２波長前の波形が不足するとともに、パワー推定に用いる１波長前の波形が不足する。したがって、波長に応じた時間遅れが発生し、危険ふらつきの判定に遅れが発生する。つまり、安全運転のリアルタイム判定には、蓄積を前提としない時系列の逐次フレーム処理での波形分解と、波形分解のフィルタリング遅れの補償とが求められる。

［実施例１の波形推定］
次に、実施例１にかかる車載装置１０が実行する波形推定について説明する。上述したように、入力波形から抽出される振動成分（平均波形）とそれ以外の振動成分（揺動波形）のうち、揺動波形を用いて波形パワーが算出される。ただし、リアルタイムに安全判定を行うためには、現在時刻より１／２波長前の波形が不足しているので、この現在時刻より１／２波長前（直近１／２波長前）の波形を推定する必要がある。そこで、実施例１では、直近１／２波長前の波形推定を実行する。

図４は、実施例１にかかる波形推定を説明する図である。まず、自己回帰モデルや自己回帰移動平均モデルなどの一般的な手法について説明する。図４の（ａ）に示すように、一般的な波形推定手法は、時系列で過去の波形形状をモデル化し、将来の波形を推定することで、揺動波形を直接推定する。

しかし、一般的な手法では、半周期前より以前の傾向しか取り込めないので、半周期前までの直近の傾向を反映できない。このため、直近半波長（直近１／２波長）の傾向が、推定に用いたそれ以前の波形と異なると、１／２波長分の推定形状と現時刻付近の波形形状とが大きくかい離してしまう。また、自己回帰移動平均モデルを用いて、入力が存在しない１／２波長内を推定した値は、定数に収束する挙動を示し、揺動を十分に表せない。このように、一般的な手法では、現時刻付近の波形を正しく推定できない。

これに対して、実施例１にかかる車載装置１０は、波形分解の方法に着目し、平均波形を推定する。具体的には、揺動波形を直接推定せずに、平均波形を推定して入力波形から差し引くことで、揺動波形を算出する。そこで、実施例１にかかる車載装置１０は、図４の（ｂ）に示すように、方式Ａ（入力波形適合）と方式Ｂ（予測適合）とを組み合わせたハイブリッド波形推定方式を採用する。

方式Ａ（入力波形適合）は、現在時刻の１／２波長前から現在時刻までの算出済みの入力波形と、現在時刻の１波長前から１／２波長前までの算出済みの平均波形とから、直近１／２波長前までの平均波形を推定する。例えば、方式Ａは、１次式で平均波形（ｙ１＝ａ×ｔ＋ｂ）を推定する。

方式Ｂ（予測適合）は、現在時刻の１／２波長前から現在時刻までの算出済みの入力波形と、現在時刻の１波長前から１／２波長前までの算出済みの平均波形と、現在時刻の１波長前から１／２波長前までの算出済みの揺動波形から予測された直近１／２波長前の揺動波形とから、直近１／２波長前までの平均波形を推定する。例えば、方式Ｂは、２次式で平均波形（ｙ２＝ｃ×ｔ^２＋ｄ×ｔ＋ｅ）を推定する。

そして、ハイブリッド波形推定方式は、方式Ａで推定された平均波形と方式Ｂで推定された平均波形とを加算して、現在時刻の１／２波長前から現在時刻までの平均波形を推定する。例えば、ハイブリッド波形推定方式は、１次式の平均波形（ｙ１＝ａ×ｔ＋ｂ）と２次式の平均波形（ｙ２＝ｃ×ｔ^２＋ｄ×ｔ＋ｅ）との和「ｙ＝０．５×ｃ×ｔ^２＋０．５×（ａ＋ｄ）×ｔ＋０．５×（ｂ＋ｅ）」を算出する。

その後、車載装置１０は、推定された平均波形を入力波形から減算することで、現在時刻の１／２波長前から現在時刻までの揺動波形を推定する。そして、車載装置１０は、推定された揺動波形から波形パワーを算出して、安全運転のリアルタイム判定を実行する。このように、車載装置１０は、揺動の増減変化で逆傾向の特性を持つ方式Ａと、定常時特性に優れる方式Ｂとを組み合わせて、揺動増減の追従性に優れた揺動波形の推定を実現する。

［機能構成］
図５は、実施例１にかかる車載装置１０の機能構成を示す機能ブロック図である。なお、自動運転制御装置２は、車両の自動運転を行うコンピュータであり、一般的な自動運転制御装置と同様の機能を有するので、ここでは詳細な説明は省略する。

図５に示すように、車載装置１０は、記憶部１１および制御部２０を有する。記憶部１１は、プログラムやデータなどを記憶する記憶装置であり、例えばメモリやハードディスクなどである。この記憶部１１は、軌跡ＤＢ１２と、波形ＤＢ１３と、揺動波形ＤＢ１４、評価値ＤＢ１５とを記憶する。軌跡ＤＢ１２は、制御部２０が生成した走行軌跡や運行情報などを走行した時間と対応付けて記憶するデータベースである。

波形ＤＢ１３は、走行軌跡の波形を長周期、中周期および短周期に分離した各成分波形（振動成分）を記憶する。波形ＤＢ１３は、制御部２０が生成した各成分波形を走行した時間と対応付けて記憶するデータベースである。

ここで、走行軌跡の波形を分離した各成分波形について説明する。図６は、走行軌跡の各成分波形を説明する図である。図６に示すように、走行軌跡の波形ａは、短周期の急操舵に関する振動成分ａ１、中周期の通常操舵に関する振動成分ａ２および長周期の路内揺動に関する振動成分ａ３に分離される。各振動成分ａ１〜ａ３は、各時刻で周期が異なるが、その周期はそれぞれ定めた一定範囲内となる。

例えば、短周期には１〜４秒、中周期には８〜１６秒、長周期には６０〜９０秒の周期波形が含まれる。周期範囲は一例で有り、周期範囲が重複しなければこれに限らずともよい。振動成分ａ１〜ａ３における各時刻の周期長は、例えば走行軌跡の波形を、抽出時刻を中心としてＤＦＴ（離散フーリエ変換）に代表される周波数領域への変換を行い、対象とする周波数範囲を代表する応答成分で決定することができる。代表応答周期（周波数）の決定は、周波数範囲で最大応答を示した周波数（周期）を用いるか、または、周波数範囲における平均値、中間値などの統計量を用いて決定してもよい。

揺動波形ＤＢ１４は、過去に算出された揺動波形を記憶するデータベースである。具体的には、揺動波形ＤＢ１４は、測定された走行軌跡から得られる各周期の成分波形に対応する揺動波形を記憶する。

評価値ＤＢ１５は、危険度の評価結果を記憶するデータベースである。例えば、評価値ＤＢ１５は、評価値を走行した時間と対応付けて記憶する。ここで、評価値は、走行軌跡に応じて評価された値であり、数値が高いほど危険度が高い運転を示す。

制御部２０は、車載装置１０全体を司る処理部であり、例えばプロセッサなどである。この制御部２０は、軌跡生成部２１、第１波形推定部２２、第２波形推定部２３、第３波形推定部２４および危険判定部２５を有する。なお、軌跡生成部２１、第１波形推定部２２、第２波形推定部２３、第３波形推定部２４および危険判定部２５は、プロセッサなどの電子回路の一例やプロセッサが実行するプロセスの一例である。

軌跡生成部２１は、ＧＰＳセンサ１から取得した位置情報にしたがって、車両の走行軌跡を生成する処理部である。また、軌跡生成部２１は、ドライブレコーダ等で撮像された走行時の画像情報等を受信し、白線検出処理などを用いて車両の走行軌跡を生成することもできる。軌跡生成部２１は、生成した走行軌跡を軌跡ＤＢ１２に格納する。

第１〜第３波形推定部２２〜２４は、走行を開始してから現時刻までの短周期、中周期および長周期の成分波形（振動成分ａ１〜ａ３）をそれぞれ走行軌跡の波形から分離する処理部である。第１波形推定部２２は、走行軌跡の波形を入力波形として、短周期の振動成分ａ１（第１揺動波形）とそれ以外の振動成分（第１平均波形）とに分離する。また、第１波形推定部２２は、短周期の振動成分ａ１の波形パワーを算出する。なお、第１波形推定部２２の処理の詳細は後述する。

続いて、第２波形推定部２３は、第１波形推定部２２が分離した第１平均波形を入力波形として、中周期の振動成分ａ２（第２揺動波形）とそれ以外の振動成分（第２平均波形）とに分離する。また、第２波形推定部２３は、中周期の振動成分ａ２の波形パワーを算出する。第３波形推定部２４は、第２波形推定部２３が分離した第２平均波形を入力波形として低周期の振動成分ａ３（第３揺動波形）に分離する。また、第３波形推定部２４は、低周期の振動成分ａ３の波形パワーを算出する。このように、第１〜第３波形推定部２２〜２４は、高周波から順に振動成分ａ１〜ａ３を分離する。

ここで、各振動成分ａ１〜ａ３の抽出時刻について説明する。図７は、振動成分ａ１〜ａ３の抽出時刻を説明する図である。なお、図７では、各振動成分ａ１〜ａ３の１波長を模式的に示している。第１〜第３波形推定部２２〜２４は、抽出時刻を中心として１波長の入力波形を周波数領域へ変換して振動成分ａ１〜ａ３を分離する。そのため、第１〜第３波形推定部２２〜２４は、現時刻から半波長前までの振動成分ａ１〜ａ３を走行軌跡の波形から直接分離することができるが、半波長前から現時刻までの振動成分ａ１〜ａ３を直接分離することができない。そのため、例えば１．５秒の短周期の場合、現時刻から０．２５秒前までの間、振動成分ａ１〜ａ３が分離されず、例えば１６秒の中周期の場合、現時刻から８秒前までの間、振動成分ａ１〜ａ３が分離されない。また、例えば９０秒の長周期の場合、現時刻から４５秒前までの間、振動成分ａ１〜ａ３が分離されない。

このように、第１〜第３波形推定部２２〜２４が走行軌跡の波形から直接振動成分ａ１〜ａ３を分離する場合、半波長に相当する時間の経過を待ってからしか振動成分ａ１〜ａ３を分離することができない。したがって、第１〜第３波形推定部２２〜２４が直接分離した振動成分ａ１〜ａ３を用いて、後述する危険判定部２５が危険ふらつき判定を行うとすると、現時刻から９０秒前までの危険ふらつきしか判定を行うことができず、リアルタイムで判定を行うことができない。そこで、本実施例では、第１〜第３波形推定部２２〜２４が現時刻から半波長前までの振動成分ａ１〜ａ３を推定することで、リアルタイムで振動成分の分離を行えるようにした。これにより、危険判定部２５がリアルタイムで危険ふらつきを判定することができる。

（波形推定部の機能ブロック）
図８は、第１波形推定部２２の機能構成を示す機能ブロック図である。なお、第２、第３波形推定部２３、２４は、入力波形および分離する振動成分の周期が異なる点を除き、第１波形推定部２２と同様の処理を行うため、説明を省略する。

図８に示すように、第１波形推定部２２は、過去波形算出部２２１、揺動波形予測部２２２、平均波形推定部２２３、平均波形算出部２２４、揺動波形算出部２２５、波形パワー算出部２２６を有する。なお、過去波形算出部２２１は、振動算出部の一例であり、揺動波形予測部２２２は、予測部の一例であり、平均波形推定部２２３および平均波形算出部２２４は、第１推定部と第２推定部の一例である。揺動波形算出部２２５は、算出部の一例である。

過去波形算出部２２１は、入力波形（走行軌跡の波形）から、現時刻の半波長前までの第１揺動波形（振動成分）と第１平均波形とを分離する処理部である。具体的には、過去波形算出部２２１は、例えば全ｎ種の波長候補Ｌｃ＝｛Ｌ１、Ｌ２、・・・、Ｌｎ｝に対して、現時刻から半波長前の時刻を中心とする入力波形についてフーリエ変換を行う。そして、過去波形算出部２２１は、フーリエ変換係数のｓｉｎ値およびｃｏｓ値のベクトルサイズが最大となる波長を、現時刻から半波長前の時刻の波長に決定する。

続いて、過去波形算出部２２１は、決定した波長範囲で、入力波形を平均し、第１平均波形を算出する。過去波形算出部２２１は、現在時刻から半波長前までの入力波形から、算出した第１平均波形を減算して第１揺動波形（振動成分）を算出する。また、過去波形算出部２２１は、現時刻の１波長前の時刻を中心とする第１揺動波形について、当該時刻で決定した波長範囲でＲＭＳ（Root Mean Square）値を波形パワーとして算出する。過去波形算出部２２１は、決定した波長と、分離した第１揺動波形（振動成分）および第１平均波形と、波形パワーとを波形ＤＢ１３に格納する。また、過去波形算出部２２１は、決定した波長と、分離した第１揺動波形（振動成分）を揺動波形ＤＢ１４に格納する。

揺動波形予測部２２２は、過去波形算出部２２１によって算出済みの過去の第１揺動波形を用いて、現時刻での第１揺動波形を予測する処理部である。具体的には、揺動波形予測部２２２は、現時刻の１波長前から１／２波長前までの期間Ｔ１の第１揺動波形を揺動波形ＤＢ１４から取得し、波形の振幅や位相を特定する。そして、揺動波形予測部２２２は、特定した振幅や位相にしたがって、現時刻の１／２波長前から現時刻までの期間Ｔ２（直近１／２波長前）の第１揺動波形を予測する。その後、揺動波形予測部２２２は、予測された期間Ｔ２の第１揺動波形を平均波形推定部２２３等に出力する。

平均波形推定部２２３は、方式Ａ（入力波形適合）と方式Ｂ（予測適合）のそれぞれについて、期間Ｔ２の第１平均波形を推定する処理部である。具体的には、平均波形推定部２２３は、各方式における期間Ｔ２の第１平均波形を表す多項式のパラメータを算出する。

（方式Ａ：入力波形適合の説明）
図９は、方式Ａによる第１平均波形の推定を説明する図である。平均波形推定部２２３は、図９に示すように、期間Ｔ２の入力波形に加え、期間Ｔ１の第１平均波形を用いて期間Ｔ２の第１平均波形を推定する。これにより、第１平均波形をより高精度に推定することができる。

ここで、平均波形推定部２２３が第１平均波形のパラメータを推定する理由について説明する。図１０は、第１平均波形の推定を説明する図である。図１０に示す第１揺動波形（振動成分）ａ１は、第１平均波形ａ４より約４倍周期が短い。そのため、現時刻から短周期の半波長前の時刻までの間（期間Ｔ２）の第１平均波形ａ４は８分の１波長程度であり、振動成分ａ１は半波長程度である。したがって、期間Ｔ２の第１平均波形ａ４は、単調増加または単調減少する多項式でモデル化することができ、式（１）に示すように、低次の多項式でモデル化することができる。

なお、式（１）において、時刻ｔでの第１平均波形の位置をｙ１_ｔとし、ａとｂはパラメータである。また、現時刻から半波長前の時刻を基準（ｔ＝０）とする。したがって、現時刻は、短周期Ｔｓの２分の１の時刻（ｔ＝Ｔｓ／２）となる。パラメータｂは、ｔ＝０のときの第１平均波形ｙ１_０であり、１／２波長前の時刻ｔ＝０での実際に平均波形位置をｍ_０とするとき、式（２）に基づいて算出される。なお、ｔ＝０のときの第１平均波形ｙ１_０は、過去波形算出部２２１が算出した第１平均波形ｍ_ｔから得られる（ｙ１_０＝ｍ_０）。

平均波形推定部２２３は、期間Ｔ２において、第１平均波形ｙ１_ｔと入力波形ｐ_ｔとの差が最小化するように式（１）の各パラメータａを算出する。また、平均波形推定部２２３は、基準時刻ｔ＝０から半波長前までの期間Ｔ１において、第１平均波形ｙ１_ｔと過去波形算出部２２１が算出した第１平均波形ｍ_ｔとの差が最小化するようにパラメータａを算出する。具体的に、平均波形推定部２２３は、式（３）に示すように、第１平均波形ｙ１_ｔと波形ｑ_ｔとの差を算出する。なお、式（４）に示すように、波形ｑ_ｔは、ｔ≧０で入力波形ｐ_ｔ、ｔ＜０で過去波形算出部２２１が算出した第１平均波形ｍ_ｔとなる波形である。

続いて、平均波形推定部２２３は、残差ｓ_ｔの二乗残差和Ｓを式（５）に基づいて算出し、二乗残差和Ｓが最小となるように式（６）を解くことで、式（７）に示すようにパラメータａを算出する。ただし、積分範囲は、波長時間をＬ（秒）とするとき、−Ｌ／２からＬ／２までとする。

上述したように、平均波形推定部２２３は、期間Ｔ２の入力波形ｐ_ｔに加え、期間Ｔ１の第１平均波形ｍ_ｔを用いて期間Ｔ２の第１平均波形ｙ１_ｔを推定する。このように、平均波形推定部２２３は、期間Ｔ２における第１平均波形と期間Ｔ２の入力波形および期間Ｔ１の第１平均波形との差が最小化する多項式モデルの係数を算出する。

（方式Ｂ：予測適合の説明）
図１１は、方式Ｂによる第１平均波形の推定を説明する図である。平均波形推定部２２３は、図１１に示すように、観測済みである期間Ｔ２の入力波形に加え、算出済みである期間Ｔ１の第１平均波形と、算出済みである期間Ｔ１の揺動波形から予測される期間Ｔ２の第１揺動波形とを用いて、期間Ｔ２の第１平均波形を推定する。これにより、第１平均波形をより高精度に推定することができる。

方式Ｂについても、方式Ａと同様、現時刻から短周期の半波長前の時刻までの間（期間Ｔ２）の第１平均波形ａ４は８分の１波長程度であり、振動成分ａ１は半波長程度である。したがって、期間Ｔ２の第１平均波形ａ４は、単調増加または単調減少する多項式でモデル化することができ、式（８）に示すように、低次の多項式でモデル化することができる。

なお、式（８）において、時刻ｔでの第１平均波形の位置をｙ２_ｔとし、ｃ、ｄ、ｅはパラメータである。また、現時刻から半波長前の時刻を基準（ｔ＝０）とする。したがって、現時刻は、短周期Ｔｓの２分の１の時刻（ｔ＝Ｔｓ／２）となる。パラメータｅは、ｔ＝０のときの第１平均波形ｙ２_０であり、１／２波長前の時刻ｔ＝０での実際に平均波形位置をｍ_０とするとき、式（９）に基づいて算出される。なお、ｔ＝０のときの第１平均波形ｙ２_０は、過去波形算出部２２１が算出した第１平均波形ｍ_ｔから得られる（ｙ２_０＝ｍ_０）。

ここで、基準時刻ｔ＝０での入力波形の振幅をＡ、位相をＢとすると、揺動波形予測部２２２によって、時刻ｔ（ｔ≧０）における揺動波形は式（１０）にように予測される。

ここで、方式Ａと同様、平均波形推定部２２３は、期間Ｔ２において、第１平均波形ｙ２_ｔと入力波形ｐ_ｔとの差が最小化するように式（８）の各パラメータｃ、ｄを算出する。また、平均波形推定部２２３は、基準時刻ｔ＝０から半波長前までの期間Ｔ１において、第１平均波形ｙ２_ｔと過去波形算出部２２１が算出した第１平均波形ｍ_ｔとの差が最小化するようにパラメータｃ、ｄを算出する。具体的に、平均波形推定部２２３は、式（１１）に示すように、第１平均波形ｙ２_ｔと波形ｑ_ｔとの差を算出する。なお、式（１２）に示すように、波形ｑ_ｔは、ｔ≧０で入力波形ｐ_ｔ−ｕ_ｔ、ｔ＜０で過去波形算出部２２１が算出した第１平均波形ｍ_ｔとなる波形である。

続いて、平均波形推定部２２３は、残差ｓ_ｔの二乗残差和Ｓを式（１３）に基づいて算出し、二乗残差和Ｓが最小となるように式（１４）、式（１５）を解くことで、式（１６）に示すようにパラメータｃ、ｄを算出する。ただし、積分範囲は、波長時間をＬ（秒）とするとき、−Ｌ／２からＬ／２までとする。

このようにして、平均波形推定部２２３は、期間Ｔ２の入力波形ｐ_ｔに加え、期間Ｔ２の予測された揺動波形ｕ_ｔ、期間Ｔ１の第１平均波形ｍ_ｔを用いて、期間Ｔ２の第１平均波形ｙ２_ｔを推定する。そして、平均波形推定部２２３は、算出した上記各パラメータを平均波形算出部２２４に出力する。

平均波形算出部２２４は、平均波形推定部２２３が算出したパラメータａ〜ｅに基づき、期間Ｔ２の第１平均波形ｙ_ｔを算出する処理部である。具体的には、平均波形算出部２２４は、方式Ａにより推定された第１平均波形ｙ１_ｔと、方式Ｂにより推定された第１平均波形ｙ２_ｔとを用いて、期間Ｔ２の第１平均波形ｙ_ｔを算出する。すなわち、平均波形算出部２２４は、上述したハイブリッド波形推定を実行する。

例えば、平均波形算出部２２４は、算出されたパラメータａとｂとが設定された式（１）と、算出されたパラメータｃとｄとｅとが設定された式（２）とを合算して、式（１７）を生成する。そして、平均波形算出部２２４は、式（１７）によって算出される、期間Ｔ２の第１平均波形ｙ_ｔを揺動波形算出部２２５に出力する。

揺動波形算出部２２５は、平均波形推定部２２３から取得した第１平均波形ｙ_ｔに基づき、短周期Ｔｓの振動成分である第１揺動波形を算出する処理部である。図１２は、第１揺動波形の算出を説明する図である。図１２に示すように、揺動波形算出部２２５は、平均波形算出部２２４によって算出された期間Ｔ２の第１平均波形ｙ_ｔを、観測された期間Ｔ２の入力波形ｐ_ｔから減算することで、期間Ｔ２における第１揺動波形ｒ_ｔを算出する。例えば、揺動波形算出部２２５は、式（１８）に基づき、期間Ｔ２における第１揺動波形ｒ_ｔを算出する。

波形パワー算出部２２６は、振動成分である第１揺動波形ｒ_ｔから波形パワーを算出する処理部である。具体的には、図１２に示すように、波形パワー算出部２２６は、期間Ｔｓ（Ｔ１＋Ｔ２）の第１揺動波形に基づいて波形パワーを算出する。すなわち、波形パワー算出部２２６は、算出済みの期間Ｔ１の第１揺動波形と、揺動波形算出部２２５によって予測された第１揺動波形とを用いて、波形パワーを算出する。

より詳細には、波形パワー算出部２２６は、式（１９）を用いて、例えば現時刻から１波長前までの振動成分ｒ_ｔのＲＭＳ値を波形パワーＰＯＷとして算出する。なお、着目波長をＬ（ｓ）としたとき、ｔ_ｃ＝Ｌ／２となる。そして、波形パワー算出部２２６は、算出した波形パワーを危険判定部２５に出力する。

なお、第２、第３波形推定部２３、２４も第１波形推定部２２と同様に、第１波形推定部２２は、過去波形算出部２２１、揺動波形予測部２２２、平均波形推定部２２３、平均波形算出部２２４、揺動波形算出部２２５、波形パワー算出部２２６を有する。例えば、第２波形推定部２３は、中周期の波形について上記処理と同様の処理を実行し、中周期の振動成分を算出する。

（各方式の特性）
次に、各方式の特性について説明する。ここで説明する特性は、入力波形の振幅変化による平均波形の推定への影響を説明する。なお、短周期、中周期および長周期のいずれにも該当するので、第１平均波形等を単に平均波形などと記載して説明する。

図１３は、方式Ａの入力波形適合の特性を説明する図である。入力波形適合の特性は、低次多項式モデルで、期間Ｔ１では実際の平均波形を使い、Ｔ０では平均波形を通ることを条件とするため、期間Ｔ２では極端な変動は抑制され、期間Ｔ１での平均波形の形状の滑らかな延長となる。

具体的には、図１３の（１）に示すように、時刻Ｔ０を境に、入力波形の振幅が拡大する場合、振幅拡大に追従できず、延長傾向では結果として小さめの平均波形が推定され、推定される振幅は実際より大きめとなる傾向がある。一方、図１３の（２）に示すように、時刻Ｔ０を境に、入力波形の振幅が縮小する場合、元々の平均波形が持続するため、延長傾向で推定した平均波形とのかい離は少なく、推定される振幅と実際の振幅との差は少ない。

図１４は、方式Ｂの予測適合の特性を説明する図である。予測適合の特性は、期間Ｔ１での実際の振幅を期間Ｔ２でも維持するように入力波形の変化に沿った平均波形が推定されるため、結果として、振幅の変化に追従しない傾向となる。

具体的には、図１４の（１）に示すように、時刻Ｔ０を境に、入力波形の振幅が拡大する場合、期間Ｔ１での振幅を維持するように平均波形が推定され、入力波形に過度に沿った平均波形となり、推定される振幅は実際より小さめとなる傾向がある。一方、図１４の（２）に示すように、時刻Ｔ０を境に、入力波形の振幅が縮小する場合、期間Ｔ１での振幅を維持するように平均波形が推定され、振幅が縮小した入力波形から離れる平均波形となり、推定される振幅は実際より大きめとなる傾向がある。

上述したように、入力波形適合と予測適合は誤差傾向が対称的で、両者の和を取れば、誤差を相殺し、揺動波形の推定精度を高めることができる。図１５は、振幅拡大時のハイブリッド波形推定の特性を説明する図である。図１６は、振幅縮小時のハイブリッド波形推定の特性を説明する図である。

振幅拡大時は、図１５に示すように、入力波形適合では、振幅拡大に追従できず振幅が大きめに推定され、予測適合では、振幅の維持が働くので振幅が小さめに推定される。ハイブリッド波形推定では、振幅拡大時であっても、両方式の相反する誤差傾向を相殺できるので、平均波形の推定精度を向上させることができる。

また、振幅縮小時は、図１６に示すように、入力波形適合では、推定される振幅と実際の振幅との差は少ないが、予測適合では、推定される振幅は実際より大きめに推定される。ハイブリッド波形推定では、振幅縮小時であっても、両方式の誤差を低減できるので、平均波形の推定精度を向上させることができる。

（振幅拡大時の各方式のシミュレーション結果）
図１７は、方式Ａのシミュレーション結果（振幅拡大時）を説明する図である。図１８は、方式Ｂのシミュレーション結果（振幅拡大時）を説明する図である。図１９は、ハイブリッド波形推定のシミュレーション結果（振幅拡大時）を説明する図である。

図１７は、方式Ａによって、過去波形算出部２２１が算出した第１平均波形ｍ_ｔを用いて期間Ｔ２の第１平均波形ｙ１_ｔを推定した場合の推定結果と、推定した期間Ｔ２の第１平均波形ｙ１_ｔを用いて期間Ｔ２の第１揺動波形ｒ１_ｔの推定結果とを示す。図１７の（１）に示すように、期間Ｔ１から期間Ｔ２に対して入力波形ｐ_ｔの振幅が拡大した場合、予測された期間Ｔ２の平均波形ｙ１_ｔとシミュレーションで得られる期間Ｔ２の正解平均波形ｙ１_ｔ´との誤差が大きい。この結果、図１７の（２）に示すように、推定された期間Ｔ２の第１揺動波形ｒ１_ｔは、シミュレーションで得られる期間Ｔ２の正解揺動波形ｒ１_ｔ´の振幅よりも大きくなる（過拡大）。

図１８は、方式Ｂによって、過去波形算出部２２１が算出した第１平均波形ｍ_ｔおよび第２揺動波形ｒ２_ｔ´´を用いて期間Ｔ２の第１平均波形ｙ２_ｔを推定した場合の推定結果と、推定した期間Ｔ２の第１平均波形ｙ２_ｔを用いて期間Ｔ２の第１揺動波形ｒ２_ｔの推定結果とを示す。図１８の（１）に示すように、期間Ｔ１から期間Ｔ２に対して入力波形ｐ_ｔの振幅が拡大した場合、予測された期間Ｔ２の平均波形ｙ２_ｔとシミュレーションで得られる期間Ｔ２の正解平均波形ｙ２_ｔ´との誤差は大きい。この結果、図１８の（２）に示すように、推定された期間Ｔ２の第１揺動波形ｒ２_ｔは、シミュレーションで得られる期間Ｔ２の正解揺動波形ｒ２_ｔ´の振幅よりも小さくなる（過縮小）。

図１９は、期間Ｔ２において、図１７で得られる第１平均波形と図１８で得られる第１平均波形とを合算した図を示す。ハイブリット波形推定では、図１９の（１）に示すように、方式Ａを用いた図１７の（１）の第１平均波形ｙ１_ｔと、方式Ｂを用いた図１８の（１）の第１平均波形ｙ２_ｔとを加算するので、方式Ａと方式Ｂの中間程度の第１平均波形ｙ_ｔが推定される。この結果、図１９の（２）に示すように、方式Ａの誤差と方式Ｂの誤差とが相殺され、推定された期間Ｔ２の第１揺動波形ｒ_ｔの振幅と、シミュレーションで得られる期間Ｔ２の正解揺動波形ｒ_ｔ´の振幅との誤差は小さくなる。

（振幅縮小時の各方式のシミュレーション結果）
図２０は、方式Ａのシミュレーション結果（振幅縮小時）を説明する図である。図２１は、方式Ｂのシミュレーション結果（振幅縮小時）を説明する図である。図２２は、ハイブリッド波形推定のシミュレーション結果（振幅縮小時）を説明する図である。

図２０は、方式Ａによって、過去波形算出部２２１が算出した第１平均波形ｍ_ｔを用いて期間Ｔ２の第１平均波形ｙ１_ｔを推定した場合の推定結果と、推定した期間Ｔ２の第１平均波形ｙ１_ｔを用いて期間Ｔ２の第１揺動波形ｒ１_ｔの推定結果とを示す。図２０の（１）に示すように、期間Ｔ１から期間Ｔ２に対して入力波形ｐ_ｔの振幅が縮小した場合、予測された期間Ｔ２の平均波形ｙ１_ｔとシミュレーションで得られる期間Ｔ２の正解平均波形ｙ１_ｔ´との誤差は小さい。この結果、図２０の（２）に示すように、推定された期間Ｔ２の第１揺動波形ｒ１_ｔと、シミュレーションで得られる期間Ｔ２の正解揺動波形ｒ１_ｔ´との振幅もほぼ同じになる。

図２１は、方式Ｂによって、過去波形算出部２２１が算出した第１平均波形ｍ_ｔおよび第２揺動波形ｒ２_ｔ´´を用いて期間Ｔ２の第１平均波形ｙ２_ｔを推定した場合の推定結果と、推定した期間Ｔ２の第１平均波形ｙ２_ｔを用いて期間Ｔ２の第１揺動波形ｒ２_ｔを推定した場合の推定結果とを示す。図２１の（１）に示すように、期間Ｔ１から期間Ｔ２に対して入力波形ｐ_ｔの振幅が縮小した場合、予測された期間Ｔ２の平均波形ｙ２_ｔとシミュレーションで得られる期間Ｔ２の正解平均波形ｙ２_ｔ´との誤差は大きい。この結果、図２１の（２）に示すように、推定された期間Ｔ２の第１揺動波形ｒ２_ｔは、シミュレーションで得られる期間Ｔ２の正解揺動波形ｒ２_ｔ´の振幅よりも大きくなる（過拡大）。

図２２は、期間Ｔ２において、図２０で得られる第１平均波形と図２１で得られる第１平均波形とを合算した図を示す。ハイブリット波形推定では、図２２の（１）に示すように、方式Ａを用いた図２０の（１）の第１平均波形ｙ１_ｔと、方式Ｂを用いた図２１の（１）の第１平均波形ｙ２_ｔとを加算するので、方式Ａと方式Ｂの中間程度の第１平均波形ｙ_ｔが推定される。この結果、図２２の（２）に示すように、方式Ｂの誤差が方式Ａの誤差によって減少するので、推定された期間Ｔ２の第１揺動波形ｒ_ｔの振幅と、シミュレーションで得られる期間Ｔ２の正解揺動波形ｒ_ｔ´の振幅との誤差も、方式Ｂの誤差に比べて減少する。

このように、方式（Ａ）と方式（Ｂ）の両方を用いることで、定常時特性に優れる方式（Ｂ）と増減変化で逆傾向の特性を持つ方式（Ａ）とを組み合わせ、増減の追従性に優れた搖動推定を実現できる。

（波形パワーのシミュレーション結果）
次に、走行を開始してから現時刻までの短周期、中周期および長周期の成分波形それぞれにおける最終的に推定される波形パワーのシミュレーション結果について説明する。図２３は、ハイブリッド波形推定の波形パワー推定のシミュレーション結果を説明する図である。

図２３では、最大波長が１．５秒程度の短周期、最大波長が１６秒程度の中周期、最大波長が９０秒程度の長周期のそれぞれについて、上述した実施例１のハイブリッド波形推定で得られる波形パワーの波形と、リアルタイム処理ができないが正確な運転判定が実行できる蓄積手法で得られる波形パワーの波形との相関（波形相関値）を示す。図２３に示すように、いずれの周期であっても、各波形パワーの推定波形と、蓄積手法で得られる波形とが高相関になっており、危険ふらつき分析に十分な精度で予測できていることを確認できる。

図２４は、各方式の波形パワー推定のシミュレーション結果を説明する図である。図２４では、揺動波形を直接推定する手法、方式（Ａ）の入力波形適合のみを用いる手法、方式（Ｂ）の予測適合のみを用いる手法、ハイブリッド波形推定を用いる手法のそれぞれについて、最終的に推定される波形パワーの波形と蓄積手法で得られる波形パワーの波形との相関およびＲＭＳ誤差を示す。図２４に示すように、ハイブリッド波形推定と蓄積手法との相関が最も高く、ハイブリッド波形推定のＲＭＳ誤差が最も小さい。つまり、各手法と比較しても、ハイブリッド波形推定の危険ふらつき分析の精度が最も高いことを確認できる。

（危険判定部の機能ブロック）
図２５は、危険判定部２５の機能構成を示す機能ブロック図である。危険判定部２５は、例えば各周期における振動成分の特徴量に基づき、車両の運転に関する危険度を判定する判定部の一例である。危険判定部２５は、操舵評価部２５１、揺動評価部２５２、危険度評価部２５３および通知部２５４を有する。

操舵評価部２５１は、運転者の運転操舵に関し、急操舵を起こした状態に基づく運転操作を評価する処理部である。操舵評価部２５１は、短周期Ｔｓの波形パワーを第１波形推定部２２から取得する。また、操舵評価部２５１は、中周期の振動成分の波形値（振幅値）を取得する。そして、操舵評価部２５１は、取得した波形パワーおよび波形値に基づき、急操舵スコア値を算出し、算出した急操舵スコア値を危険度評価部２５３に出力する。

具体的には、操舵評価部２５１は、まず評価関数ｆ１（ｘ）を用いて、急操舵波形パワースコア値ｓ１を算出する。例えば、短周期Ｔｓの波形パワーをｐ１とすると、操舵評価部２５１はｓ１＝ｆ１（ｐ１）を算出する。なお、評価関数ｆ１（ｘ）は、式（２０−１）に示すように波形パワーをスコア値に線形変換する、あるいは式（２０−２）に示すように波形パワーをスコア値に非線形変換する関数である。なお、式（２０−１）のα１およびβ１は、定数であり、スコア値に変換する波形パワーの種別ごとに設定される。

また、操舵評価部２５１は、中周期の振動成分の波形値（通常操舵揺動波形値）ｗ２を用いて通常操舵揺動率ｔｖを算出する。具体的には、式（２１）に示すように、波形値ｗ２の絶対値を閾値距離ｂａｓｅで除算した値および下限値のうち小さい値を通常操舵揺動率ｔｖとする。なお、閾値距離ｂａｓｅは、予め設定される値であり、例えば２００ｍｍである。また、式（２１）では下限値を「１」としている。続いて、操舵評価部２５１は、急操舵波形パワースコア値ｓ１と通常操舵揺動率ｔｖとを乗算して、急操舵スコア値ｓｃ１＝ｓ１＊ｔｖを算出する。

揺動評価部２５２は、運転者の運転操舵に関し、走行軌跡の全体揺動の大きさを評価する処理部である。揺動評価部２５２は、第２、第３波形推定部２３、２４から取得した中周期および長周期の波形パワーに基づいて全体揺動スコア値を算出し、算出した全体揺動スコア値を危険度評価部２５３に出力する。

具体的には、揺動評価部２５２は、まず評価関数ｆ１（ｘ）を用いて、通常操舵波形パワースコア値Ｓ２および路内揺動波形パワースコア値Ｓ３を算出する。例えば、中周期の波形パワーをｐ２、長周期の波形パワーをｐ３とすると、揺動評価部２５２はＳ２＝ｆ１（ｐ２）およびＳ３＝ｆ１（ｐ３）を算出する。続いて、揺動評価部２５２は、式（２２）を用いてパワースコア高騰率ｇを算出する。

そして、揺動評価部２５２は、通常操舵波形パワースコア値Ｓ２と路内揺動波形パワースコア値Ｓ３との加算結果にパワースコア高騰率ｇを乗算して全体揺動スコア値ｓｃ２＝（Ｓ２＋Ｓ３）＊ｇを算出する。

危険度評価部２５３は、操舵評価部２５１の評価結果と揺動評価部２５２の評価結果とに基づき、運転者の運転の危険度を評価する処理部である。具体的には、まず、走行軌跡から車両が走行している横変位ｍ１を取得する。続いて、危険度評価部２５３は、評価関数ｆ２（ｘ）を用いて、横変位危険スコア値ｃｒを算出する。例えば、危険度評価部２５３は、ｃｒ＝ｆ２（ｍ１）を算出する。なお、評価関数ｆ２（ｘ）は、非線形変換を行う場合は、式（２３−１）に示すように評価関数ｆ１（ｘ）と同じ関数であり、線形変換を行う場合は、式（２３−２）に示す一次関数である。なお、式（２３−２）のα２およびβ２は定数である。

続いて、危険度評価部２５３は、急操舵スコア値ｓｃ１と全体揺動スコア値ｓｃ２との加算結果に横変位危険スコア値ｃｒを乗算して危険度スコア値ｓｃ３＝（ｓｃ１＋ｓｃ２）＊ｃｒを算出する。

そして、危険度評価部２５３は、算出した危険度スコア値ｓｃ３を通知部２５４に出力する。また、対象となった走行軌跡の時間と危険度スコア値ｓｃ３とを対応付けて、評価値ＤＢ１５に格納する。なお、危険度スコア値ｓｃ３に加え、急操舵スコア値ｓｃ１および全体揺動スコア値ｓｃ２を評価値ＤＢ１５に格納するようにしてもよい。

通知部２５４は、一定値を超える危険度スコア値ｓｃ３を危険度が高いと判定して、自動運転制御装置２に通知する処理部である。例えば、通知部２５４は、危険度評価部２５３から取得した危険度スコア値ｓｃ３が一定値を超える場合、危険ふらつきが発生したことを自動運転制御装置２に通知する。自動運転制御装置２は、通知部２５４からの通知を受けると、自動運転モードに切り替える。また、通知部２５４は、危険ふらつきが発生していることを示すメッセージをディスプレイ等に表示したり、音声で通知したりすることで、運転者に危険を通知するようにしてもよい。

［処理の流れ］
次に、各処理の流れについて説明する。ここでは、全体的な処理の流れ、波形推定処理の流れ、危険判定処理の流れについて説明する。

（全体的な処理の流れ）
図２６は、処理の流れを示すフローチャートである。図２６に示すように、第１波形推定部２２が短周期の波形推定処理を実行し（Ｓ１０１）、第２波形推定部２３が中周期の波形推定処理を実行する（Ｓ１０２）。続いて、第３波形推定部２４が長周期の波形推定処理を実行する（Ｓ１０３）。その後、危険判定部２５が、危険判定処理を実行し（Ｓ１０４）、処理を終了する。

（波形推定処理の流れ）
図２７は、波形推定処理の流れを示すフローチャートである。ここでは、第１波形推定部２２が行う短周期の波形推定処理について説明する。なお、第２、第３波形推定部２３、２４が実行する波形推定処理も短周期の波形推定処理と同様であるため、説明を省略する。

図２７に示すように、第１波形推定部２２は過去波形算出を実行する。具体的には、過去波形算出部２２１は、波長を算出し（Ｓ２０１）、第１平均波形を算出し（Ｓ２０２）、振動成分を算出する（Ｓ２０３）。続いて、過去波形算出部２２１は、波形パワーを算出し（Ｓ２０４）、算出した各波形等を波形ＤＢ１３等に格納する（Ｓ２０５）。

続いて、第１波形推定部２２は現在波形算出を実行する。具体的には、揺動波形予測部２２２または平均波形推定部２２３は、方式Ａによる第１平均波形を推定し（Ｓ２０６）、方式Ｂによる第１平均波形を推定する（Ｓ２０７）。

続いて、平均波形算出部２２４は、各方式で推定された各第１平均波形を合算する（Ｓ２０８）。その後、揺動波形算出部２２５は、合算された第１平均波形と入力波形とから、第１揺動波形を推定する（Ｓ２０９）。

そして、波形パワー算出部２２６は、Ｓ２０９で算出した第１揺動波形（振動成分）の波形パワーを算出する（Ｓ２１０）。その後、第１波形推定部２２は、各波形等を後段の処理部に出力し（Ｓ２１１）、処理を終了する。

（危険判定処理の流れ）
図２８は、危険判定処理の流れを示すフローチャートである。図２８に示すように、操舵評価部２５１は急操舵の評価値（急操舵スコア値ｓｃ１）を算出する（Ｓ３０１）。続いて、揺動評価部２５２は全体揺動の評価値（全体揺動スコア値ｓｃ２）を算出する（Ｓ３０２）。その後、危険度評価部２５３は、Ｓ３０１で算出された急操舵の評価値とＳ３０２で算出された全体揺動の評価値とに基づいて危険度の評価値（危険度スコア値ｓｃ３）を算出する（Ｓ３０３）。そして、通知部２５４は、危険度の評価値が閾値以上であれば（Ｓ３０４：Ｙｅｓ）、評価結果を自動運転制御装置２に出力し（Ｓ３０５）、評価値が閾値未満であれば（Ｓ３０４：Ｎｏ）、処理を終了する。

［効果］
上述したように、実施例１にかかる車載装置１０は、車道を走行する車両の走行軌跡に応じた入力波形と、入力波形から算出される算出済みの第１周波数未満の第１振動成分（平行波形）とを用いて、現時刻における第１振動成分を推定する。車載装置１０は、入力波形と、算出済みの第１振動成分と、入力波形から算出される算出済みの第１周波数以上の第２振動成分（揺動波形）から予測される現時刻における第２振動成分とを用いて、現時刻における第１振動成分を推定する。車載装置１０は、推定された各第１振動成分に基づいて、入力波形から現時刻における第２振動成分を算出する。

つまり、実施例１にかかる車載装置１０は、操舵行動に関連した波形分解の方法に着目する。そして、車載装置１０は、直流成分に相当する、変化が緩やかで低次パラメータで安定推定可能な平均波形を、直近１／２波長間でも観測される入力波形と、直近１／２波長の揺動波形から推定した推定揺動波形を拘束として利用し推定する。その後、車載装置１０は、「入力波形」−「平均波形」の算術演算で揺動波形を算出することで、現時刻の波形パワー値を高精度に推定することができる。

この結果、実施例１にかかる車載装置１０は、現時点の走行状態として、危険ふらつきの有無をリアルタイムに推定することができる。実施例１にかかる車載装置１０は、ドライブレコーダなど車載器を用いて、白線検出技術で得られる車線内の走行軌跡からリアルタイムに走行状態を判定し運転者にフィードバックすることができる。

さて、これまで本発明の実施例について説明したが、本発明は上述した実施例以外にも、種々の異なる形態にて実施されてよいものである。

［判定手法］
上記実施例で説明した危険判定部２５による評価値の算出は、具体的に記載した算出例に限定されるものではない。例えば、危険判定部２５が短周期の振動成分から路内揺動の位相等も用いて評価値を算出するようにしてもよい。

例えば図２５に示す危険判定部２５の操舵評価部２５１が、車両の横変位、路内揺動の位相および短周期Ｔｓの振動成分ｒｔの波形パワーを用いて運転者の急操舵の検出および評価を実行するようにしてもよい。

具体的には、操舵評価部２５１は、走行軌跡から車両が走行している横変位を取得する。また、操舵評価部２５１は、短周期の振動成分から路内揺動の位相を算出する。つまり、操舵評価部２５１は、長期的視野で運転者のハンドル操作（操舵）のタイミングを示す路内揺動の位相（図２のｄの変化）を算出する。ここで算出される位相は、ハンドル操作を行ったタイミングを示すものであることから、白線の実線や破線で示す道路端（揺動端）ほど危険度が高い操舵を示す。なお、操舵評価部２５１は、様々な公知の手法を用いて位相を算出することができる。

図２９は、操舵評価部２５１の処理を説明する図である。図２９に示すように、操舵評価部２５１は、横変位にフィルタをかけて正規化した値を算出する。同様に、操舵評価部２５１は、短周期の波形パワー（急操舵パワー）に対してパワーフィルタをかけて正規化し、路内揺動の位相に対し位相フィルタをかけて正規化する。そして、操舵評価部２５１は、正規化された各値を加算して、急操舵のスコア値を算出し、危険度評価部２５３に出力する。

ここで、操舵評価部２５１が用いる各フィルタを説明する。図３０は、急操舵の評価に使用するフィルタを説明する図である。図３０に示すフィルタは、サチュレートするフィルタ、すなわち大小の値でともに飽和するフィルタである。図３０の（１）に示すフィルタは、横変位量に対するフィルタであり、左右のうちより白線に近い側の距離により応答値を得るフィルタであり、白線との距離が短い場合に高い応答値を示すフィルタである。例えば、白線との距離（白線距離）が近い距離（ＴＨ１１）より長くなると、応答値が徐々に小さく変化し、閾値ＴＨ１２を超えると一定値に収束する。

図３０の（２）に示すフィルタは、路内揺動の位相に対するフィルタであり、路内揺動の位相θからｓｉｎθを算出し、より位相が波端の場合に高い応答値を示すフィルタである。例えば、波端の閾値ＴＨ２１を超えると応答値が徐々に大きくなり、閾値ＴＨ２２を超えると一定値に収束する。

図３０の（３）に示すフィルタは、急操舵パワーに対するフィルタであり、閾値ＴＨ３１を超えた場合に強く反応するフィルタである。例えば、閾値ＴＨ３１よりもパワーが大きい閾値ＴＨ３２の方が応答最大値に近い値が返される。なお、閾値ＴＨ１１〜ＴＨ３２の各閾値はそれぞれ個別の値で事前に定めておく。

そして、操舵評価部２５１は、フィルタ処理後の各値を積算することでスコア値を算出して急操舵を判定する。フィルタ処理後の各特徴量の積算は、すべてのフィルタ処理後の値が高い場合に高い評価値となる。図３１は、急操舵の判定を説明する図である。

図３１のａ５は走行軌跡の高周波成分、図３１のａ６は低周波成分を表しており、横変位量のフィルタ処理後の値が高くなるのは、図３１のＲｃの領域であり、路内揺動の位相のフィルタ処理後の値が高くなるのは図３１のＲｄの領域である。このため、急な操舵修正の評価結果がより明確に現れるのは、図３１のＲｅ１、Ｒｅ２の２領域である。この例では、Ｒｅ１、Ｒｅ２のうち、より高周波パワーが高いＲｅ２がより危険度が高いタイミングであると判定される。

つまり、操舵評価部２５１は、走行軌跡の中から、路内揺動の位相を取得することで、位相端（Ｒｄ）として４か所が特定できる。続いて、操舵評価部２５１は、特定した４か所のうち、車線内の横位置（ａ５）に基づいて車両が白線と近い位置となった地点Ｒｅ１とＲｅ２とを抽出する。その後、操舵評価部２５１は、地点Ｒｅ１とＲｅ２のうち、短周期の急操舵パワーが閾値以上である地点Ｒｅ２を特定する。この結果、操舵評価部２５１は、走行軌跡のうちＲｅ２で危険運転が発生したと評価する。このように、操舵評価部２５１は、路内揺動の位相と車線内の横位置から、逸脱を避けるための操舵修正タイミングに絞り込み、集中度低下に伴う修正操舵の急激性の評価を行う。

また、図２５に示す揺動評価部２５２が、長周期の振動成分の波形パワーと中周期の振動成分の波形パワーとを用いて、運転者の揺動を評価するようにしてもよい。つまり、揺動評価部２５２は、路内揺動の波形パワーと通常操舵の波形パワーとを加算して、危険度を示すスコア値を算出する。

図３２は、揺動評価部２５２の処理を説明する図である。図３２に示すように、揺動評価部２５２は、路内揺動の波形パワーにフィルタをかけて波形パワーを正規化し、同様に、通常操舵の波形パワーにフィルタをかけて波形パワーを正規化する。そして、揺動評価部２５２は、フィルタ処理後の各波形パワーを加算してスコア値を算出して、危険度評価部２５３に出力する。このようにすることで、揺動評価部２５２は、いずれか一方の波形パワーが高いときにスコア値が高くなることを除去し、両方の波形パワーが高いときにスコア値が高くなるように制御する。

図３３は、全体揺動の評価を説明する図である。図３３は、路内揺動の波形パワーと通常操舵の波形パワーの両パワーの時系列変化の例を示している。図３３の（１）は、路内揺動および通常操舵の２つの波形を示しており、ａ１１は、一方だけの波形パワーが高い状態であり、ｂ１１は、両方の波形パワーが高い状態を表している。

このような波形を単純に加算した状態を図３３の（５）に示す。図３３の（５）では、一方だけ非常に高いａ１５の状態も、両方とも少し高いｂ１５の状態も同様のスコア値となり、両パワーがともに高まっているかどうかの判別が難しい。つまり、このような処理は、癖がある運転者の通常操舵についても危険と判定するので、通常操舵とは異なる意識的な急操舵すなわち本当に危険な状態を特定することができない。

これに対して、揺動評価部２５２は、図３３の（１）に示す２つの波形に対して、（４）に示すフィルタをかける。このフィルタは、通常操舵と路内揺動でパラメータは異なるが、基本的な性質は同じである。すなわち低パワーの閾値ＴＨ１までは、最小値近傍値を応答し、高パワーの閾値ＴＨ２より最高値近傍値を応答し、閾値ＴＨ１、ＴＨ２間ではパワーの大きさに従い応答する性質を持つ。

そして、揺動評価部２５２が、図３３の（１）の波形に図３３の（４）のフィルタによってフィルタ処理することで、両方の波形パワーが高い状態（ｂ１１）の両パワーの応答値が高くなる（ｂ１２）。その後、揺動評価部２５２は、得られたフィルタ処理後の両パワーの単純加算を行うことで、図３３の（３）に示す波形を得ることができる。なお、適宜さらにパワーの比を重み係数にするなどしてもよい。図３３の（３）では、図３３の（２）で引き上げられたｂ１２が、加算の結果、ａ１３の領域に比べ高いスコアとなり、事前に定めた全体揺動の危険性基準のｃ１を超える。

このように、図３３の（４）に示すフィルタ処理後の波形パワーを評価することで、通常操舵と路内揺動の両波形パワーの大きさに加え、ともに高まっている状態を評価することが可能となり、ロバストに継続的な集中力低下状態を特定可能となる。なお、各フィルタの応答値を同範囲、例えば０〜１．０に設定することで、異なるパワーレンジの高まりを正規化して評価することが可能である。

さらに、図２５に示す危険度評価部２５３が、操舵評価部２５１により算出されたスコア値と揺動評価部２５２により算出されたスコア値とを加算した評価値を算出するようにしてもよい。

なお、揺動評価と操舵評価はいずれを先に実行してもよく、処理する順番は任意に変更することができる。また、揺動評価と操舵評価のうちいずれかの評価だけを実行してもよく、両方を実行してもよい。

［数値］
上記実施例で説明した各数値や閾値はあくまで例示であり、具体的に記載した値に限定されるものではない。

［各方式］
上記実施例では、方式Ａと方式Ｂとを組み合わせたハイブリッド波形推定について説明したが、これに限定されるものではなく、例えば方式Ｂのみなどのように、各方式だけを用いて揺動波形の推定を実行することもできる。

［システム］
また、図５に示した各装置の各構成は、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、任意の単位で分散または統合して構成することができる。例えば、第１〜第３波形推定部２２〜２４を統合することができる。さらに、各装置にて行なわれる各処理機能は、その全部または任意の一部が、ＣＰＵ（Central Processing Unit）および当該ＣＰＵにて解析実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現され得る。この他、上記文書中や図面中で示した処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて任意に変更することができる。

［ハードウェア構成］
上記車載装置１０は、例えば、次のようなハードウェア構成を有するコンピュータにより実現することができる。図３４は、ハードウェア構成例を説明する図である。図３４に示すように、車載装置１０は、通信インタフェース１０ａ、ＨＤＤ（Hard Disk Drive）１０ｂ、メモリ１０ｃ、プロセッサ１０ｄを有する。

通信インタフェース１０ａの一例としては、ネットワークインタフェースカードなどである。ＨＤＤ１０ｂは、図３等に示した各種ＤＢを記憶する記憶装置である。

メモリ１０ｃの一例としては、ＳＤＲＡＭ（Synchronous Dynamic Random Access Memory）等のＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）、フラッシュメモリ等が挙げられる。プロセッサ１０ｄの一例としては、ＣＰＵ、ＤＳＰ（Digital Signal Processor）、ＦＰＧＡ（Field Programmable Gate Array）、ＰＬＤ（Programmable Logic Device）等が挙げられる。

また、車載装置１０は、プログラムを読み出して実行することで波形推定方法を実行する情報処理装置として動作する。つまり、車載装置１０は、軌跡生成部２１、第１波形推定部２２、第２波形推定部２３、第３波形推定部２４および危険判定部２５と同様の機能を実行するプログラムを実行する。この結果、車載装置１０は、軌跡生成部２１、第１波形推定部２２、第２波形推定部２３、第３波形推定部２４および危険判定部２５と同様の機能を実行するプロセスを実行することができる。なお、この他の実施例でいうプログラムは、車載装置１０によって実行されることに限定されるものではない。例えば、他のコンピュータまたはサーバがプログラムを実行する場合や、これらが協働してプログラムを実行するような場合にも、本発明を同様に適用することができる。

このプログラムは、インターネットなどのネットワークを介して配布することができる。また、このプログラムは、ハードディスク、フレキシブルディスク（ＦＤ）、ＣＤ−ＲＯＭ、ＭＯ（Magneto−Optical disk）、ＤＶＤ（Digital Versatile Disc）などのコンピュータで読み取り可能な記録媒体に記録され、コンピュータによって記録媒体から読み出されることによって実行することができる。

以上の各実施例を含む実施形態に関し、さらに以下の付記を開示する。

（付記１）車道を走行する車両の走行軌跡に応じた入力波形と、前記入力波形から算出される算出済みの第１周波数未満の第１振動成分とを用いて、現時刻における第１振動成分を推定する第１推定部と、
前記入力波形と、算出済みの前記第１振動成分と、前記入力波形から算出される算出済みの第１周波数以上の第２振動成分から予測される前記現時刻における第２振動成分とを用いて、前記現時刻における第１振動成分を推定する第２推定部と、
推定された各前記第１振動成分に基づいて、前記入力波形から前記現時刻における第２振動成分を算出する算出部と
を有することを特徴とする波形推定装置。

（付記２）前記現時刻以前の前記入力波形から、前記現時刻の一波長前から半波長前までの第１期間の前記第１振動成分を算出する振動算出部をさらに有し、
前記第１推定部は、前記半波長前から前記現時刻までの第２期間の前記入力波形と、前記振動算出部によって算出された前記第１期間の前記第１振動成分とを用いて、前記第２期間の前記第１振動成分を推定することを特徴とする付記１に記載の波形推定装置。

（付記３）前記第１期間の前記入力波形から、前記第１期間の前記第１振動成分を減算して、前記第１期間の前記第２振動成分を算出し、前記第１期間の前記第２振動成分の振幅および位相に基づいて、前記第２期間の前記第２振動成分を予測する予測部をさらに有し、
前記第２推定部は、前記第２期間の前記入力波形と、前記第１期間の前記第１振動成分と、前記予測部によって予測された前記第２期間の前記第２振動成分とから、前記第２期間の前記第１振動成分を推定することを特徴とする付記２に記載の波形推定装置。

（付記４）前記算出部は、前記第１推定部によって推定された前記第２期間の前記第１振動成分と、前記第２推定部によって推定された前記第２期間の前記第１振動成分とを合算して、前記第２期間の前記第１振動成分を算出し、算出した前記第２期間の前記第１振動成分を前記第２期間の前記入力波形から減算して、前記第２期間の前記第２振動成分を算出することを特徴とする付記３に記載の波形推定装置。

（付記５）前記第２振動成分の特徴量に基づき、前記車両の運転に関する危険度を判定する判定部をさらに有することを特徴とする付記１〜４のいずれか一つに記載の波形推定装置。

（付記６）コンピュータが、
車道を走行する車両の走行軌跡に応じた入力波形と、前記入力波形から算出される算出済みの第１周波数未満の第１振動成分とを用いて、現時刻における第１振動成分を推定し、
前記入力波形と、算出済みの前記第１振動成分と、前記入力波形から算出される算出済みの第１周波数以上の第２振動成分から予測される前記現時刻における第２振動成分とを用いて、前記現時刻における第１振動成分を推定し、
推定された各前記第１振動成分に基づいて、前記入力波形から前記現時刻における第２振動成分を算出する
処理を実行することを特徴とする波形推定方法。

（付記７）コンピュータに、
車道を走行する車両の走行軌跡に応じた入力波形と、前記入力波形から算出される算出済みの第１周波数未満の第１振動成分とを用いて、現時刻における第１振動成分を推定し、
前記入力波形と、算出済みの前記第１振動成分と、前記入力波形から算出される算出済みの第１周波数以上の第２振動成分から予測される前記現時刻における第２振動成分とを用いて、前記現時刻における第１振動成分を推定し、
推定された各前記第１振動成分に基づいて、前記入力波形から前記現時刻における第２振動成分を算出する
処理を実行させることを特徴とする波形推定プログラム。

（付記８）車両に搭載される車載装置および自動運転制御装置を含む運転制御システムにおいて、
前記車載装置は、
車道を走行する車両の走行軌跡に応じた入力波形と、前記入力波形から算出される算出済みの第１周波数未満の第１振動成分とを用いて、現時刻における第１振動成分を推定する第１推定部と、
前記入力波形と、算出済みの前記第１振動成分と、前記入力波形から算出される算出済みの第１周波数以上の第２振動成分から予測される前記現時刻における第２振動成分とを用いて、前記現時刻における第１振動成分を推定する第２推定部と、
推定された各前記第１振動成分に基づいて、前記入力波形から前記現時刻における第２振動成分を算出する算出部と、
前記第２振動成分の特徴量に基づき、前記車両の運転に関する危険度を判定する判定部と、を有し、
前記自動運転制御装置は、
前記車載装置によって判定された前記危険度に基づいて、前記車両の運転を制御する運転制御部を有することを特徴とする運転制御システム。

２ 自動運転制御装置
１０ 車載装置
２１ 軌跡生成部
２２ 第１波形推定部
２３ 第２波形推定部
２４ 第３波形推定部
２５ 危険判定部
２２１ 過去波形算出部
２２２ 揺動波形予測部
２２３ 平均波形推定部
２２４ 平均波形算出部
２２５ 揺動波形算出部
２２６ 波形パワー算出部

Claims (7)

1. 車道を走行する車両の走行軌跡に応じた入力波形と、前記入力波形から算出される算出済みの第１周波数未満の第１振動成分とを用いて、現時刻における第１振動成分を推定する第１推定部と、
   前記入力波形と、算出済みの前記第１振動成分と、前記入力波形から算出される算出済みの第１周波数以上の第２振動成分から予測される前記現時刻における第２振動成分とを用いて、前記現時刻における第１振動成分を推定する第２推定部と、
   推定された各前記第１振動成分に基づいて、前記入力波形から前記現時刻における第２振動成分を算出する算出部と
   を有することを特徴とする波形推定装置。
2. 前記現時刻以前の前記入力波形から、前記現時刻の一波長前から半波長前までの第１期間の前記第１振動成分を算出する振動算出部をさらに有し、
   前記第１推定部は、前記半波長前から前記現時刻までの第２期間の前記入力波形と、前記振動算出部によって算出された前記第１期間の前記第１振動成分とを用いて、前記第２期間の前記第１振動成分を推定することを特徴とする請求項１に記載の波形推定装置。
3. 前記第１期間の前記入力波形から、前記第１期間の前記第１振動成分を減算して、前記第１期間の前記第２振動成分を算出し、前記第１期間の前記第２振動成分の振幅および位相に基づいて、前記第２期間の前記第２振動成分を予測する予測部をさらに有し、
   前記第２推定部は、前記第２期間の前記入力波形と、前記第１期間の前記第１振動成分と、前記予測部によって予測された前記第２期間の前記第２振動成分とから、前記第２期間の前記第１振動成分を推定することを特徴とする請求項２に記載の波形推定装置。
4. 前記算出部は、前記第１推定部によって推定された前記第２期間の前記第１振動成分と、前記第２推定部によって推定された前記第２期間の前記第１振動成分とを合算して、前記第２期間の前記第１振動成分を算出し、算出した前記第２期間の前記第１振動成分を前記第２期間の前記入力波形から減算して、前記第２期間の前記第２振動成分を算出することを特徴とする請求項３に記載の波形推定装置。
5. 前記第２振動成分の特徴量に基づき、前記車両の運転に関する危険度を判定する判定部をさらに有することを特徴とする請求項１〜４のいずれか一項に記載の波形推定装置。
6. コンピュータが、
   車道を走行する車両の走行軌跡に応じた入力波形と、前記入力波形から算出される算出済みの第１周波数未満の第１振動成分とを用いて、現時刻における第１振動成分を推定し、
   前記入力波形と、算出済みの前記第１振動成分と、前記入力波形から算出される算出済みの第１周波数以上の第２振動成分から予測される前記現時刻における第２振動成分とを用いて、前記現時刻における第１振動成分を推定し、
   推定された各前記第１振動成分に基づいて、前記入力波形から前記現時刻における第２振動成分を算出する
   処理を実行することを特徴とする波形推定方法。
7. コンピュータに、
   車道を走行する車両の走行軌跡に応じた入力波形と、前記入力波形から算出される算出済みの第１周波数未満の第１振動成分とを用いて、現時刻における第１振動成分を推定し、
   前記入力波形と、算出済みの前記第１振動成分と、前記入力波形から算出される算出済みの第１周波数以上の第２振動成分から予測される前記現時刻における第２振動成分とを用いて、前記現時刻における第１振動成分を推定し、
   推定された各前記第１振動成分に基づいて、前記入力波形から前記現時刻における第２振動成分を算出する
   処理を実行させることを特徴とする波形推定プログラム。

Images