(57)【要約】【課題】波形の信号処理にかかる時間を短縮すること。【解決手段】波形推定装置は、推定部と、算出部とを有する。推定部は、車道を走行する車両の走行軌跡に応じた入力波形において、現時刻から第１周波数の半波長前までの期間の入力波形を用いて、期間における第１周波数未満の第１振動成分を推定する。算出部は、入力波形から第１振動成分を減算して、期間における第１周波数以上の第２振動成分を算出する。【選択図】図３
2018184092.tif 000002
【特許請求の範囲】
【請求項１】
車道を走行する車両の走行軌跡に応じた入力波形において、現時刻から第１周波数の半波長前までの期間の前記入力波形を用いて、前記期間における前記第１周波数未満の第１振動成分を推定する推定部と、
前記第１振動成分に基づいて、前記入力波形から前記期間における前記第１周波数以上の第２振動成分を算出する算出部と、
を有することを特徴とする波形推定装置。
【請求項２】
前記算出部は、前記入力波形から前記第１振動成分を減算して前記第２振動成分を算出することを特徴とする請求項１に記載の波形推定装置。
【請求項３】
前記現時刻までの前記入力波形を用いて、前記第１周波数の半波長より以前の前記第１振動成分を算出する振動算出部をさらに有し、
前記推定部は、前記第１周波数の半波長より以前の前記第１振動成分を用いて、前記期間における前記第１振動成分を推定することを特徴とする請求項１または２に記載の波形推定装置。
【請求項４】
前記推定部は、前記期間における前記第１振動成分を前記期間における前記入力波形との一致度に基づいて推定することを特徴とする請求項１〜３のいずれか一項に記載の波形推定装置。
【請求項５】
前記期間における前記第１振動成分と前記期間における前記入力波形との差が最小化する多項式モデルの係数を算出する係数算出部をさらに有し、
前記推定部は、前記多項式モデルの前記係数から前記期間における前記第１振動成分を推定することを特徴とする請求項１〜４のいずれか一項に記載の波形推定装置。
【請求項６】
現時刻から前記第１周波数より低い第２周波数の半波長前までの第２期間の前記第１振動成分を用いて、前記第２期間における前記第２周波数未満の第３振動成分を推定する第２推定部と、
前記第１振動成分から前記第３振動成分を減算して、前記第２期間における前記第２周波数以上前記第１周波数未満の第４振動成分を算出する第２算出部と、
をさらに有することを特徴とする請求項１〜５のいずれか一項に記載の波形推定装置。
【請求項７】
前記第２振動成分の特徴量に基づき、前記車両の運転に関する危険度を判定する判定部をさらに有することを特徴とする請求項１〜６のいずれか一項に記載の波形推定装置。
【請求項８】
コンピュータが、
車道を走行する車両の走行軌跡に応じた入力波形において、現時刻から第１周波数の半波長前までの期間の前記入力波形を用いて、前記期間における前記第１周波数未満の第１振動成分を推定し、
前記入力波形から前記第１振動成分を減算して、前記期間における前記第１周波数以上の第２振動成分を算出する
処理を実行することを特徴とする波形推定方法。
【請求項９】
コンピュータに、
車道を走行する車両の走行軌跡に応じた入力波形において、現時刻から第１周波数の半波長前までの期間の前記入力波形を用いて、前記期間における前記第１周波数未満の第１振動成分を推定し、
前記入力波形から前記第１振動成分を減算して、前記期間における前記第１周波数以上の第２振動成分を算出する
処理を実行させることを特徴とする波形推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、波形推定装置、波形推定方法および波形推定プログラムに関する。
【背景技術】
【０００２】
従来から、運転への集中力低下で発生する、走行中の危険なふらつき挙動を検出してフィードバックすることで、安全運転を支援する装置やシステムが知られている。例えば、操舵角や横変位の時間変化による波形を高周波成分と低周波成分に分類し、低周波成分を高周波成分で除算したパワー比が高いときに低覚醒状態と判定する技術が知られている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２０００−１８５５７５号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかしながら、上記の技術では、波形の信号処理に時間がかかる。例えば、上記技術では、覚醒度を精度良く推定するために、５０秒から８０秒程度の長い時間、変位量データの蓄積が必要となり、波形分類のための信号処理に時間がかかってしまう。
【０００５】
１つの側面では、波形の信号処理にかかる時間を短縮することができる波形推定装置、波形推定方法および波形推定プログラムを提供することを目的とする。
【課題を解決するための手段】
【０００６】
第１の案では、波形推定装置は、推定部と、算出部とを有する。推定部は、車道を走行する車両の走行軌跡に応じた入力波形において、現時刻から第１周波数の半波長前までの期間の入力波形を用いて、期間における第１周波数未満の第１振動成分を推定する。算出部は、入力波形から第１振動成分を減算して、期間における第１周波数以上の第２振動成分を算出する。
【発明の効果】
【０００７】
一実施形態によれば、波形の信号処理にかかる時間を短縮することができる。
【図面の簡単な説明】
【０００８】
【図１】図１は、実施例１にかかる走行制御システムの機能構成を示す機能ブロック図である。
【図２】図２は、車両の走行軌跡を説明する図である。
【図３】図３は、実施例１にかかる車載装置の機能構成を示す機能ブロック図である。
【図４】図４は、走行軌跡の各成分波形を説明する図である。
【図５】図５は、振動成分の抽出時刻を説明する図である。
【図６】図６は、第１波形推定部の機能構成を示す機能ブロック図である。
【図７】図７は、第１平均波形の推定を説明する図である。
【図８】図８は、第１平均波形の推定を説明する図である。
【図９】図９は、短周期の振動成分の算出を説明する図である。
【図１０】図１０は、第１平均波形推定のシミュレーション結果の一例を示す図である。
【図１１】図１１は、第１平均波形推定のシミュレーション結果の一例を示す図である。
【図１２】図１２は、危険判定部の機能構成を示す機能ブロック図である。
【図１３】図１３は、処理の流れを示すフローチャートである。
【図１４】図１４は、波形推定処理の流れを示すフローチャートである。
【図１５】図１５は、危険判定処理の流れを示すフローチャートである。
【図１６】図１６は、操舵評価部の処理を説明する図である。
【図１７】図１７は、急操舵の評価に使用するフィルタを説明する図である。
【図１８】図１８は、急操舵の判定を説明する図である。
【図１９】図１９は、揺動評価部の処理を説明する図である。
【図２０】図２０は、全体揺動の評価を説明する図である。
【図２１】図２１は、ハードウェア構成例を説明する図である。
【発明を実施するための形態】
【０００９】
以下に、本願の開示する波形推定装置、波形推定方法および波形推定プログラムの実施例を図面に基づいて詳細に説明する。なお、この実施例によりこの発明が限定されるものではない。
【実施例１】
【００１０】
［全体構成］
図１は、実施例１にかかる走行制御システムの機能構成を示す機能ブロック図である。図１に示すように、走行制御システムは、ＧＰＳセンサ１と、車載装置１０と、自動運転制御装置２と、アクチュエータ３とを有する。走行制御システムは例えば車両に搭載される。
【００１１】
ＧＰＳ（Global  Positioning  System）センサ１は、車両に搭載されるセンサであり、車両の位置情報を取得する。ＧＰＳセンサ１は、取得した位置情報を車載装置１０および自動運転制御装置２に出力する。
【００１２】
車載装置１０は、車両に搭載されるコンピュータであり、ＧＰＳセンサ１から位置情報を取得し、車両の走行軌跡などを生成する。また、車載装置１０は、例えば白線検出処理などの公知の手法により、車両の走行軌跡を生成することもできる。また、車載装置１０は、走行軌跡を生成するために使用される車両の運行情報などを様々な公知の手法やＧＰＳセンサ１等により取得することもできる。
【００１３】
また、車載装置１０は、車両の走行軌跡や運行情報に基づき、運転者の運転を評価する。例えば、車載装置１０は、走行軌跡の波形を、無意識の舵取りの大きさが表れる長周期、意識的な舵取りの大きさが表れる中周期、突発的な舵取りの大きさが表れる短周期の３種の成分波形に分離する。このとき、車載装置１０は、各周期の半波長前から現時刻までの各成分波形を推定することで、走行軌跡の波形を現時刻までの各成分波形に分離することができる。そのため、車載装置１０は、波形を各成分波形に分離する処理（信号処理）にかかる時間を短縮することができる。
【００１４】
続いて、車載装置１０は、各波形の特徴量から、全体揺動および急操舵を評価する。車載装置１０は、走行軌跡の波形を現時刻までの各成分波形に分離することができるため、各評価をリアルタイムで行うことができる。そして、車載装置１０は、リアルタイムでの評価結果からふらつきの危険度を算出し、閾値を超えた場合に危険ふらつきと判定するので、危険な車両のふらつきを検出するまでの時間を短縮することができる。
【００１５】
ここで、走行軌跡と検出ポイントとの関係について説明する。図２は、車両の走行軌跡を説明する図である。図２に示すように、車両は、白線の破線と白線の実線との間の斜線を走行する。車両が破線または実線を外れて走行すると、他車線を走行する他車両などとの衝突などの危険が発生する。
【００１６】
例えば、図２に示す走行軌跡ａにおいて、危険ふらつきとは、全体揺動（図２のｂ）が高まっており、かつ、急操舵軌跡（図２のｃ）が道路端（図２のｄ）かつ揺動端で発生している状態である。ここで、全体揺動の高まりは車線維持集中力の低下、急操舵軌跡は反応遅れの突発的対応、道路端での急操舵は逸脱を回避する操舵タイミングを示し、車載装置１０はこれを評価することで危険度を判定する。車載装置１０は、判定結果を自動運転制御装置２に出力する。
【００１７】
自動運転制御装置２は、車両に搭載されるコンピュータであり、車両の自動運転を制御する。自動運転制御装置２は、例えば車載装置１０が危険ふらつきと判定した場合に、車両の運転モードを運転者が手動運転を行う手動運転モードから、車両を自動的に走行させる自動運転モードに切り替える。また、自動運転制御装置２は、例えば運転者からの指示にしたがって運転モードを切り替えることもできる。自動運転制御装置２は、自動運転モードに切り替えると、ＧＰＳセンサ１からの位置情報等に基づき、アクチュエータ３を制御して舵角や車速を制御することで、車両が自動的に走行するように自動運転制御を行う。また、自動運転制御装置２は、例えば白線検出処理などの公知の手法により、車両の自動運転制御を行うこともできる。
【００１８】
アクチュエータ３は、例えばアクセルアクチュエータ、ブレーキアクチュエータおよび操舵アクチュエータ等を含む。アクセルアクチュエータは、エンジンのスロットルバルブの開度を調整する。アクセルアクチュエータは、自動運転制御装置２からの指示に従い、スロットルバルブの開度を調整する。ブレーキアクチュエータは、各車輪のホイールシリンダそれぞれに油圧を加えることで、各車輪と共に回転するディスクやドラム等の回転を抑制する。ブレーキアクチュエータは、自動運転制御装置２からの指示に従い、油圧を調整する。操舵アクチュエータは、ステアリング機構に操舵トルクを付与する。操舵アクチュエータは、自動運転制御装置２からの指示に従い、操舵トルクを調整する。
【００１９】
［機能構成］
図３は、実施例１にかかる車載装置１０の機能構成を示す機能ブロック図である。なお、自動運転制御装置２は、車両の自動運転を行うコンピュータであり、一般的な自動運転制御装置と同様の機能を有するので、ここでは詳細な説明は省略する。
【００２０】
図３に示すように、車載装置１０は、記憶部１２および制御部２０を有する。記憶部１２は、プログラムやデータなどを記憶する記憶装置であり、例えばメモリやハードディスクなどである。この記憶部１２は、軌跡ＤＢ１３と、波形ＤＢ１４と、評価値ＤＢ１５とを記憶する。軌跡ＤＢ１３は、制御部２０が生成した走行軌跡や運行情報などを走行した時間と対応付けて記憶するデータベースである。
【００２１】
波形ＤＢ１４は、走行軌跡の波形を長周期、中周期および短周期に分離した各成分波形（振動成分）を記憶する。波形ＤＢ１４は、制御部２０が生成した各成分波形を走行した時間と対応付けて記憶するデータベースである。
【００２２】
ここで、走行軌跡の波形を分離した各成分波形について説明する。図４は、走行軌跡の各成分波形を説明する図である。図４に示すように、走行軌跡の波形ａは、短周期の急操舵に関する振動成分ａ１、中周期の通常操舵に関する振動成分ａ２および長周期の路内揺動に関する振動成分ａ３に分離される。各振動成分ａ１〜ａ３は、各時刻で周期が異なるが、その周期はそれぞれ定めた一定範囲内となる。
【００２３】
例えば、短周期には１〜４秒、中周期には８〜１６秒、長周期には６０〜９０秒の周期波形が含まれる。周期範囲は一例で有り、周期範囲が重複しなければこれに限らずともよい。振動成分ａ１〜ａ３における各時刻の周期長は、例えば走行軌跡の波形を、抽出時刻を中心としてＤＦＴ（離散フーリエ変換）に代表される周波数領域への変換を行い、対象とする周波数範囲を代表する応答成分で決定することができる。代表応答周期（周波数）の決定は、周波数範囲で最大応答を示した周波数（周期）を用いるか、または、周波数範囲における平均値、中間値などの統計量を用いて決定してもよい。
【００２４】
図３に戻る。評価値ＤＢ１５は、危険度の評価結果を記憶するデータベースである。例えば、評価値ＤＢ１５は、評価値を走行した時間と対応付けて記憶する。ここで、評価値は、走行軌跡に応じて評価された値であり、数値が高いほど危険度が高い運転を示す。
【００２５】
制御部２０は、車載装置１０全体を司る処理部であり、例えばプロセッサなどである。この制御部２０は、軌跡生成部２１、第１波形推定部２２、第２波形推定部２３、第３波形推定部２４および危険判定部２５を有する。なお、軌跡生成部２１、第１波形推定部２２、第２波形推定部２３、第３波形推定部２４および危険判定部２５は、プロセッサなどの電子回路の一例やプロセッサが実行するプロセスの一例である。
【００２６】
軌跡生成部２１は、ＧＰＳセンサ１から取得した位置情報にしたがって、車両の走行軌跡を生成する処理部である。また、軌跡生成部２１は、ドライブレコーダ等で撮像された走行時の画像情報等を受信し、白線検出処理などを用いて車両の走行軌跡を生成することもできる。軌跡生成部２１は、生成した走行軌跡を軌跡ＤＢ１３に格納する。
【００２７】
第１〜第３波形推定部２２〜２４は、走行を開始してから現時刻までの短周期、中周期および長周期の成分波形（振動成分ａ１〜ａ３）をそれぞれ走行軌跡の波形から分離する処理部である。第１波形推定部２２は、走行軌跡の波形を入力波形として、短周期の振動成分ａ１とそれ以外の振動成分（第１平均波形）とに分離する。また、第１波形推定部２２は、短周期の振動成分ａ１の波形パワーを算出する。なお、第１波形推定部２２の処理の詳細は、図６を用いて後述する。
【００２８】
続いて、第２波形推定部２３は、第１波形推定部２２が分離した第１平均波形を入力波形として、中周期の振動成分ａ２とそれ以外の振動成分（第２平均波形）とに分離する。また、第２波形推定部２３は、中周期の振動成分ａ２の波形パワーを算出する。第３波形推定部２４は、第２波形推定部２３が分離した第２平均波形を入力波形として低周期の振動成分ａ３に分離する。また、第３波形推定部２４は、低周期の振動成分ａ３の波形パワーを算出する。このように、第１〜第３波形推定部２２〜２４は、高周波から順に振動成分ａ１〜ａ３を分離する。
【００２９】
ここで、各振動成分ａ１〜ａ３の抽出時刻について説明する。図５は、振動成分ａ１〜ａ３の抽出時刻を説明する図である。なお、図５では、各振動成分ａ１〜ａ３の１波長を模式的に示している。第１〜第３波形推定部２２〜２４は、抽出時刻を中心として１波長の入力波形を周波数領域へ変換して振動成分ａ１〜ａ３を分離する。そのため、第１〜第３波形推定部２２〜２４は、現時刻から半波長前までの振動成分ａ１〜ａ３を走行軌跡の波形から直接分離することができるが、半波長前から現時刻までの振動成分ａ１〜ａ３を直接分離することができない。そのため、例えば１．５秒の短周期の場合、現時刻から０．２５秒前までの間、振動成分ａ１〜ａ３が分離されず、例えば１６秒の中周期の場合、現時刻から８秒前までの間、振動成分ａ１〜ａ３が分離されない。また、例えば９０秒の長周期の場合、現時刻から４５秒前までの間、振動成分ａ１〜ａ３が分離されない。
【００３０】
このように、第１〜第３波形推定部２２〜２４が走行軌跡の波形から直接振動成分ａ１〜ａ３を分離する場合、半波長に相当する時間の経過を待ってからしか振動成分ａ１〜ａ３を分離することができない。したがって、第１〜第３波形推定部２２〜２４が直接分離した振動成分ａ１〜ａ３を用いて、後述する危険判定部２５が危険ふらつき判定を行うとすると、現時刻から９０秒前までの危険ふらつきしか判定を行うことができず、リアルタイムで判定を行うことができない。そこで、本実施例では、第１〜第３波形推定部２２〜２４が現時刻から半波長前までの振動成分ａ１〜ａ３を推定することで、リアルタイムで振動成分の分離を行えるようにした。これにより、危険判定部２５がリアルタイムで危険ふらつきを判定することができる。
【００３１】
図６は、第１波形推定部２２の機能構成を示す機能ブロック図である。なお、第２、第３波形推定部２３、２４は、入力波形および分離する振動成分の周期が異なる点を除き、第１波形推定部２２と同様の処理を行うため、説明を省略する。
【００３２】
図６に示すように、第１波形推定部２２は、過去波形算出部２１１、パラメータ推定部２１２、平均波形算出部２１３、振動成分算出部２１４および波形パワー算出部２１５を有する。
【００３３】
過去波形算出部２１１は、入力波形（走行軌跡の波形）から、現時刻の半波長前までの振動成分と第１平均波形とを分離する処理部である。すなわち、過去波形算出部２１１は、現時刻までの入力波形を用いて、半波長より以前の第１平均波形を算出する振動算出部の一例である。具体的には、過去波形算出部２１１は、例えば全ｎ種の波長候補Ｌｃ＝｛Ｌ１、Ｌ２、・・・、Ｌｎ｝に対して、現時刻から半波長前の時刻を中心とする入力波形についてフーリエ変換を行う。そして、過去波形算出部２１１は、フーリエ変換係数のｓｉｎ値およびｃｏｓ値のベクトルサイズが最大となる波長を、現時刻から半波長前の時刻の波長に決定する。
【００３４】
続いて、過去波形算出部２１１は、決定した波長範囲で、入力波形を平均し、第１平均波形を算出する。過去波形算出部２１１は、現在時刻から半波長前までの入力波形から算出した第１平均波形を減算して振動成分を算出する。また、過去波形算出部２１１は、現時刻の１波長前の時刻を中心とする振動成分について、当該時刻で決定した波長範囲でＲＭＳ（Root  Mean  Square）値を波形パワーとして算出する。過去波形算出部２１１は、決定した波長と、分離した振動成分および第１平均波形と、波形パワーとを波形ＤＢ１４に格納する。
【００３５】
パラメータ推定部２１２は、第１平均波形を表す多項式の各パラメータを推定する処理部である。ここで、パラメータ推定部２１２が第１平均波形のパラメータを推定する理由について説明する。図７は、第１平均波形の推定を説明する図である。図７に示す振動成分ａ１は、第１平均波形ａ４より約４倍周期が短い。そのため、現時刻から短周期の半波長前の時刻までの間（期間Ｔ１）の第１平均波形ａ４は８分の１波長程度であり、振動成分ａ１は半波長程度である。したがって、期間Ｔ１の第１平均波形ａ４は、単調増加または単調減少する多項式でモデル化することができ、式（１）に示すように、振動成分ａ１を表す多項式より低次の多項式（二次多項式）でモデル化することができる。
【００３６】
【数１】
2018184092.tif 000003
【００３７】
なお、式（１）のｙｔは時刻ｔにおける第１平均波形であり、α、β、γはそれぞれパラメータである。また、現時刻から半波長前の時刻を基準（ｔ＝０）とする。したがって、現時刻は、短周期Ｔｓの２分の１の時刻（ｔ＝Ｔｓ／２）となる。パラメータγは、ｔ＝０のときの第１平均波形ｙ０であり、式（２）に基づいて算出される。なお、ｔ＝０のときの第１平均波形ｙ０は、過去波形算出部２１１が算出した第１平均波形ｑｔから得られる（ｙ０＝ｑ０）。
【００３８】
【数２】
2018184092.tif 000004
【００３９】
パラメータ推定部２１２は、期間Ｔ１において、第１平均波形ｙｔと入力波形ｐｔとの差が最小化するように式（１）の各パラメータα、βを算出する。また、パラメータ推定部２１２は、基準時刻ｔ＝０から半波長前までの期間Ｔ２において、第１平均波形ｙｔと過去波形算出部２１１が算出した第１平均波形ｑｔとの差が最小化するようにパラメータα、βを算出する。具体的に、パラメータ推定部２１２は、式（３）に示すように、第１平均波形ｙｔと波形ｍｔとの差を算出する。なお、式（４）に示すように、波形ｍｔは、ｔ≧０で入力波形ｐｔ、ｔ＜０で過去波形算出部２１１が算出した第１平均波形ｑｔとなる波形である。
【００４０】
【数３】
2018184092.tif 000005
【００４１】
続いて、パラメータ推定部２１２は、残差ｓｔの二乗残差和Ｓを式（５）に基づいて算出し、二乗残差和Ｓが最小となるように式（６）、式（７）を解くことで、式（８）に示すようにパラメータα、βを算出する。
【００４２】
【数４】
2018184092.tif 000006
【００４３】
【数５】
2018184092.tif 000007
【００４４】
【数６】
2018184092.tif 000008
【００４５】
【数７】
2018184092.tif 000009
【００４６】
図８は、第１平均波形の推定を説明する図である。パラメータ推定部２１２は、図８に示すように、期間Ｔ１の入力波形ｐｔに加え、期間Ｔ２の第１平均波形ｑｔを用いて期間Ｔ１の第１平均波形ｙｔを推定する。これにより、第１平均波形ｙｔをより高精度に推定することができる。このように、パラメータ推定部２１２は、期間Ｔ１における第１平均波形ｙｔと期間Ｔ１の入力波形ｐｔおよび期間Ｔ２の第１平均波形ｑｔとの差が最小かする多項式モデルの係数を算出する係数算出部の一例である。
【００４７】
図６に戻り、パラメータ推定部２１２は、算出したパラメータα〜γを平均波形算出部２１３に出力する。
【００４８】
平均波形算出部２１３は、パラメータ推定部２１２が算出したパラメータα〜γに基づき、期間Ｔ１の第１平均波形ｙｔを算出する処理部である。すなわち、平均波形算出部２１３は、期間Ｔ１の入力波形ｐｔを用いて、期間Ｔ１の第１平均波形ｙｔを推定する推定部の一例である。また、平均波形算出部２１３は、期間Ｔ１以前の第１平均波形ｑｔと期間Ｔ１の第１平均波形ｙｔを現時刻までの第１平均波形として、第２波形推定部２３に出力する。現時刻までの第１平均波形が第２波形推定部２３の入力波形となる。また、平均波形算出部２１３は、期間Ｔ１の第１平均波形ｙｔを振動成分算出部２１４に出力する。
【００４９】
振動成分算出部２１４は、平均波形算出部２１３から取得した第１平均波形ｙｔに基づき短周期Ｔｓの振動成分を算出する処理部である。すなわち、振動成分算出部２１４は、第１平均波形ｙｔに基づいて、期間Ｔ１における短周期Ｔｓの振動成分を算出する算出部の一例である。図９は、短周期Ｔｓの振動成分の算出を説明する図である。図９に示すように、振動成分算出部２１４は、入力波形ｐｔから平均波形算出部２１３が推定した第１平均波形ｙｔを減算して期間Ｔ１の振動成分を算出する。具体的に、振動成分算出部２１４は、式（９）に基づき、期間Ｔ１の振動成分ｒｔを算出する。
【００５０】
【数８】
2018184092.tif 000010
【００５１】
振動成分算出部２１４は、算出した振動成分ｒｔを危険判定部２５および波形パワー算出部２１５に出力する。
【００５２】
波形パワー算出部２１５は、振動成分ｒｔから波形パワーを算出する処理部である。波形パワー算出部２１５は、例えば現時刻から１波長前までの振動成分ｒｔのＲＭＳ値を波形パワーとして算出する（図９参照）。波形パワー算出部２１５は算出した波形パワーを危険判定部２５に出力する。
【００５３】
図１０および図１１は、第１平均波形推定のシミュレーション結果の一例を示す図である。図１０は、過去波形算出部２１１が算出した第１平均波形ｑｔを用いて期間Ｔ１の第１平均波形ｙｔを推定した場合の推定結果を示す。図１１は、過去波形算出部２１１が算出した第１平均波形ｑｔおよび入力波形ｐｔを用いて期間Ｔ１の第１平均波形ｙｔを推定した場合の推定結果を示す。図１０（１）および図１１（１）は推定した第１平均波形ｙｔと実際の第１平均波形ｙ１ｔとを示す。また、図１０（２）および図１１（２）は算出した第１平均波形ｙｔに基づいて算出した振動成分ｒｔと実際の振動成分ｒ１ｔとを示す。
【００５４】
図１０（１）に示すように、第１平均波形ｑｔを用いて期間Ｔ１の第１平均波形ｙｔを推定した場合、推定した第１平均波形ｙｔと実際の第１平均波形ｙ１ｔとが異なっている。そのため、図１０（２）に示すように、算出した振動成分ｒｔと実際の振動成分ｒ１ｔとが異なっている。これは、期間Ｔ１での振動成分ｒｔの傾向が期間Ｔ１以前の振動成分ｒｔと異なるためである。このように、振動成分ｒｔの傾向が変わると、算出した振動成分ｒｔに含まれる誤差が大きくなる。
【００５５】
一方、図１１（１）に示すように、第１平均波形ｑｔおよび入力波形ｐｔを用いて期間Ｔ１の第１平均波形ｙｔを推定した場合、第１波形推定部２２は、実際の第１平均波形ｙ１ｔとほぼ同じ波形の第１平均波形ｙｔを推定することができる。そのため、図１１（２）に示すように、第１波形推定部２２は、実際の振動成分ｒ１ｔとほぼ同じ波形の振動成分ｒｔを算出することができる。このように、第１平均波形ｑｔを用いて振動成分ｒｔを算出することで期間Ｔ１以前の振動成分ｒｔの傾向に応じた振動成分ｒｔを算出することができ、振動成分ｒｔに含まれる誤差を小さくすることができる。
【００５６】
なお、第２、第３波形推定部２３、２４も第１波形推定部２２と同様に、過去波形算出部２１１、パラメータ推定部２１２、平均波形算出部２１３、振動成分算出部２１４および波形パワー算出部２１５を有する。例えば、第２波形推定部２３の平均波形算出部２１３は、現時刻から中周期の半波長前までの第１平均波形ｙｔを入力波形として、第２期間における中周期の第２平均波形を推定する第２推定部として動作する。また、例えば第２波形推定部２３の振動成分算出部２１４は、第１平均波形ｙｔから第２平均波形を減算して中周期の振動成分を算出する第２算出部として動作する。
【００５７】
図１２は、危険判定部２５の機能構成を示す機能ブロック図である。危険判定部２５は、例えば各周期における振動成分の特徴量に基づき、車両の運転に関する危険度を判定する判定部の一例である。危険判定部２５は、操舵評価部２５１、揺動評価部２５２、危険度評価部２５３および通知部２５４を有する。
【００５８】
操舵評価部２５１は、運転者の運転操舵に関し、急操舵を起こした状態に基づく運転操作を評価する処理部である。操舵評価部２５１は、短周期Ｔｓの波形パワーを第１波形推定部２２から取得する。また、操舵評価部２５１は、中周期の振動成分の波形値（振幅値）を取得する。そして、操舵評価部２５１は、取得した波形パワーおよび波形値に基づき、急操舵スコア値を算出し、算出した急操舵スコア値を危険度評価部２５３に出力する。
【００５９】
具体的には、操舵評価部２５１は、まず評価関数ｆ１（ｘ）を用いて、急操舵波形パワースコア値ｓ１を算出する。例えば、短周期Ｔｓの波形パワーをｐ１とすると、操舵評価部２５１はｓ１＝ｆ１（ｐ１）を算出する。なお、評価関数ｆ１（ｘ）は、式（１０）に示すように波形パワーをスコア値に線形変換する、あるいは式（１１）に示すように波形パワーをスコア値に非線形変換する関数である。なお、式（１０）のα１およびβ１は、定数であり、スコア値に変換する波形パワーの種別ごとに設定される。
【００６０】
【数９】
2018184092.tif 000011
【００６１】
また、操舵評価部２５１は、中周期の振動成分の波形値（通常操舵揺動波形値）ｗ２を用いて通常操舵揺動率ｔｖを算出する。具体的には、式（１２）に示すように、波形値ｗ２の絶対値を閾値距離ｂａｓｅで除算した値および下限値のうち小さい値を通常操舵揺動率ｔｖとする。なお、閾値距離ｂａｓｅは、予め設定される値であり、例えば２００ｍｍである。また、式（１２）では下限値を「１」としている。続いて、操舵評価部２５１は、急操舵波形パワースコア値ｓ１と通常操舵揺動率ｔｖとを乗算して、急操舵スコア値ｓｃ１＝ｓ１＊ｔｖを算出する。
【００６２】
【数１０】
2018184092.tif 000012
【００６３】
揺動評価部２５２は、運転者の運転操舵に関し、走行軌跡の全体揺動の大きさを評価する処理部である。揺動評価部２５２は、第２、第３波形推定部２３、２４から取得した中周期および長周期の波形パワーに基づいて全体揺動スコア値を算出し、算出した全体揺動スコア値を危険度評価部２５３に出力する。
【００６４】
具体的には、揺動評価部２５２は、まず評価関数ｆ１（ｘ）を用いて、通常操舵波形パワースコア値Ｓ２および路内揺動波形パワースコア値Ｓ３を算出する。例えば、中周期の波形パワーをｐ２、長周期の波形パワーをｐ３とすると、揺動評価部２５２はＳ２＝ｆ１（ｐ２）およびＳ３＝ｆ１（ｐ３）を算出する。続いて、揺動評価部２５２は、式（１３）を用いてパワースコア高騰率ｇを算出する。
【００６５】
【数１１】
2018184092.tif 000013
【００６６】
そして、揺動評価部２５２は、通常操舵波形パワースコア値Ｓ２と路内揺動波形パワースコア値Ｓ３との加算結果にパワースコア高騰率ｇを乗算して全体揺動スコア値ｓｃ２＝（Ｓ２＋Ｓ３）＊ｇを算出する。
【００６７】
危険度評価部２５３は、操舵評価部２５１の評価結果と揺動評価部２５２の評価結果とに基づき、運転者の運転の危険度を評価する処理部である。具体的には、まず、走行軌跡から車両が走行している横変位ｍ１を取得する。続いて、危険度評価部２５３は、評価関数ｆ２（ｘ）を用いて、横変位危険スコア値ｃｒを算出する。例えば、危険度評価部２５３は、ｃｒ＝ｆ２（ｍ１）を算出する。なお、評価関数ｆ２（ｘ）は、非線形変換を行う場合は、式（１４）に示すように評価関数ｆ１（ｘ）と同じ関数であり、線形変換を行う場合は、式（１５）に示す一次関数である。なお、式（１５）のα２およびβ２は定数である。
【００６８】
【数１２】
2018184092.tif 000014
【００６９】
続いて、危険度評価部２５３は、急操舵スコア値ｓｃ１と全体揺動スコア値ｓｃ２との加算結果に横変位危険スコア値ｃｒを乗算して危険度スコア値ｓｃ３＝（ｓｃ１＋ｓｃ２）＊ｃｒを算出する。
【００７０】
そして、危険度評価部２５３は、算出した危険度スコア値ｓｃ３を通知部２５４に出力する。また、対象となった走行軌跡の時間と危険度スコア値ｓｃ３とを対応付けて、評価値ＤＢ１５に格納する。なお、危険度スコア値ｓｃ３に加え、急操舵スコア値ｓｃ１および全体揺動スコア値ｓｃ２を評価値ＤＢ１５に格納するようにしてもよい。
【００７１】
通知部２５４は、一定値を超える危険度スコア値ｓｃ３を危険度が高いと判定して、自動運転制御装置２に通知する処理部である。例えば、通知部２５４は、危険度評価部２５３から取得した危険度スコア値ｓｃ３が一定値を超える場合、危険ふらつきが発生したことを自動運転制御装置２に通知する。自動運転制御装置２は、通知部２５４からの通知を受けると、自動運転モードに切り替える。また、通知部２５４は、危険ふらつきが発生していることを示すメッセージをディスプレイ等に表示したり、音声で通知したりすることで、運転者に危険を通知するようにしてもよい。
【００７２】
［処理の流れ］
図１３は、処理の流れを示すフローチャートである。図１３に示すように、第１波形推定部２２が短周期の波形推定処理を実行し（Ｓ１０１）、第２波形推定部２３が中周期の波形推定処理を実行する（Ｓ１０２）。続いて、第３波形推定部２４が長周期の波形推定処理を実行する（Ｓ１０３）。その後、危険判定部２５が、危険判定処理を実行し（Ｓ１０４）、処理を終了する。
【００７３】
図１４は、波形推定処理の流れを示すフローチャートである。ここでは、第１波形推定部２２が行う短周期の波形推定処理について説明する。なお、第２、第３波形推定部２３、２４が実行する波形推定処理も短周期の波形推定処理と同様であるため、説明を省略する。
【００７４】
図１４に示すように、第１波形推定部２２は過去波形算出を実行する。具体的には、過去波形算出部２１１は、波長を算出し（Ｓ２０１）、第１平均波形を算出し（Ｓ２０２）、振動成分を算出する（Ｓ２０３）。続いて、過去波形算出部２１１は、波形パワーを算出し（Ｓ２０４）、算出した各波形等を波形ＤＢ１４に格納する（Ｓ２０５）。
【００７５】
続いて、第１波形推定部２２は現在波形算出を実行する。具体的には、パラメータ推定部２１２は、期間Ｔ１における第１平均波形のパラメータを推定する（Ｓ２０６）。続いて、平均波形算出部２１３は、第１平均波形を算出する（Ｓ２０７）。振動成分算出部２１４は、Ｓ２０７で算出した第１平均波形から振動成分を算出する（Ｓ２０８）。そして、波形パワー算出部２１５がＳ２０８で算出した振動成分の波形パワーを算出する（Ｓ２０９）。その後、第１波形推定部２２は、各波形等を後段の処理部に出力し（Ｓ２１０）、処理を終了する。
【００７６】
図１５は、危険判定処理の流れを示すフローチャートである。図１５に示すように、操舵評価部２５１は急操舵の評価値（急操舵スコア値ｓｃ１）を算出する（Ｓ３０１）。続いて、揺動評価部２５２は全体揺動の評価値（全体揺動スコア値ｓｃ２）を算出する（Ｓ３０２）。その後、危険度評価部２５３は、Ｓ３０１で算出された急操舵の評価値とＳ３０２で算出された全体揺動の評価値とに基づいて危険度の評価値（危険度スコア値ｓｃ３）を算出する（Ｓ３０３）。そして、通知部２５４は、危険度の評価値が閾値以上であれば（Ｓ３０４：Ｙｅｓ）、評価結果を自動運転制御装置２に出力し（Ｓ３０５）、評価値が閾値未満であれば（Ｓ３０４：Ｎｏ）、処理を終了する。
【００７７】
［効果］
上記実施例によれば、走行軌跡に応じた入力波形に応じて、現時刻から振動成分の半波長前までの第１平均波形を推定し、振動成分を算出する。平均波形の推定では、入力波形に加え、現時刻から半波長前より以前の第１平均波形を用いて平均波形を推定することで、現時刻までの第１平均波形をより高精度に推定可能である。これにより、走行軌跡から各周期の振動成分をリアルタイムで分離することができる。
【００７８】
また、リアルタイムで分離した各周期の振動成分に基づき、全体揺動および急操舵を評価し、危険ふらつきを検出する。これにより、集中力低下で車線内維持が徐々に困難（潜在的危険）となり生じる、逸脱に至る前の危険なふらつきを、リアルタイムで精度よく検出できる。すなわち、走行軌跡から各周期の振動成分をリアルタイムで分離することができるため、危険な車両のふらつきの検出時間を短縮することができる。
【００７９】
従来技術では、ある程度長い時間データを蓄積して運転者の覚醒状態を判定する装置が存在する。これに対して、車載装置１０は、リアルタイムで危険なふらつきの検出が可能であり、検出時に自動運転制御装置２に通知することで、より早い段階から車線逸脱を予防することができる。
【実施例２】
【００８０】
さて、これまで本発明の実施例について説明したが、本発明は上述した実施例以外にも、種々の異なる形態にて実施されてよいものである。
【００８１】
［判定手法］
上記実施例で説明した危険判定部２５による評価値の算出は、具体的に記載した算出例に限定されるものではない。例えば、危険判定部２５が短周期の振動成分から路内揺動の位相等も用いて評価値を算出するようにしてもよい。
【００８２】
例えば図１２に示す危険判定部２５の操舵評価部２５１が、車両の横変位、路内揺動の位相および短周期Ｔｓの振動成分ｒｔの波形パワーを用いて運転者の急操舵の検出および評価を実行するようにしてもよい。
【００８３】
具体的には、操舵評価部２５１は、走行軌跡から車両が走行している横変位を取得する。また、操舵評価部２５１は、短周期の振動成分から路内揺動の位相を算出する。つまり、操舵評価部２５１は、長期的視野で運転者のハンドル操作（操舵）のタイミングを示す路内揺動の位相（図２のｄの変化）を算出する。ここで算出される位相は、ハンドル操作を行ったタイミングを示すものであることから、白線の実線や破線で示す道路端（揺動端）ほど危険度が高い操舵を示す。なお、操舵評価部２５１は、様々な公知の手法を用いて位相を算出することができる。
【００８４】
図１６は、操舵評価部２５１の処理を説明する図である。図１６に示すように、操舵評価部２５１は、横変位にフィルタをかけて正規化した値を算出する。同様に、操舵評価部２５１は、短周期の波形パワー（急操舵パワー）に対してパワーフィルタをかけて正規化し、路内揺動の位相に対し位相フィルタをかけて正規化する。そして、操舵評価部２５１は、正規化された各値を加算して、急操舵のスコア値を算出し、危険度評価部２５３に出力する。
【００８５】
ここで、操舵評価部２５１が用いる各フィルタを説明する。図１７は、急操舵の評価に使用するフィルタを説明する図である。図１７に示すフィルタは、サチュレートするフィルタ、すなわち大小の値でともに飽和するフィルタである。図１７の（１）に示すフィルタは、横変位量に対するフィルタであり、左右のうちより白線に近い側の距離により応答値を得るフィルタであり、白線との距離が短い場合に高い応答値を示すフィルタである。例えば、白線との距離（白線距離）が近い距離（ＴＨ１１）より長くなると、応答値が徐々に小さく変化し、閾値ＴＨ１２を超えると一定値に収束する。
【００８６】
図１７の（２）に示すフィルタは、路内揺動の位相に対するフィルタであり、路内揺動の位相θからｓｉｎθを算出し、より位相が波端の場合に高い応答値を示すフィルタである。例えば、波端の閾値ＴＨ２１を超えると応答値が徐々に大きくなり、閾値ＴＨ２２を超えると一定値に収束する。
【００８７】
図１７の（３）に示すフィルタは、急操舵パワーに対するフィルタであり、閾値ＴＨ３１を超えた場合に強く反応するフィルタである。例えば、閾値ＴＨ３１よりもパワーが大きい閾値ＴＨ３２の方が応答最大値に近い値が返される。なお、閾値ＴＨ１１〜ＴＨ３２の各閾値はそれぞれ個別の値で事前に定めておく。
【００８８】
そして、操舵評価部２５１は、フィルタ処理後の各値を積算することでスコア値を算出して急操舵を判定する。フィルタ処理後の各特徴量の積算は、すべてのフィルタ処理後の値が高い場合に高い評価値となる。図１８は、急操舵の判定を説明する図である。
【００８９】
図１８のａ５は走行軌跡の高周波成分、図１８のａ６は低周波成分を表しており、横変位量のフィルタ処理後の値が高くなるのは、図１８のＲｃの領域であり、路内揺動の位相のフィルタ処理後の値が高くなるのは図１８のＲｄの領域である。このため、急な操舵修正の評価結果がより明確に現れるのは、図１８のＲｅ１、Ｒｅ２の２領域である。この例では、Ｒｅ１、Ｒｅ２のうち、より高周波パワーが高いＲｅ２がより危険度が高いタイミングであると判定される。
【００９０】
つまり、操舵評価部２５１は、走行軌跡の中から、路内揺動の位相を取得することで、位相端（Ｒｄ）として４か所が特定できる。続いて、操舵評価部２５１は、特定した４か所のうち、車線内の横位置（ａ５）に基づいて車両が白線と近い位置となった地点Ｒｅ１とＲｅ２とを抽出する。その後、操舵評価部２５１は、地点Ｒｅ１とＲｅ２のうち、短周期の急操舵パワーが閾値以上である地点Ｒｅ２を特定する。この結果、操舵評価部２５１は、走行軌跡のうちＲｅ２で危険運転が発生したと評価する。このように、操舵評価部２５１は、路内揺動の位相と車線内の横位置から、逸脱を避けるための操舵修正タイミングに絞り込み、集中度低下に伴う修正操舵の急激性の評価を行う。
【００９１】
また、図１２に示す揺動評価部２５２が、長周期の振動成分の波形パワーと中周期の振動成分の波形パワーとを用いて、運転者の揺動を評価するようにしてもよい。つまり、揺動評価部２５２は、路内揺動の波形パワーと通常操舵の波形パワーとを加算して、危険度を示すスコア値を算出する。
【００９２】
図１９は、揺動評価部２５２の処理を説明する図である。図１９に示すように、揺動評価部２５２は、路内揺動の波形パワーにフィルタをかけて波形パワーを正規化し、同様に、通常操舵の波形パワーにフィルタをかけて波形パワーを正規化する。そして、揺動評価部２５２は、フィルタ処理後の各波形パワーを加算してスコア値を算出して、危険度評価部２５３に出力する。このようにすることで、揺動評価部２５２は、いずれか一方の波形パワーが高いときにスコア値が高くなることを除去し、両方の波形パワーが高いときにスコア値が高くなるように制御する。
【００９３】
図２０は、全体揺動の評価を説明する図である。図２０は、路内揺動の波形パワーと通常操舵の波形パワーの両パワーの時系列変化の例を示している。図２０の（１）は、路内揺動および通常操舵の２つの波形を示しており、ａ１１は、一方だけの波形パワーが高い状態であり、ｂ１１は、両方の波形パワーが高い状態を表している。
【００９４】
このような波形を単純に加算した状態を図２０の（５）に示す。図２０の（５）では、一方だけ非常に高いａ１５の状態も、両方とも少し高いｂ１５の状態も同様のスコア値となり、両パワーがともに高まっているかどうかの判別が難しい。つまり、このような処理は、癖がある運転者の通常操舵についても危険と判定するので、通常操舵とは異なる意識的な急操舵すなわち本当に危険な状態を特定することができない。
【００９５】
これに対して、揺動評価部２５２は、図２０の（１）に示す２つの波形に対して、（４）に示すフィルタをかける。このフィルタは、通常操舵と路内揺動でパラメータは異なるが、基本的な性質は同じである。すなわち低パワーの閾値ＴＨ１までは、最小値近傍値を応答し、高パワーの閾値ＴＨ２より最高値近傍値を応答し、閾値ＴＨ１、ＴＨ２間ではパワーの大きさに従い応答する性質を持つ。
【００９６】
そして、揺動評価部２５２が、図２０の（１）の波形に図２０の（４）のフィルタによってフィルタ処理することで、両方の波形パワーが高い状態（ｂ１１）の両パワーの応答値が高くなる（ｂ１２）。その後、揺動評価部２５２は、得られたフィルタ処理後の両パワーの単純加算を行うことで、図２０の（３）に示す波形を得ることができる。なお、適宜さらにパワーの比を重み係数にするなどしてもよい。図２０の（３）では、図２０の（２）で引き上げられたｂ１２が、加算の結果、ａ１３の領域に比べ高いスコアとなり、事前に定めた全体揺動の危険性基準のｃ１を超える。
【００９７】
このように、図２０の（４）に示すフィルタ処理後の波形パワーを評価することで、通常操舵と路内揺動の両波形パワーの大きさに加え、ともに高まっている状態を評価することが可能となり、ロバストに継続的な集中力低下状態を特定可能となる。なお、各フィルタの応答値を同範囲、例えば０〜１．０に設定することで、異なるパワーレンジの高まりを正規化して評価することが可能である。
【００９８】
さらに、図１２に示す危険度評価部２５３が、操舵評価部２５１により算出されたスコア値と揺動評価部２５２により算出されたスコア値とを加算した評価値を算出するようにしてもよい。
【００９９】
なお、揺動評価と操舵評価はいずれを先に実行してもよく、処理する順番は任意に変更することができる。また、揺動評価と操舵評価のうちいずれかの評価だけを実行してもよく、両方を実行してもよい。
【０１００】
［数値］
上記実施例で説明した各数値や閾値はあくまで例示であり、具体的に記載した値に限定されるものではない。
【０１０１】
［システム］
また、図３に示した各装置の各構成は、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、任意の単位で分散または統合して構成することができる。例えば、第１〜第３波形推定部２２〜２４を統合することができる。さらに、各装置にて行なわれる各処理機能は、その全部または任意の一部が、ＣＰＵ（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）などのコンピュータで読み取り可能な記録媒体に記録され、コンピュータによって記録媒体から読み出されることによって実行することができる。
【符号の説明】
【０１０８】
２  自動運転制御装置
１０  車載装置
２１  軌跡生成部
２２  第１波形推定部
２３  第２波形推定部
２４  第３波形推定部
２５  危険判定部
２１１  過去波形算出部
２１２  パラメータ推定部
２１３  平均波形算出部
２１４  振動成分算出部
２１５  波形パワー算出部
【図１】
2018184092.tif 000015
【図２】
2018184092.tif 000016
【図３】
2018184092.tif 000017
【図４】
2018184092.tif 000018
【図５】
2018184092.tif 000019
【図６】
2018184092.tif 000020
【図７】
2018184092.tif 000021
【図８】
2018184092.tif 000022
【図９】
2018184092.tif 000023
【図１０】
2018184092.tif 000024
【図１１】
2018184092.tif 000025
【図１２】
2018184092.tif 000026
【図１３】
2018184092.tif 000027
【図１４】
2018184092.tif 000028
【図１５】
2018184092.tif 000029
【図１６】
2018184092.tif 000030
【図１７】
2018184092.tif 000031
【図１８】
2018184092.tif 000032
【図１９】
2018184092.tif 000033
【図２０】
2018184092.tif 000034
【図２１】
2018184092.tif 000035