(57)【特許請求の範囲】
【請求項１】
回転機械の運転状態を診断する診断装置であって、
時系列で前記回転機械を計測した計測値から、自己回帰診断モデルの特性を規定するモデルパラメータを学習するための第１運転データを生成し、かつ、前記回転機械の運転状態を表す特徴量を計算するための第２運転データを生成するデータ前処理部と、
前記第１運転データを用いて、前記モデルパラメータを学習するモデル学習部と、
前記第２運転データと、学習された前記モデルパラメータと、前記自己回帰診断モデルとを用いて、前記特徴量を計算する特徴量計算部と、
前記特徴量に基づき、前記回転機械の運転状態を診断する状態診断部と、
を備え、
前記データ前処理部は、前記計測値が含まれる前記回転機械の起動から停止までの作動期間をバッチとして特定するバッチ情報を有する前記第１運転データを生成し、
前記モデル学習部は、前記バッチ情報を用いて前記計測値の時系列の連続性を確保しつつ、前記モデルパラメータを学習する、
診断装置。
【請求項２】
前記モデル学習部は、前記モデルパラメータを学習するための学習行列を作成する行列作成部を含み、
前記行列作成部は、前記学習行列の同一行内には同一の前記バッチ内の前記計測値のみが含まれるように、前記学習行列を作成する、
請求項１に記載の診断装置。
【請求項３】
前記行列作成部は、複数の前記バッチ内の前記計測値が含まれるように、前記学習行列を作成する、
請求項２に記載の診断装置。
【請求項４】
前記回転機械は、空気圧縮機を含む、
請求項１〜３のいずれか１項に記載の診断装置。
【請求項５】
回転機械の運転状態を診断する診断装置に用いられる診断方法であって、
時系列で前記回転機械を計測した計測値から、自己回帰診断モデルの特性を規定するモデルパラメータを学習するための第１運転データを生成し、かつ、前記回転機械の運転状態を表す特徴量を計算するための第２運転データを生成するデータ前処理工程と、
前記第１運転データを用いて、前記モデルパラメータを学習するモデル学習工程と、
前記第２運転データと、学習された前記モデルパラメータと、前記自己回帰診断モデルとを用いて、前記特徴量を計算する特徴量計算工程と、
前記特徴量に基づき、前記回転機械の運転状態を診断する状態診断工程と、
を備え、
前記データ前処理工程では、前記計測値が含まれる前記回転機械の起動から停止までの作動期間をバッチとして特定するバッチ情報を有する前記第１運転データが生成され、
前記モデル学習工程では、前記バッチ情報を用いて前記計測値の時系列の連続性を確保しつつ、前記モデルパラメータが学習される、
診断方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、回転機械の運転状態を診断する技術に関するものである。
【背景技術】
【０００２】
従来、省人化、ヒューマンエラーの削減などを目的に、機械製品又はプラントの稼働状態、生産ラインの品質特性などに関する正常又は異常を自動的に診断する様々な手法が考案されている。これらの手法では、基本的に、以下の手順に従って診断が行われる。
【０００３】
まず、診断対象が過去に正常に稼働していた時の時系列データに対してデータマイニング、機械学習、統計処理等の分析手法を適用して、診断に用いる特徴量を計算する特性モデルを学習する。そして、新たに観測される時系列データに対して該特性モデルが計算した特徴量の、正常時の特徴量の確率分布からの外れ度を計算し、その外れ度が予め決定された閾値を超えた場合に異常と判断する。
【０００４】
上述の特性モデルの学習法としては、重回帰分析や主成分分析等の統計学習法なども用いられるが、時系列データを扱う点との親和性の高さから、特許文献１に記載の自己回帰（ＡＲ）モデル、ベクトル自己回帰（ＶＡＲ）モデル等の時系列データ分析手法がしばしば用いられる。
【０００５】
特許文献１に記載の技術は、時系列データに対する自己回帰モデル学習装置並びにそれを用いた外れ値および変化点の検出装置に関する。特に、特許文献１に記載の技術は、データ解析技術およびデータマイニング技術に関わり、順次入力される離散値変量と連続値変量との両者または一方で記述されたデータに対して、その外れ値スコアおよび変化点スコアを計算して、外れ値および変化点を精度よく検出する検出装置に関する。
【０００６】
特許文献１に記載の技術による時系列データに対する外れ値および変化点の検出装置には、実数ベクトル値のデータ列を順次読み込みながら該データ列の発生する確率分布を、自己回帰モデルを用いて学習する装置として、次のような自己回帰モデル学習装置を用いることができる。すなわち、自己回帰モデルの十分統計量を、新たに読み込んだデータを用いて過去のデータを忘却しながら更新するデータ更新装置と、該データ更新装置が更新した十分統計量を読み込み、それを用いて自己回帰モデルのパラメータを計算するパラメータ計算装置とを含むものである。このように、過去のデータを忘却しながら更新する処理を行うので、時系列データを処理するのに適すると共に、その精度を上げることができる。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開２００４−５４３７０号公報
【発明の概要】
【発明が解決しようとする課題】
【０００８】
上記特許文献１に記載の技術では、自己回帰モデルを用いた時系列モデリングによって、前記のメカニズムに従った正常又は異常の診断が自動的に実行される。その結果、診断処理の省人化、ヒューマンエラーの削減が期待できる。
【０００９】
しかしながら、本診断手法の適用対象を、空気圧縮機、ガスタービン、蒸気タービンなどの回転機械に限定した場合については、上記特許文献１に記載の技術では、回転機械の運用の特徴が十分に考慮されていない。このため、回転機械の運転状態を精度良く診断できるように改善することが求められている。
【００１０】
本発明は、上記問題を解決するもので、回転機械の運転状態を精度良く診断することができる診断装置及び診断方法を提供することを目的とする。
【課題を解決するための手段】
【００１１】
本発明の第１態様は、
回転機械の運転状態を診断する診断装置であって、
時系列で前記回転機械を計測した計測値から、自己回帰診断モデルの特性を規定するモデルパラメータを学習するための第１運転データを生成し、かつ、前記回転機械の運転状態を表す特徴量を計算するための第２運転データを生成するデータ前処理部と、
前記第１運転データを用いて、前記モデルパラメータを学習するモデル学習部と、
前記第２運転データと、学習された前記モデルパラメータと、前記自己回帰診断モデルとを用いて、前記特徴量を計算する特徴量計算部と、
前記特徴量に基づき、前記回転機械の運転状態を診断する状態診断部と、
を備え、
前記データ前処理部は、前記計測値が含まれる前記回転機械の起動から停止までの作動期間をバッチとして特定するバッチ情報を有する前記第１運転データを生成し、
前記モデル学習部は、前記バッチ情報を用いて前記計測値の時系列の連続性を確保しつつ、前記モデルパラメータを学習するものである。
【００１２】
本発明の第２態様は、
回転機械の運転状態を診断する診断装置に用いられる診断方法であって、
時系列で前記回転機械を計測した計測値から、自己回帰診断モデルの特性を規定するモデルパラメータを学習するための第１運転データを生成し、かつ、前記回転機械の運転状態を表す特徴量を計算するための第２運転データを生成するデータ前処理工程と、
前記第１運転データを用いて、前記モデルパラメータを学習するモデル学習工程と、
前記第２運転データと、学習された前記モデルパラメータと、前記自己回帰診断モデルとを用いて、前記特徴量を計算する特徴量計算工程と、
前記特徴量に基づき、前記回転機械の運転状態を診断する状態診断工程と、
を備え、
前記データ前処理工程では、前記計測値が含まれる前記回転機械の起動から停止までの作動期間をバッチとして特定するバッチ情報を有する前記第１運転データが生成され、
前記モデル学習工程では、前記バッチ情報を用いて前記計測値の時系列の連続性を確保しつつ、前記モデルパラメータが学習されるものである。
【００１３】
この第１及び第２態様では、時系列で回転機械を計測した計測値から、自己回帰診断モデルの特性を規定するモデルパラメータを学習するための第１運転データが生成される。また、時系列で回転機械を計測した計測値から、回転機械の運転状態を表す特徴量を計算するための第２運転データが生成される。第１運転データを用いて、モデルパラメータが学習される。第２運転データと、学習されたモデルパラメータと、自己回帰診断モデルとを用いて、特徴量が計算される。特徴量に基づき、回転機械の運転状態が診断される。
【００１４】
また、第１及び第２態様では、計測値が含まれる回転機械の起動から停止までの作動期間をバッチとして特定するバッチ情報を有する第１運転データが生成される。モデルパラメータは、バッチ情報を用いて、計測値の時系列の連続性を確保しつつ、学習される。したがって、起動から停止までの作動期間が短期間となるような回転機械の運転状態を診断するための診断モデルを、時系列の自己回帰手法を用いて構築する場合であっても、学習データの時系列連続性を確保することができる。このため、十分なモデル精度を得ることができる。その結果、回転機械の運転状態を精度良く診断することができる。
【００１５】
上記第１態様において、例えば、前記モデル学習部は、前記モデルパラメータを学習するための学習行列を作成する行列作成部を含んでもよい。前記行列作成部は、前記学習行列の同一行内には同一の前記バッチ内の前記計測値のみが含まれるように、前記学習行列を作成してもよい。
【００１６】
本態様では、モデルパラメータを学習するための学習行列の同一行内には同一のバッチ内の計測値が含まれるように、学習行列が作成される。したがって、本態様によれば、学習データの時系列の連続性を確保することができる。このため、十分なモデル精度を得ることができる。
【００１７】
上記第１態様において、例えば、前記行列作成部は、複数の前記バッチ内の前記計測値が含まれるように、前記学習行列を作成してもよい。
【００１８】
本態様では、複数のバッチ内の計測値が含まれるように、学習行列が作成される。したがって、本態様によれば、学習汎化のためのデータ量を十分に確保して、学習を実行することができる。このため、十分なモデル精度を得ることができる。
【００１９】
上記第１態様において、例えば、前記回転機械は、空気圧縮機を含んでもよい。
【発明の効果】
【００２０】
本発明によれば、計測値が含まれる回転機械の起動から停止までの作動期間をバッチとして特定するバッチ情報を有する第１運転データが生成され、モデルパラメータは、バッチ情報を用いて、計測値の時系列の連続性を確保しつつ、学習されることから、起動から停止までの作動期間が短期間となるような回転機械の運転状態を診断するための診断モデルを、時系列の自己回帰手法を用いて構築する場合であっても、学習データの時系列連続性を確保することができるため、十分なモデル精度を得ることができるので、回転機械の運転状態を精度良く診断することができる。
【図面の簡単な説明】
【００２１】
【図１】本発明の一実施形態の診断装置の構成を概略的に示すブロック図である。
【図２】診断装置に含まれるモデル学習部の構成を概略的に示すブロック図である。
【図３】前処理運転データの一例を概略的に示す図である。
【図４】前処理運転データの一例を概略的に示す図である。
【図５】図１に示される診断装置の動作を概略的に示すフローチャートである。
【図６】図５のフローチャートの診断モデルの学習処理のサブルーチンを概略的に示すフローチャートである。
【図７】学習行列の更新を説明する図である。
【図８】図７と異なる学習行列の更新を説明する図である。
【図９】図６のループによる学習行列の更新を説明する図である。
【図１０】本実施形態の学習行列の行とバッチとの関係を説明する図である。
【図１１】比較例の学習行列の行とバッチとの関係を説明する図である。
【発明を実施するための形態】
【００２２】
（本発明の基礎となった知見）
空気圧縮機、ガスタービン、蒸気タービンなどの回転機械は、使用方針にもよるが、一般的に、起動から停止までの作動時間が短期間（数時間〜数日間）で運用されることが多い。このため、回転機械の運転状態を診断するためには、短期間で起動及び停止が繰り返される運用を前提条件として考慮する必要がある。このような前提条件を考慮しようとすると、上記特許文献１に記載の技術では、以下の問題が生じる虞がある。
【００２３】
すなわち、短期間での起動及び停止の繰り返しにより、時系列データの連続性が保証される期間におけるデータ量が少なくなる。上記特許文献１に記載の技術では、１回のモデル学習に使用できるデータは連続するデータが１セットのみのため、学習に使用するデータが不足することから、所望の精度の診断結果が得られない可能性がある。
【００２４】
以上の考察により、本発明者は、回転機械の運転状態を精度良く診断することができる発明を想到するに至った。
【００２５】
すなわち、本開示の診断装置および診断方法では、短期間で起動及び停止が繰り返される運用を前提条件として有する回転機械の運転状態を診断する診断モデルを学習する際に、短期間の連続の時系列データが複数セット集約される。これによって、アルゴリズム上、時系列データの連続性を損なうことなく、診断モデルの学習を可能にしている。
【００２６】
（実施の形態）
以下、本発明の一実施の形態について、図面を参照しながら説明する。なお、各図面において、同じ構成要素については同じ符号が用いられている。
【００２７】
（構成）
図１は、本発明の一実施形態の診断装置１００の構成を概略的に示すブロック図である。図２は、診断装置１００に含まれるモデル学習部１０４の構成を概略的に示すブロック図である。
【００２８】
診断装置１００は、空気圧縮機１（回転機械の一例）の運転状態を診断する。図１に示されるように、診断装置１００は、空気圧縮機１からオンラインで信号を受信する。
【００２９】
診断装置１００は、例えばコンピュータであり、中央演算処理装置（ＣＰＵ）、プログラム及びデータ等を保存するメモリ、周辺回路、表示部１２０、及び入力部１３０を含む。メモリは、半導体メモリ、ハードディスク等で構成される。メモリは、運転データベース（ＤＢ）１０３、モデル情報ＤＢ１０５、診断条件ＤＢ１０７、出力ログＤＢ１０９を含む。
【００３０】
メモリに保存されたプログラムを実行することにより、ＣＰＵは、データ前処理部１０２、モデル学習部１０４、特徴量計算部１０６、状態診断部１０８として機能し、学習行列作成部１０４ａ（図２）、パラメータ計算部１０４ｂ（図２）として機能する。周辺回路は、データ入力インターフェース（Ｉ／Ｆ）１０１と、データ出力Ｉ／Ｆ１１０とを含む。
【００３１】
表示部１２０は、例えば液晶表示パネルを含む。代替的に、表示部１２０は、有機ＥＬ（Ｅｌｅｃｔｒｏ  Ｌｕｍｉｎｅｓｃｅｎｃｅ）パネルを含んでもよい。入力部１３０は、キーボード及びマウスを含んでもよい。表示部１２０がタッチパネル式の場合には、入力部１３０は、表示部１２０のタッチパネルであってもよい。
【００３２】
図１に示される診断装置１００は、概略次のような動作を行う。まず、診断装置１００のデータ前処理部１０２は、空気圧縮機１からオンラインで入力された計測信号を診断モデルの学習に適切な形式のデータを生成し、運転ＤＢ１０３に保存する。本実施形態では例えば、診断モデルとして、自己回帰モデル又はベクトル自己回帰モデルを用いる。次に、診断モデルの構築が必要なタイミングと判断された場合は、モデル学習部１０４は、運転ＤＢ１０３に保存されたデータを基に、診断モデルの特性を規定するモデルパラメータを学習する。以上のようにして得られたモデルパラメータは、モデル情報ＤＢ１０５に保存される。
【００３３】
次に、学習した診断モデルを用いた空気圧縮機１の運転状態の診断に関して、診断装置１００は以下の動作を行う。まず、データ前処理部１０２は、診断に用いる特徴量の計算に適切な形式のデータを生成する。このデータと、モデル情報ＤＢ１０５に保存されたモデルパラメータとを用いて、特徴量計算部１０６は、診断用の特徴量を計算する。次に、状態診断部１０８は、計算された特徴量を用いて、空気圧縮機１の運転状態が正常であるか異常であるかを判定する。判定結果は、表示部１２０に表示される。
【００３４】
以下、診断装置１００の具体的な構成が説明される。データ入力Ｉ／Ｆ１０１は、空気圧縮機１から出力される空気圧縮機１の各部の計測信号１０を予め定められた通信周期（例えば１秒）で受信する。
【００３５】
計測信号１０は、本実施形態では例えば、空気圧縮機１への吸込み空気圧、空気圧縮機１からの吐出空気圧、空気圧縮機１に潤滑用油を供給するための給油フィルタの前後の差圧、空気圧縮機１から流出する使用済み潤滑用油の不純物を分離する分離フィルタの前後の差圧、給油温度、油回収器の出口温度、空気圧縮機１の主電動機の軸受温度、空気圧縮機１の主電動機の巻線温度、空気圧縮機１のケーシング振動の振幅、温度によって変動する空気圧縮機１のロータ軸の軸方向長さの１０種類の信号を含む。
【００３６】
データ入力Ｉ／Ｆ１０１は、１０種類の計測信号１０を例えばパラレルデータとして受信する。データ入力Ｉ／Ｆ１０１は、パラレルデータとして受信した１０種類の計測信号１０を、診断装置１００の内部で処理可能な例えばシリアルデータの計測信号１１に変換する。データ入力Ｉ／Ｆ１０１は、計測信号１１をデータ前処理部１０２に出力する。
【００３７】
データ前処理部１０２は、入力された計測信号１１を診断モデルの学習に適切な形式に変換し、前処理運転データ１２（第１運転データの一例）を生成する。データ前処理部１０２は、生成した前処理運転データ１２を運転ＤＢ１０３に出力する。
【００３８】
図３は、前処理運転データ１２の一例を概略的に示す図である。前処理運転データ１２は、図３に示されるように、日付時刻１２０１、バッチＩＤ１２０２、時系列データＩＤ１２０３、第１〜第１０信号の計測データ１２０４を含む。
【００３９】
日付時刻１２０１は、計測データ１２０４の各レコードが計測された日付および時刻情報を表す。バッチＩＤ１２０２は、計測データ１２０４の時系列連続性を基準に分類されるバッチを特定する識別子である。本実施形態では、空気圧縮機１における起動から停止までの一連の連続運転中に得られた計測信号がひとつのバッチと定義される。バッチＩＤ１２０２として、時系列順に整数の添字ｉを用いて識別子が管理される。
【００４０】
時系列データＩＤ１２０３は、予め定められた通信周期毎に入力される計測信号１１に時系列順に付与される識別子である。時系列データＩＤ１２０３として、整数の添字ｊを用いて識別子が管理される。計測データ１２０４は、計測信号１０に含まれる計測値を表す。本実施形態では、上述のように１０種類の計測信号１０が空気圧縮機１から入力されるので、前処理運転データ１２は、第１〜第１０信号の１０個の計測データ１２０４を含む。なお、本実施形態では、バッチＩＤ１２０２の最大値はＩ、時系列データＩＤ１２０３の最大値（つまりバッチｉに属するデータの最大値）はＪｉとする。
【００４１】
図１に戻って、データ前処理部１０２は、入力された計測信号１１を空気圧縮機１の診断に用いる特徴量の計算に適切な形式に変換し、前処理運転データ１６（第２運転データの一例）を生成する。データ前処理部１０２は、生成した前処理運転データ１６を特徴量計算部１０６に出力する。
【００４２】
図４は、前処理運転データ１６の一例を概略的に示す図である。前処理運転データ１６は、図４に示されるように、日付時刻１２０１、計測データ１２０４を含む。すなわち、前処理運転データ１６は、前処理運転データ１２からバッチＩＤ１２０２及び時系列データＩＤ１２０３を除いた構成になっている。
【００４３】
図１に戻って、データ前処理部１０２は、データ入力Ｉ／Ｆ１０１が新たな計測信号１０を受信するたびに、前処理運転データ１２，１６を順次更新する。運転ＤＢ１０３は、データ前処理部１０２から出力された前処理運転データ１２を保存する。
【００４４】
モデル学習部１０４は、前処理運転データ１２と、モデル情報ＤＢ１０５に保存されているモデル構築情報１５とに基づき、モデルパラメータ１４を学習する。学習されたモデルパラメータ１４は、モデル情報ＤＢ１０５に保存される。
【００４５】
モデル構築情報１５は、自己回帰モデルを構成するために必要なモデル次数を含む。モデル次数は、参照する時系列データのステップ数である。つまり、モデル次数は、どれだけの過去のデータまでさかのぼって考慮するかを表す。
【００４６】
モデル構築情報１５は、出力変数を特定する情報を含む。モデル構築情報１５は、入力変数の組合せを特定する情報を含む。自己回帰モデルを用いる場合には、通常、入力変数の組合せは、出力変数を含む。自己回帰モデルとして、ベクトル自己回帰（ＶＡＲ）モデルを適用する場合には、モデル構築情報１５に含まれるモデルの入力変数の組合せは、複数の入力変数を含む。
【００４７】
出力変数及び入力変数は、それぞれ、計測信号１０として入力される１０種類の信号のなかから設定される。モデル構築情報１５は、出力変数として、例えば、空気圧縮機１からの吐出空気圧を含んでいてもよい。モデル情報ＤＢ１０５は、予め定められた出力変数と、予め定められた入力変数の組合せとを含むモデル構築情報１５を保存してもよい。
【００４８】
図２において、モデル学習部１０４の学習行列作成部１０４ａ（行列作成部の一例）は、運転ＤＢ１０３から前処理運転データ１２を読み取る。学習行列作成部１０４ａは、モデル情報ＤＢ１０５からモデル構築情報１５を読み取る。学習行列作成部１０４ａは、前処理運転データ１２に対し、モデル構築情報１５を適用し、データの時系列連続性が確保されるように、学習用の行列（以下、「学習行列」と称される）１３を作成する。学習行列作成部１０４ａは、作成した学習行列１３をパラメータ計算部１０４ｂに出力する。
【００４９】
パラメータ計算部１０４ｂは、モデル情報ＤＢ１０５からモデル構築情報１５を読み取る。パラメータ計算部１０４ｂは、入力された学習行列１３に対して、モデル構築情報１５を適用し、モデルパラメータ１４を計算する。パラメータ計算部１０４ｂは、計算したモデルパラメータ１４を、モデルパラメータ１４の学習結果として、モデル情報ＤＢ１０５に保存する。
【００５０】
診断条件ＤＢ１０７は、特徴量計算および状態診断の実行条件パラメータ１８を保存する。具体的な実行条件パラメータ１８は、状態診断に特徴量の閾値判定を用いる場合の閾値パラメータ、確率検定を用いる場合の有意水準パラメータ等を含む。
【００５１】
特徴量計算部１０６は、モデルパラメータ１４、前処理運転データ１６およびモデル構築情報１５を用いて、現時点での計測信号に対するモデル特性評価を実行する。特徴量計算部１０６は、そのモデル特性評価の結果を基に、診断用の実行条件パラメータ１８を用いて、診断用の指標である特徴量１９を計算する。
【００５２】
状態診断部１０８は、診断用の実行条件パラメータ１８と、特徴量１９とを用いて、空気圧縮機１の運転状態を診断し、診断結果２０を出力ログＤＢ１０９及びデータ出力Ｉ／Ｆ１１０に出力する。
【００５３】
出力ログＤＢ１０９は、特徴量計算部１０６により計算された特徴量１９と、状態診断部１０８の診断結果２０とを保存する。データ出力Ｉ／Ｆ１１０は、診断結果２０を含む、表示部１２０で表示される形式の表示データ２１を生成し、生成した表示データ２１を表示部１２０に出力する。表示部１２０は、表示データ２１を表示する。
【００５４】
（動作）
図５は、図１に示される診断装置１００の動作を概略的に示すフローチャートである。診断装置１００は、空気圧縮機１から計測信号１０を受信し、空気圧縮機１の運転状態を診断し、診断結果を表示部１２０に表示する。計測信号の通信タイミングおよび診断結果の更新タイミングは、それぞれ、予め定められた通信周期（例えば１秒）および診断周期（例えば２秒）に依存する。一般に、通信周期＜診断周期となるのが望ましい。したがって、診断装置１００の図５に示される動作は、通信周期毎に実行される。なお、通信周期、診断周期は、診断装置１００を含むプラント計装システムの通信速度、演算処理速度に依存して決定されてもよい。以下、図５の各ステップが説明される。
【００５５】
ステップＳ２０００において、診断装置１００のデータ入力Ｉ／Ｆ１０１は、空気圧縮機１から計測信号１０を受信する。データ入力Ｉ／Ｆ１０１は、計測信号１０を診断装置１００において処理可能な形式である計測信号１１に変換する。データ入力Ｉ／Ｆ１０１は、計測信号１１をデータ前処理部１０２に出力する。
【００５６】
ステップＳ２０１０（データ前処理工程の一例）において、データ前処理部１０２は、計測信号１１を診断モデルの学習に必要な形式に整形し、前処理運転データ１２を生成する。データ前処理部１０２は、生成した前処理運転データ１２を運転ＤＢ１０３に出力する。データ前処理部１０２は、計測信号１１を、診断に用いる特徴量の計算に必要な形式に整形し、前処理運転データ１６を生成する。データ前処理部１０２は、生成した前処理運転データ１６を特徴量計算部１０６に出力する。
【００５７】
ステップＳ２０２０において、運転ＤＢ１０３は、データ前処理部１０２から入力された前処理運転データ１２を保存する。
【００５８】
ステップＳ２０３０において、診断装置１００のＣＰＵは、診断モデルの学習を実行するか否かを、予め設定された実行条件を基に判定する。本実施形態では例えば、診断装置１００のＣＰＵは、一定周期（週単位、月単位など）の経過時に実行条件を満たすと判定して診断モデルの学習を実行する。なお、ユーザが入力部１３０を用いて所定の入力を行うと、診断装置１００のＣＰＵは、実行条件を満たすと判定して診断モデルの学習を実行してもよい。また、十分なモデル精度を獲得するために最低限必要なデータ量を閾値として設定しておいてもよい。この実施形態では、診断装置１００のＣＰＵは、運転ＤＢ１０３に蓄積されたデータ量が閾値に満たない場合には診断モデルの学習を実行しないと判定する。
【００５９】
判定の結果、診断モデルの学習を実行する場合には（ステップＳ２０３０でＹＥＳ）、処理はステップＳ２０４０に進む。一方、学習を実行しない場合には（ステップＳ２０３０でＮＯ）、処理はステップＳ２０６０に進む。
【００６０】
ステップＳ２０４０（モデル学習工程の一例）において、モデル学習部１０４は、運転ＤＢ１０３に保存されている前処理運転データ１２と、モデル情報ＤＢ１０５に保存されているモデル構築条件１５とを用いて、診断モデルを学習する。なお、この診断モデルの学習の詳細な処理内容については、図６のフローチャートを用いて次に説明される。
【００６１】
図６は、図５のフローチャートの診断モデルの学習処理（ステップＳ２０４０）のサブルーチンを概略的に示すフローチャートである。図７は、学習行列Ａの更新を説明する図である。図８は、学習行列Ｂの更新を説明する図である。
【００６２】
図６のステップＳ４０００において、モデル学習部１０４の学習行列作成部１０４ａは、運転ＤＢ１０３から前処理運転データ１２を読み込み、モデル情報ＤＢ１０５からモデル構築情報１５を読み込む。モデル構築情報１５は、上述のように、診断モデルを構成するために必要なモデル次数Ｎと、出力変数を特定する情報とを含む。また、モデル構築情報１５は、入力変数の組合せ情報Ｍを含む。なお、診断モデルとして、ベクトル自己回帰（ＶＡＲ）モデルを適用する場合には、入力変数の組合せ情報Ｍの大きさ｜Ｍ｜、つまり入力変数の個数は２以上の整数であり、ＶＡＲモデル以外の場合には、入力変数の組合せ情報Ｍの大きさ｜Ｍ｜は１である。
【００６３】
ステップＳ４０１０において、学習行列作成部１０４ａは、診断モデルの学習に使用する２種類の学習行列Ａ、Ｂを初期化する。つまり、学習行列作成部１０４ａは、学習行列Ａ、Ｂとして、空の行列を作成する。
【００６４】
ステップＳ４０２０において、学習行列作成部１０４ａは、後述する学習行列Ａ、Ｂの作成処理において使用する、バッチＩＤ１２０２（図３）のカウンタｉを初期化する。つまり、学習行列作成部１０４ａは、ｉ＝１とする。
【００６５】
ステップＳ４０３０において、学習行列作成部１０４ａは、後述する学習行列Ａ、Ｂの作成処理において使用する、時系列データＩＤ１２０３（図３）のカウンタｊを初期化する。つまり、学習行列作成部１０４ａは、ｊ＝１とする。
【００６６】
ステップＳ４０４０において、学習行列作成部１０４ａは、前処理運転データ１２のバッチＩＤ１２０２（図３）が「ｉ」のバッチに属する、［ｊ，ｊ＋Ｎ−１］の範囲の計測データ１２０４（図３）を用いて、学習行列Ａに追加する行ベクトル候補Ｘｉ，ｊを作成する。ここで、行ベクトル候補Ｘｉ，ｊは、以下の式（１）で表される。式（１）の右辺の要素であるベクトルｘｉ，ｊは、以下の式（２）で表される。
【００６７】
【数１】
0006573846.tif 000002
【００６８】
【数２】
0006573846.tif 000003
【００６９】
ここで、式（２）のベクトルの要素ｘｉ，ｊ，ｍは、バッチＩＤ１２０２が「ｉ」のバッチに属する時系列データＩＤ１２０３が「ｊ」の時刻における、ｍ番目の入力変数（つまり入力変数に含まれる計測信号のｍ番目のデータ）を意味する。すなわち、ｍは、診断モデルの入力変数を表す添字である。上述のように、診断モデルとして自己回帰モデルの一種のＶＡＲモデルを用いる場合には、入力変数の組合せ情報Ｍの大きさ｜Ｍ｜は２以上の整数となる。一方、ＶＡＲモデル以外の自己回帰モデルを用いる場合には．入力変数の組合せ情報Ｍの大きさ｜Ｍ｜は１となり、ｍも１となる。
【００７０】
ステップＳ４０５０において、学習行列作成部１０４ａは、カウンタｊがＪｉ−Ｎ＋１以下であるか否かを判定する。ここで、バッチＩＤ１２０２が「ｉ」のバッチに属する時系列データの個数（バッチｉのレコード数）をＪｉとする。カウンタｊがＪｉ−Ｎ＋１以下であれば（ステップＳ４０５０でＹＥＳ）、学習行列作成部１０４ａは、処理をステップＳ４０６０に進める。一方、カウンタｊがＪｉ−Ｎ＋１以下でなければ（ステップＳ４０５０でＮＯ）、学習行列作成部１０４ａは、処理をステップＳ４０９０に進める。
【００７１】
ステップＳ４０６０において、学習行列作成部１０４ａは、ステップＳ４０４０で作成した行ベクトル候補Ｘｉ，ｊを学習行列Ａに追加して学習行列Ａを更新する。学習行列Ａの更新処理は、図７に示されるように、更新前の学習行列７００の末尾に行ベクトル候補Ｘｉ，ｊを追加することにより行われる。
【００７２】
ステップＳ４０７０において、学習行列作成部１０４ａは、学習行列Ｂを更新する。学習行列Ｂの更新処理は、図８に示されるように、更新前の学習行列８００の末尾に、バッチＩＤ１２０２が「ｉ」のバッチに属するｊ＋Ｎ番目の出力変数ｙの値ｙｉ，ｊ＋Ｎを追加することにより行われる。ここで、出力変数ｙは、行ベクトル候補Ｘｉ，ｊの各要素と同様に前処理運転データ１２に含まれ、モデル構築情報１５により特定される、診断モデルの出力変数となる計測データ１２０４（図３）である。
【００７３】
ステップＳ４０８０において、学習行列作成部１０４ａは、カウンタｊを１加算して、処理をステップＳ４０４０に戻す。
【００７４】
ステップＳ４０９０において、学習行列作成部１０４ａは、カウンタｉがバッチ数Ｉ以下であるか否かを判定する。カウンタｉがバッチ数Ｉ以下であれば（ステップＳ４０９０でＹＥＳ）、学習行列作成部１０４ａは、処理をステップＳ４１００に進める。一方、カウンタｉがバッチ数Ｉ以下でなければ（ステップＳ４０９０でＮＯ）、学習行列作成部１０４ａは、処理をステップＳ４１１０に進める。
【００７５】
ステップＳ４１００において、学習行列作成部１０４ａは、カウンタｉを１加算して、処理をステップＳ４０３０に戻す。
【００７６】
上記ステップＳ４０３０〜Ｓ４１００のループを繰り返し実行することにより、最終的に得られた学習行列Ａは、以下の式（３）のように表記される。学習行列Ｂは、以下の式（４）のように表記される。式（４）から分かるように、学習行列Ｂは、厳密に言うと、ベクトルである。
【００７７】
【数３】
0006573846.tif 000004
【００７８】
【数４】
0006573846.tif 000005
【００７９】
ステップＳ４１１０において、モデル学習部１０４のパラメータ計算部１０４ｂは、得られた学習行列Ａ、Ｂを含む学習行列１３を用いて、モデルパラメータθを以下の公知の式（５）を用いて最小二乗法により計算して、図６の処理を終了する。
【００８０】
【数５】
0006573846.tif 000006
【００８１】
図９は、図６のステップＳ４０４０〜Ｓ４０８０のループによる学習行列Ａの更新を説明する図である。図１０は、本実施形態の学習行列Ａの行とバッチとの関係を説明する図である。図１１は、比較例の学習行列Ａ１の行とバッチとの関係を説明する図である。
【００８２】
図９の例では、バッチＩＤ１２０２が「ｉ」のバッチに属する時系列データ数Ｊｉは１００であり、モデル次数Ｎは１０である。図６のステップＳ４０４０において、［ｊ，ｊ＋Ｎ−１］の範囲の時系列データを用いて、行ベクトル候補が作成される。
【００８３】
図９において、カウンタｊがｊ＝１からｊ＝Ｊｉ−Ｎ＋１＝９１までの間、図６のステップＳ４０５０でＹＥＳと判断される。したがって、バッチＩＤ１２０２が「ｉ」のバッチに属する時系列データによって、９１行の学習行列が得られる。
【００８４】
ｊ＝Ｊｉ−Ｎ＋２＝９２になると、図９に示されるように、Ｎ（＝１０）個目の時系列データは、このバッチのデータではない。このとき、図６のステップＳ４０５０でＮＯと判断される。したがって、バッチＩＤ１２０２が「ｉ」のバッチに属する時系列データを用いる学習行列Ａの更新は終了する。
【００８５】
このように、図６のステップＳ４０４０〜Ｓ４０８０のループの処理によって、本実施形態の式（３）に示される学習行列Ａにおいて、図１０に示される同一バッチ内の連続する時系列データ１０００は、行単位で区切られる。
【００８６】
これによって、本実施形態によれば、データの連続性が保証され、モデルパラメータを適切に学習することができる。また、データの連続性が保証されるため、所望のモデル精度を得るために必要十分なサイズの学習データを、データの不連続性を考慮せずに用いることができる。これは、特に、短期間で起動と停止とを繰り返す空気圧縮機１のような回転機械の計測信号１０（前処理運転データ１２）を用いる場合でも、所望の精度を有する診断モデルを学習できることを意味する。
【００８７】
起動及び停止が頻繁に繰り返される空気圧縮機１以外の、起動及び停止が頻繁に繰り返されない機械では、大量の連続する時系列データを容易に得ることができる。このため、バッチの違いを考慮する必要が無い。
【００８８】
これに対して、本実施形態の空気圧縮機１のように、起動及び停止が頻繁に繰り返される回転機械では、バッチの相違を考慮しないと、図１１に示されるように、同一行内に異なるバッチに属するデータが含まれることとなる。この場合、図１１の位置１０１０でデータが不連続になるため、モデルパラメータを適切に学習できない。
【００８９】
図５に戻って、ステップＳ２０５０において、パラメータ計算部１０４ｂは、計算されたモデルパラメータ１４（つまり式（５）のモデルパラメータθ）を、診断モデルの学習結果として、モデル情報ＤＢ１０５に保存する。
【００９０】
ステップＳ２０６０において、特徴量計算部１０６は、前回の診断を実行した時から予め設定された診断周期を経過しているか否かを判定する。診断周期が経過していれば（ステップＳ２０６０でＹＥＳ）、処理はステップＳ２０７０に進められる。一方、診断周期が経過していなければ（ステップＳ２０６０でＮＯ）、処理はステップＳ２０００に戻る。なお、診断装置１００の初回起動後であって、診断モデルの学習結果が保存されている場合には、処理は無条件にステップＳ２０７０に進められてもよい。
【００９１】
ステップＳ２０７０（特徴量計算工程の一例）において、特徴量計算部１０６は、データ前処理部１０２から入力される前処理運転データ１６と、モデル情報ＤＢ１０５に保存されているモデル構築情報１５及びモデルパラメータ１４（つまり式（５）のモデルパラメータθ）とを用いて、現時点での計測信号に対するモデル特性評価を実行する。具体的には、特徴量計算部１０６は、時刻ｔにおける診断モデルの出力変数の予測値ｙ＊（ｔ）を、式（６）、（７）、（８）によって求める。
【００９２】
【数６】
0006573846.tif 000007
【００９３】
【数７】
0006573846.tif 000008
【００９４】
【数８】
0006573846.tif 000009
【００９５】
ここで、式（６）のベクトルＺ（ｔ）は、時刻ｔにおける自己回帰ベクトルを表す。式（７）のベクトルｘ（ｔ−ｎ）は、時刻ｔ−ｎにおける入力変数のベクトルを表す。式（８）のベクトルの要素ｘｍ（ｔ−ｎ）は、時刻ｔ−ｎにおける入力変数ｍの値を表す。また、ｎは、時刻ｎから過去の時刻を参照するための引数であり、正の整数で与えられる。
【００９６】
特徴量計算部１０６は、その特性評価結果を基に、診断条件ＤＢ１０７に保存されている実行条件パラメータ１８を用いて、診断用の指標である特徴量１９を計算する。具体的には、特徴量計算部１０６は、上記式（６）、（７）、（８）によって求めた出力変数の予測値ｙ＊（ｔ）と、時刻ｔにおける出力変数の実測値との誤差（モデル予測誤差）と、実行条件パラメータ１８とを基に、特徴量１９を計算する。特徴量１９を計算するときの実行条件パラメータ１８は、例えば該モデルの学習時に使用した正常時のデータから求めた、出力変数の平均値および標準偏差などが含まれる。
【００９７】
特徴量計算部１０６は、出力変数の予測値ｙ＊（ｔ）と、時刻ｔにおける出力変数の実測値との誤差（モデル予測誤差）を特徴量１９としてもよい。特徴量計算部１０６は、上記モデル予測誤差を用いて、別の特徴量１９を計算してもよい。
【００９８】
ステップＳ２０８０（状態診断工程の一例）において、状態診断部１０８は、特徴量計算部１０６により計算された特徴量１９と、診断条件ＤＢ１０７に保存された実行条件パラメータ１８とを用いて、空気圧縮機１の運転状態を診断する。状態診断部１０８は、診断結果２０をデータ出力Ｉ／Ｆ１１０に出力する。
【００９９】
特徴量１９が例えば上記モデル予測誤差の場合には、実行条件パラメータ１８は、閾値である。状態診断部１０８は、モデル予測誤差が閾値より大きい場合には、空気圧縮機１は異常であると診断する。状態診断部１０８は、モデル予測誤差が閾値以下の場合には、空気圧縮機１は正常であると診断する。
【０１００】
ステップＳ２０９０において、データ出力Ｉ／Ｆ１１０は、診断結果２０を含む表示データ２１を表示部１２０に適合する形式で生成する。データ出力Ｉ／Ｆ１１０は、生成した表示データ２１を表示部１２０に出力する。表示部１２０は、入力された表示データ２１を表示する。
【０１０１】
ステップＳ２１００において、特徴量計算部１０６は、特徴量１９を出力ログＤＢ１０９に保存する。状態診断部１０８は、診断結果２０を出力ログＤＢ１０９に保存する。
【０１０２】
ステップＳ２１１０において、診断装置１００は、動作を終了するか否かを判定する。終了する場合は（ステップＳ２１１０でＹＥＳ）、図５の動作を終了する。一方、終了しない場合は（ステップＳ２１１０でＮＯ）、処理はステップＳ２０００に戻る。
【０１０３】
ここで、本実施形態の診断装置１００の診断モデルとしてＶＡＲモデルが適用される場合の具体例が説明される。この場合の診断モデルおよび学習行列は、以下のように説明される。
【０１０４】
空気圧縮機１の計測信号１０の計測値が前処理運転データ１２に含まれる。このため、式（２）に示したベクトルｘｉ，ｊを構成する信号の種類は、上述のように、空気圧縮機１への吸込み空気圧、空気圧縮機１からの吐出空気圧、空気圧縮機１に潤滑用油を供給するための給油フィルタの前後の差圧、空気圧縮機１から流出する使用済み潤滑用油の不純物を分離する分離フィルタの前後の差圧、給油温度、油回収器の出口温度、空気圧縮機１の主電動機の軸受温度、空気圧縮機１の主電動機の巻線温度、空気圧縮機１のケーシング振動の振幅、温度によって変動する空気圧縮機１のロータ軸の軸方向長さを含む１０種類である。
【０１０５】
すなわち、式（１）に示したベクトル候補Ｘｉ，ｊは、上記１０種類の各信号からなるベクトルの、バッチｉに属する時刻ｊからｊ＋Ｎ−１までのＮ個の連続時系列データから構成され、最終的に学習行列Ａを構成する。
【０１０６】
同様に、学習行列Ｂにおいては、空気圧縮機１の計測信号１０のうちから１項目が出力変数ｙとして選択され、学習行列Ｂを構成する。
【０１０７】
また、このような学習行列Ａ、Ｂを用いて式（５）によって求められたモデルパラメータθは、本実施形態においては以下の意味を有するものと解釈できる。すなわち、空気圧縮機１のある計測信号（診断モデルの出力変数）の将来の計測値を予測する際に用いられる、該計測信号を含む他の計測信号から任意に選んだ計測信号群（診断モデルの入力変数）の、過去のＮ時刻前までの計測値と、該出力変数との因果関係の強さを表す。例えば、診断モデルの出力変数として空気圧縮機１の吐出空気圧を選択し、診断モデルの入力変数として、空気圧縮機１の吐出空気圧および吸込空気圧を選択し、モデル次数Ｎを２とした場合に、式（６）は、以下の式（９）のように書き下される。
【０１０８】
【数９】
0006573846.tif 000010
【０１０９】
ここで、式（９）において、ｐ＊ｏｕｔ（ｔ）は、吐出空気圧の時刻ｔにおける予測値である。ｐｏｕｔ（ｔ−１）、ｐｏｕｔ（ｔ−２）は、それぞれ、時刻ｔから１時刻前及び２時刻前の吐出空気圧の計測値である。ｐｉｎ（ｔ−１）、ｐｉｎ（ｔ−２）は、それぞれ、時刻ｔから１時刻前及び２時刻前の吸込み空気圧の計測値である。
【０１１０】
モデルパラメータθ１ｏｕｔ、θ１ｉｎ、θ２ｏｕｔ、θ２ｉｎは、それぞれ、４種類の信号ｐｏｕｔ（ｔ−１）、ｐｏｕｔ（ｔ−２）、ｐｉｎ（ｔ−１）、ｐｉｎ（ｔ−２）の過去の計測値に対する重み付け係数となっている。すなわち、モデルパラメータは、診断モデルの出力変数を予測するための、各入力変数の出力変数との因果関係の強さと言い換えることができる。
【０１１１】
以上説明されたように、本実施形態によれば、診断装置１００の学習方法によって、空気圧縮機１の計測信号から選んだ任意の計測信号を、該計測信号を含む計測信号群の過去の計測値を用いて、その因果関係を考慮した式（９）で示される線形演算によって予測する診断モデルが提供される。
【０１１２】
（その他）
（１）上記実施形態では、診断装置１００は、空気圧縮機１の運転状態を診断しているが、空気圧縮機に限られず、一般的な回転機械を診断対象としてもよい。
【０１１３】
診断装置１００は、例えば回転機械の一例であるガスタービンの運転状態を診断してもよい。この場合、計測信号として、燃料の流量、動翼の振動の振幅、燃焼器に送り込まれる燃焼用空気の圧力、燃焼器における燃焼ガスの温度、燃焼ガスの流量、タービン軸の軸受温度、発電電力などを含んでもよい。
【０１１４】
また、診断装置１００は、例えば回転機械の一例である蒸気タービンの運転状態を診断してもよい。この場合、計測信号として、蒸気流量、蒸気温度、蒸気圧力、ケーシング振動の振幅、タービン軸の軸受温度、発電電力などを含んでもよい。
【０１１５】
（２）上記実施形態では、空気圧縮機１から入力される計測信号１０は、空気圧縮機１への吸込み空気圧、空気圧縮機１からの吐出空気圧、空気圧縮機１に潤滑用油を供給するための給油フィルタの前後の差圧、空気圧縮機１から流出する使用済み潤滑用油の不純物を分離する分離フィルタの前後の差圧、給油温度、油回収器の出口温度、空気圧縮機１の主電動機の軸受温度、空気圧縮機１の主電動機の巻線温度、空気圧縮機１のケーシング振動の振幅、温度によって変動する空気圧縮機１のロータ軸の軸方向長さを含む１０種類としている。しかし、計測信号１０は、これら１０種類を全て含まなくてもよい。計測信号１０は、これらのうち少なくとも１つを含んでいればよい。
【０１１６】
（３）上記実施形態において、空気圧縮機１の空気圧縮方式は、特定の方式に限られない。例えば、ピストン往復運動によりシリンダ容積を変化させ圧縮するレシプロ方式でもよく、ケーシング内のスクリュローター回転によって圧縮するスクリュ方式でもよく、インペラの遠心力によって圧縮するターボ方式でもよい。
【０１１７】
（４）上記実施形態では、診断モデルとして、自己回帰（ＡＲ）モデル又はベクトル自己回帰（ＶＡＲ）モデルの例を示したが、これに限られない。診断モデルとして、自己回帰移動平均（ＡＲＭＡ）モデル、外部入力付自己回帰（ＡＲＸ）モデル、外部入力付自己回帰移動平均（ＡＲＭＡＸ）モデル、自己回帰和分移動平均（ＡＲＩＭＡ）モデルのいずれかを用いてもよい。
【０１１８】
（５）上記実施形態において、ユーザは、入力部１３０を用いて、計測信号１０として入力される１０種類の信号のなかから、所望の信号を出力変数として選択してモデル構築条件１５を生成し、生成したモデル構築条件１５をモデル情報ＤＢ１０５に登録可能に構成してもよい。また、ユーザは、入力部１３０を用いて、計測信号１０として入力される１０種類の信号のなかから、所望の信号を入力変数として選択してモデル構築条件１５を生成し、生成したモデル構築条件１５をモデル情報ＤＢ１０５に登録可能に構成してもよい。この実施形態によれば、空気圧縮機１の機能、仕様、経時的な運転特性などの変化に対応して、モデル構築条件１５を適宜調整することにより、診断モデルの精度を最適に維持することができる。
【０１１９】
（６）上記実施形態において、ユーザは、入力部１３０を用いて、任意の実行条件パラメータ１８を診断条件ＤＢ１０７に登録可能に構成してもよい。この実施形態によれば、空気圧縮機１の機能、仕様、経時的な運転特性などの変化に対応して、実行条件パラメータ１８を適宜調整することにより、診断性能を最適に維持することができる。
【０１２０】
（７）上記実施形態において、表示部１２０が音出力機能を備え、データ出力Ｉ／Ｆ１１０は、異常であることを示す診断結果２０が状態診断部１０８から入力された場合には、その旨を表す警告音を表示部１２０から出力させてもよい。この実施形態によれば、異常であるとの診断結果をユーザが視覚的に加えて聴覚的に確認することができ、空気圧縮機１の異常に対して適切に対処することが可能となる。
【０１２１】
（８）上記実施形態では、診断装置１００は、運転ＤＢ１０３、モデル情報ＤＢ１０５、診断条件ＤＢ１０７、出力ログＤＢ１０９を含んでいるが、これに限られない。例えば診断装置１００の外部に記憶装置を設け、この記憶装置が運転ＤＢ１０３、モデル情報ＤＢ１０５、診断条件ＤＢ１０７、出力ログＤＢ１０９を含むように構成してもよい。
【０１２２】
（９）本発明においては、診断装置１００内の処理は、上述の専用のハードウェアにより実現されるもの以外に、その機能を実現するためのプログラムを診断装置１００にて読取可能な記録媒体に記録し、この記録媒体に記録されたプログラムを診断装置１００に読み込ませ、実行するものであっても良い。
【０１２３】
診断装置１００にて読取可能な記録媒体とは、フロッピーディスク、光磁気ディスク、ＣＤ−ＲＯＭなどの移設可能な記録媒体の他、診断装置１００に内蔵されたハードディスクドライブ（ＨＤＤ）などの記憶部を指す。さらに、診断装置１００にて読取可能な記録媒体は、通信網を介してプログラムを送信する場合のように、短時間の間、動的にプログラムを保持するもの（伝送媒体もしくは伝送波）、その場合のサーバとなる診断装置１００内部の揮発性メモリのように、一定時間プログラムを保持しているものを含む。
【符号の説明】
【０１２４】
１  空気圧縮機
１０２  データ前処理部
１０４  モデル学習部
１０４ａ  学習行列作成部
１０６  特徴量計算部
１０８  状態診断部
【図１】
0006573846.tif 000011
【図２】
0006573846.tif 000012
【図３】
0006573846.tif 000013
【図４】
0006573846.tif 000014
【図５】
0006573846.tif 000015
【図６】
0006573846.tif 000016
【図７】
0006573846.tif 000017
【図８】
0006573846.tif 000018
【図９】
0006573846.tif 000019
【図１０】
0006573846.tif 000020
【図１１】
0006573846.tif 000021