(57)【特許請求の範囲】
【請求項１】
船体の参照方位及び参照位置を出力する軌道計画部と、該軌道計画部からの参照方位及び参照位置と、センサで検出された船体の方位及び位置から方位誤差及び位置誤差を演算する軌道航路誤差演算部と、該方位誤差及び位置誤差から船体の方位と位置とを参照方位及び参照位置に追跡させるべく指令舵角を出力するフィードバック制御部とを備えた船舶用自動操舵装置において、
前記フィードバック制御部は、センサで検出された方位ψ及び位置（ｘ、ｙ）並びに指令舵角δFBから推定方位誤差ψe^、推定位置誤差（ｘe^，ｙe^）、及び推定潮流ベクトル（ｄx^，ｄy^）の推定を行う推定器と、推定器から出力される推定誤差に対してフィードバックゲインを作用して指令舵角δFBを出力するフィードバックゲイン器とを備えており、該推定器は、方位誤差の推定を行う方位制御系推定手段と、航路誤差の推定を行う航路誤差系推定手段とからなり、
航路誤差系推定手段は、制御対象モデルの状態量を推定するための特性多項式として、
【数１】
0005042905.tif 000108
とし、前記フィードバックゲイン制御器を含む閉ループ系の特性多項式
【数２】
0005042905.tif 000109
の固有周波数ωyを航路制御系固有周波数としたときに、前記ωtを航路制御系固有周波数ωyのρt（＞１）倍に設定し、ωtoを航路制御系固有周波数ωyのρto（＜１）倍に設定することを特徴とする船舶用自動操舵装置。
【請求項２】
前記減衰係数ζt、ζtoは、１に設定することを特徴とする請求項１記載の船舶用自動操舵装置。
【請求項３】
航路誤差系推定手段は、前記λTCSの極配置を実現する推定ゲインを用いることを特徴とする請求項１または２記載の船舶用自動操舵装置。
【請求項４】
航路誤差系推定手段は、位置誤差ｘe、ｙe及び潮流ベクトルｄx、ｄyからなる航路誤差の状態量ベクトルＸt＝［ｘe、ｙe、ｄx、ｄy］Tに対する推定ゲイン行列Ｋtを
【数３】
0005042905.tif 000110
としたときに、ｋ22=ｋ11、ｋ32=-ｋ41、ｋ42=ｋ31，ｋ12＝ｋ21＝０、ｋ31 sinψR＝ｋ41 cosψRと設定し（ψRは参照方位である）、
【数４】
0005042905.tif 000111
と設定することを特徴とする請求項３記載の船舶用自動操舵装置。
【請求項５】
船首方向の推定位置誤差ｘeは、推定速度誤差ｘ・^に、Ｕ（cosψe^−１）（Ｕは船速）を修正したものを積分したものであることを特徴とする請求項１ないし４のいずれか１項に記載の船舶用自動操舵装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、航路制御系の船舶用自動操舵装置に関し、方位制御系を基本にした航路制御系のフィードバック制御系を持つ船舶用自動操舵装置に関する。
【背景技術】
【０００２】
船舶用自動操舵装置は、舵角を制御して設定方位に船首方位を追従させる方位制御系（HCS： Heading Control System）と、計画航路に船体位置を追跡させる航路制御系（TCS：Track Control System）とに分けられる。マイクロチップの高機能化、衛星測位システム（GNSS：Global Navigation Satellite System）の小型化・低コスト化・高精度化により位置情報が簡単に得られるようになったことに伴い、航路制御系の要求が高まっている。方位制御系の場合は、適時針路修正を必要とするが、航路制御系の場合は操船が不要で使い勝手がよい。
【０００３】
一般的に航路制御系の船舶用自動操舵装置は、図１に示すように軌道計画部１２、軌道航路誤差演算部１４、フィードバック制御部１６及び加算器１７を備える。計画航路に基づき軌道計画部１２が出力する参照方位ψRと参照位置ｘR、ｙRと、センサから検出される方位ψと位置ｘ、ｙとの誤差を軌道航路誤差演算部１４で求め、その誤差から保針時にフィードバック制御部１６が、船体の方位と位置とを追跡させるべくフィードバック舵角δFBを出力する。フィードバック制御部１６は、図２に示すように、推定器１８と、フィードバックゲイン器２０とからなる。
【０００４】
航路制御系では、航路誤差を直接制御できるアクチュエータを利用しないで船首方位を媒介にして間接的に制御するので劣駆動系である。
【０００５】
このフィードバック制御部の設計にあたっては、次数が大きいために、非特許文献１，２では、ＬＱＧベースで代数リカッチ式に帰着させている。このような方法による設計は全体を見通して行なわれるので解の信頼性が高い半面、閉ループ特性は定性的な把握になり、制御対象の特性を把握しないで、直接的に解を求めているので、制御対象に最適な解であるかは不明である。
【０００６】
例えば、非特許文献２では、船体は操舵によって船体横方向に対水の反力による横滑り速度を発生させるが、従来提案されている設計では、この横滑り速度を積極的に考慮して設計に反映させていない。そのため、横滑り速度が閉ループ系を不安定にして発散する原因になるおそれがある、という問題がある。
【０００７】
一方、本願発明者は、特許文献１、非特許文献３で方位制御系の設計方法を提案している。航路制御系は、基本的に、方位制御系に船体と潮流による横方向位置制御を加えたものであるから、航路制御系の性能は、かなりの部分で方位制御系の性能（方位安定性，外乱除去性など）に依存すると考えられる。そこで、本願発明者は、特許文献１で開示する方位制御系のフィードバック制御系を基礎にして、航路制御系のフィードバック制御系の設計を行うことにより、閉ループの定量的把握ができる船舶用自動操舵装置を非特許文献４及び特願２００７-５４５９９号で提案しており、また、旋回時に計画旋回に軌跡を乗せることができる船舶用自動操舵装置を特願２００７−５４６６７号で提案している。
【０００８】
しかしながら、これらの提案では、潮流の横方向（sway）成分のみの推定を行って、この横方向成分のみで制御を行っているために、潮流の船首方向（surge）成分を考慮に入れていないことによる誤差が生じる。特に、旋回時には船体に印加する潮流成分は船首方位によって変化するので、横方向（sway）の潮流成分の推定に基づく制御だけでは航路誤差が過渡的に生じる、という問題がある。
【０００９】
他方、特許文献２では、船位情報を用いて、潮流に関する状態量を推定することができる潮流推定方法及び装置を提案しており、航走体の位置に関する位置情報を受け付け、受け付けた位置情報に基づいて、航走体の運動を数学モデルで表現した潮流モデルと、航走体に作用する流体力を数学モデルで表現した流体力モデルとを含む運動モデルから状態量を推定することを提案するが、この提案では、具体的なオブザーバゲインの設定が定かでない。
【００１０】
【非特許文献１】Fossen, T. I.,「Marine Control Systems」, Marine Cybernetics AS, 392/394(2002)
【非特許文献２】Zwierzewicz, Z,「On the Ship Guidance Automatic System Design via LQG-Integral Control」, Manoeuvring and Control of Marine Craft 2003, IFAC.
【非特許文献３】羽根，「不確かさパラメータの下で閉ループ系代表根の減衰係数を設定する推定器設計」，第３４回制御理論シンポジウム，２１１/２１４（２００５）
【非特許文献４】羽根，米田，「方位制御系を基本にした航路制御系の設計」，第７回計測自動制御学会制御部門大会，２００７，８４−１−１
【特許文献１】特開２００６−４４４１１号公報
【特許文献２】特開２００５−１７２６１８号公報
【発明の開示】
【発明が解決しようとする課題】
【００１１】
そこで、本発明の目的は、潮流のベクトルを簡単な計算で推定することができる船舶用自動操舵装置を提供することである。
【課題を解決するための手段】
【００１２】
かかる目的を達成するために、本発明の請求項１に記載の発明は、船体の参照方位及び参照位置を出力する軌道計画部と、該軌道計画部からの参照方位及び参照位置と、センサで検出された船体の方位及び位置から方位誤差及び位置誤差を演算する軌道航路誤差演算部と、該方位誤差及び位置誤差から船体の方位と位置とを参照方位及び参照位置に追跡させるべく指令舵角を出力するフィードバック制御部とを備えた船舶用自動操舵装置において、
前記フィードバック制御部は、センサで検出された方位ψ及び位置（ｘ、ｙ）並びに指令舵角δFBから推定方位誤差ψe^、推定位置誤差（ｘe^，ｙe^）、及び推定潮流ベクトル（ｄx^，ｄy^）の推定を行う推定器と、推定器から出力される推定誤差に対してフィードバックゲインを作用して指令舵角δFBを出力するフィードバックゲイン器とを備えており、該推定器は、方位誤差の推定を行う方位制御系推定手段と、航路誤差の推定を行う航路誤差系推定手段とからなり、
航路誤差系推定手段は、制御対象モデルの状態量を推定するための特性多項式として、
【００１３】
【数１】
0005042905.tif 000002
とし、前記フィードバックゲイン制御器を含む閉ループ系の特性多項式
【００１４】
【数２】
0005042905.tif 000003
の固有周波数ωyを航路制御系固有周波数としたときに、前記ωtを航路制御系固有周波数ωyのρt（＞１）倍に設定し、ωtoを航路制御系固有周波数ωyのρto（＜１）倍に設定することを特徴とする。
【００１５】
請求項２に記載の発明は、請求項１記載の前記減衰係数ζt、ζtoが、１に設定されることを特徴とする。
【００１６】
請求項３に記載の発明は、請求項１または２記載の航路誤差系推定手段が、前記λTCSの極配置を実現する推定ゲインを用いることを特徴とする。
【００１７】
請求項４に記載の発明は、請求項３記載の航路誤差系推定手段が、位置誤差ｘe、ｙe及び潮流ベクトルｄx、ｄyからなる航路誤差の状態量ベクトルＸt＝［ｘe、ｙe、ｄx、ｄy］Tに対する推定ゲイン行列Ｋtを
【００１８】
【数３】
0005042905.tif 000004
としたときに、ｋ22=ｋ11、ｋ32=-ｋ41、ｋ42=ｋ31，ｋ12＝ｋ21＝０、ｋ31 sinψR＝ｋ41 cosψRと設定し（ψRは参照方位である）、
【００１９】
【数４】
0005042905.tif 000005
と設定することを特徴とする。
【００２０】
請求項５に記載の発明は、請求項１ないし４のいずれか１項に記載の船首方向の推定位置誤差ｘeが、推定速度誤差ｘ・^に、Ｕ（cosψe^−１）（Ｕは船速）を修正したものを積分したものであることを特徴とする。
【発明の効果】
【００２１】
本発明によれば、推定器を、舵角を制御して設定方位に船首方位を追従させる方位制御系を基礎とした方位制御系推定手段と、航路誤差系の航路誤差系推定手段とに分離して構成する。即ち、推定器の方位制御のみを行うときの制御対象モデルの状態量を推定するための特性多項式と、航路誤差系の推定位置誤差（ｘe^，ｙe^）、推定潮流ベクトル（ｄx^，ｄy^）を推定するための特性多項式とを分けることによって、それぞれの推定ゲインの決定が可能になる。これによって、既存の方位制御系についての推定器の構成をそのまま使用し、航路誤差系推定手段の構成を追加することで、推定器を構成することができるので、設計の手間を省き、コスト低減を図ることができる。
【００２２】
また、航路誤差系推定手段において、極配置による設計思想及び座標関係を利用した推定ゲインの関係から簡単な計算によって決まる推定ゲインを設定することで、潮流ベクトルの推定を簡単に行うことができる。
【００２３】
航路制御系の中に方位制御系がそのまま存在しているために、航路制御と方位制御との切り替えが円滑になり、切替時にその連続性を担保することができる。
【発明を実施するための最良の形態】
【００２４】
以下、図面を用いて本発明の実施の形態を説明する。
【００２５】
図１は、船舶用自動操舵装置と制御対象の全体のブロック図である。既述のように、船舶用自動操舵装置１０は、計画航路に船体位置を追跡させるために舵を制御する装置であり、軌道計画部１２、軌道航路誤差演算部１４、フィードバック制御部１６、加算器１７及び各パラメータを同定する図示しない同定器を備えている。誘導システム２２から計画航路及びセンサ類２６のスピードログから船速Ｕが軌道計画部１２に入力され、軌道計画部１２からは参照方位ψR、参照位置ｘR、ｙR及び旋回時にはフィードフォワード舵角δFFが出力される。
【００２６】
軌道航路誤差演算部１４には、ジャイロコンパスからの船首方位ψ、ＧＰＳ等の衛星測位システム（ＧＮＳＳ）からの船体の位置（ｘ，ｙ）といったセンサ類２６からの検出信号が入力され、軌道航路誤差演算部１４は、前記参照方位ψR、参照位置ｘR、ｙRとの比較を行い方位誤差と位置誤差（以下、航路誤差とも称する）（方位誤差と位置誤差を合わせて軌道誤差とも称する）を出力する。
【００２７】
船舶用自動操舵装置１０の閉ループ系は図１に示すように船体モデルと外乱モデルとからなる制御対象２４と、フィードバック制御部１６とから構成される。
【００２８】
フィードバック制御部１６は、図２に示すように、推定器１８とフィードバックゲイン器２０とからなる。軌道航路誤差演算部１４からの方位誤差及び位置誤差は、推定器１８に入力される。推定器１８において、方位誤差及び位置誤差から、方位誤差を含めた方位誤差系の状態量、位置誤差、潮流ベクトルを推定する。推定器１８は、方位誤差の推定を行う方位制御系推定器１８Ａと、航路誤差の推定を行う航路誤差系推定器１８Ｂとからなり、こうして推定器１８によって、外乱が除去された推定航路誤差である推定方位誤差ψe^及び横方向の推定位置誤差ｙe^が算出される。
【００２９】
フィードバックゲイン器２０でこの推定誤差に対してフィードバックゲインを乗じて指令舵角であるところのフィードバック舵角δFBを操舵機に出力する。
【００３０】
操舵機は指令舵角に比例した舵角を動かすため、船体は舵角によって旋回角速度を生じ方位、位置が変化する。旋回角速度の発生と共に、斜航角（横方向速度）が発生する。
【００３１】
以下、上記構成の詳細について説明する。
【００３２】
１．運動方程式
１．１  座標系
航路制御系で用いる座標系は、図３に示すように、以下の座標系から構成する。
・対地系（ＸＯＹ）：地球固定の緯度経度座標系で、ＧＮＳＳからの位置出力（ｘ、ｙ）に相当する。
・船体系（ＸBＧＹB）：船体固定の運動座標系で、船体の重心を原点とし、船首方位をＸB 軸とし、船体運動を定める。
・参照系（ＸRＯRＹR）：誘導システム２２により生成され指定された計画航路から定まる移動座標系である。
尚、座標系の回転極性は右ネジ方向を正とし、Ｚ軸方向は重力方向を正とする。座標系はＸ軸、Ｙ軸の２次元を用いる。
【００３３】
１．２  船体の運動方程式
前記制御対象である船体モデルを定めるために、船体の運動方程式を導出する。船体の運動方程式は、前進方向を除き、横方向と方位軸回りとの運動を扱うので、
【００３４】
【数５】
0005042905.tif 000006
【００３５】
を用いる。ここでＭx ，Ｍy はそれぞれｘ，ｙ方向の付加質量を含んだ質量を、Ｉz はｚ軸まわりの付加慣性モーメントを含んだ慣性モーメントを，Ｙ,Ｎはそれぞれｙ方向の流体力、ｚ軸まわりの流体モーメントを示し、添字は対応する変数を意味する。変数Ｕ,ｖ,ｒ,δはそれぞれ前進速度、横流れ（横滑り）速度、旋回（回頭）角速度と舵角とを示す。上式を整理すると
【００３６】
【数６】
0005042905.tif 000007
になり、ラプラス変換すると（ｓ はラプラス演算子を示す）
【００３７】
【数７】
0005042905.tif 000008
のようになる。ここで、横滑り速度と旋回角速度とが舵角を入力とした関係で結ばれる。(5)式を整理して、舵加速度δ・・(t)≒０とし、
【００３８】
【数８】
0005042905.tif 000009
を得る。ここで、
【００３９】
【数９】
0005042905.tif 000010
とおいている。
【００４０】
(6)式より、横方向（sway）の横滑り速度とヨー軸回りの旋回角速度との運動方程式は同一の形でかつ舵角による入力係数が異なるだけである。船体特性は舵を切ることによって、方位軸回りに旋回角速度を発生させ、同時に船体の横方向に対水の反力による横滑り速度を発生させる。(6)式において、実用的見地からｓ2 項による影響は無視できるためにｓ2 項を省略する。すると、(6)式は
【００４１】
【数１０】
0005042905.tif 000011
になる。ここで、
【００４２】
【数１１】
0005042905.tif 000012
であり、Ｔs = Ｔ1＋Ｔ2である。ここでＫs  は旋回力ゲイン、Ｋv  は横滑りゲイン、Ｔs3、Ｔv3は時定数である。
【００４３】
また、船舶用自動操舵装置は、舵角を通して制御量を入力するので、制御変数は角度単位の方が都合がよい場合があるため、
【００４４】
【数１２】
0005042905.tif 000013
を用いて横滑り速度ｖを斜航角βに変換すると、
【００４５】
【数１３】
0005042905.tif 000014
として導出される。ここで、Ｋβ は横滑り角ゲインまたは斜航角ゲイン（旋回力ゲインＫs と異符号となる）で、Ｔβ3は時定数で
【００４６】
【数１４】
0005042905.tif 000015
を示す。
【００４７】
船体モデルの横滑り速度ｖおよび斜航角β は旋回角速度ｒと共通の一次遅れ要素をもつので、(13)式からｒを基準にすると
【００４８】
【数１５】
0005042905.tif 000016
【００４９】
になる。これよりｖ， β はｒ にそれぞれゲインＫv ／Ｋs，Ｋβ／Ｋs によってほぼ比例することがわかる。Ｔv3 ＞Ｔβ3，Ｔs3＞Ｔβ3 の関係があるので，船体によりＴv3，Ｔβ3 が省略可能である一方で、Ｋv およびＫβ が必須なパラメータになる。
【００５０】
(1)式、(2)式はテーラー展開で１次の項までを考慮しているが、大きな変針角や旋回角速度に対応するためには、非線形項までを考慮する必要がある。その場合には、(9)式の運動方程式は、
【００５１】
【数１６】
0005042905.tif 000017
とするとよい。ここでαsは非線形係数であり、ゲイン、時定数と同様前記同定器によって与えられる。
【００５２】
１．３  外乱成分
外乱成分は、舵角オフセット成分、波浪成分と潮流成分を考える。舵角オフセットは風などに起因する方位軸まわりに作用する船体モーメントを舵角換算にしたものである。波浪は白色ノイズを入力した狭帯域フィルタの出力を方位相当したものである。潮流成分は対地系の速度成分である。式にまとめると、
【００５３】
【数１７】
0005042905.tif 000018
【００５４】
になる。ここでδo はほぼ一定値をした舵角オフセット成分を、ψw は波浪成分を、ν は白色ノイズＮ(0，1) を、Ｇw は波浪モデルの伝達関数を、Ｋw ，ζw ，ωw はそれぞれゲイン、減衰係数と中心周波数とを、ｄx，ｄy は対地系のそれぞれ北向き，東向きの潮流速度を、Ｕd ，ψd は潮流速度の大きさと方位とを、それぞれ示す。
【００５５】
制御対象は船体運動と外乱成分とからなる。外乱成分はほぼ一定値とする潮流成分と舵角オフセット成分，および平均値ゼロの有色モデルの波浪成分からなる。制御対象の構成を図５に示す。ただし同図では非線形項は省く。
【００５６】
１．４  対地速度と対地位置
船体の位置は、対地系で定義され、船体運動と潮流速度とからの対地速度を積分することで求めることができる。船体速度は対水速度で、潮流速度は対地速度でそれぞれ表す。
【００５７】
船体運動は舵角δを取ると角速度ｒ＝ψ・を生じるが、同時に横流れ速度ｖ＝ｙ・も生じる。その様子を、図４を用いて説明する。同図のように定常に旋回しているときまたは参照系が接線方位に一致しているとき、船体系は接線方位に対して斜航角βだけ内側に傾いて釣り合う。このとき船体速度の対地速度成分は
【００５８】
【数１８】
0005042905.tif 000019
になる。ここでｕx，ｖy は対地系のそれぞれ北向き，東向きの船体速度、ｕ，ｖは船体系の速度を示す。船体速度と斜航角との関係は、
【００５９】
【数１９】
0005042905.tif 000020
を用いる。
これより、対地系の速度は、船体速度と潮流速度との和になるから、
【００６０】
【数２０】
0005042905.tif 000021
になる。ここでｘ、ｙは船体位置でそれぞれ北向き、東向きを示す。なおｄx，ｄy は風が船体の上部構造物を押すことによる速度成分も含まれているとする。
位置と方位とはそれぞれ速度ｘ・、ｙ・と角速度ｒとを積分して
【００６１】
【数２１】
0005042905.tif 000022
得る。ここで（０）は初期値を意味する。
【００６２】
１．５  軌道航路誤差
軌道航路誤差演算部１４で求められる位置誤差及び方位誤差（航路誤差）は、参照系に対する船体系の位置誤差及び方位誤差であり、図３より、
【００６３】
【数２２】
0005042905.tif 000023
【００６４】
になる。ここでψ，ｘ，ｙは船首方位と対地位置とを、ψR，ｘR，ｙR は参照方位と参照位置とを、ψｅ，ｘｅ，ｙｅは方位誤差と航路誤差とをそれぞれ示す。ｘｅ，ｙｅは船体位置Ｇから参照系のＹＲ，ＸＲ 軸にそれぞれ垂線を下ろした距離に相当する。
【００６５】
制御対象の誤差モデルを構成するために、航路誤差の速度誤差を導入する。そこで、参照速度ｕR，ｖR を次のようにおく。
【００６６】
【数２３】
0005042905.tif 000024
【００６７】
参照系における船体系の速度誤差は、前記参照速度からの横流れ速度との偏差を参照系に座標変換したものとなる。
よって、速度誤差は、ψｅ＝ψ−ψR（(25)式より）の関係及び(21)式を用いて、
【００６８】
【数２４】
0005042905.tif 000025
になる。
一方、潮流速度による速度誤差は対地系成分を参照系成分に変換すると
【００６９】
【数２５】
0005042905.tif 000026
になる。
【００７０】
ここでｕd，ｖd は参照系の潮流速度成分を、ｄx，ｄyは対地系の潮流速度成分をそれぞれ示す。よって参照系速度誤差ｘe・，ｙe・は
【００７１】
【数２６】
0005042905.tif 000027
で与えられる。
【００７２】
図６は、制御対象のモデルを偏差系で表したものである。
潮流成分による斜航角βd は
【００７３】
【数２７】
0005042905.tif 000028
になる。数値例を挙げると、β ＝ ０ とすれば、Ｕ＝ 20[knot]，ψR ＝ 140[deg]，Ｕd ＝ 5[knot]，ψd ＝ 0[deg] のとき
【００７４】
【数２８】
0005042905.tif 000029
より、βd ＝ −11.24[deg] およびψ ＝ψR−βd ＝ 128.8[deg] になる。
【００７５】
２．保針制御系
２．１  制御対象
航路保針系の設計で用いる制御対象は参照系に対する誤差を扱う偏差系で定めることとする。船体運動は舵角を入力とし、旋回角速度と斜航角とを出力とするので、そのまま偏差系として利用される。
旋回運動と横滑り運動とは(9)式及び図５を参照すると
【００７６】
【数２９】
0005042905.tif 000030
【００７７】
になる。ここでｓはラプラス演算子を、ψ は船首方位を、δ は舵角を、δo は方位軸まわりの外乱モーメントを舵角換算したオフセットを、Ｋs，Ｋv，Ｔs，Ｔs3，Ｔv3は操縦性指数でそれぞれ旋回力ゲイン，横すべりゲインと３つの時定数とを示す。ｕ，ｖはそれぞれ船首方向（surge），横方向（sway）の船体速度を示す。非線形項αsｒ・3は省く。
【００７８】
船体系の速度誤差成分ｕｅ，ｖｅは(27)式から得られるが、ψeを微小項として近似し、(22)式を用いると、
【００７９】
【数３０】
0005042905.tif 000031
になる。ここで、ｕ≒Ｕ、ψeｖ≒０としている。
これにより、参照系の制御対象は、図５、(29)式、(31)式、(32)式、(33)式から、
【００８０】
【数３１】
0005042905.tif 000032
になる。ここで、ｒx、ｖx は状態変数であり、（図５（ｂ）、（ｃ）参照）、旋回角速度と横滑り速度とは
【００８１】
【数３２】
0005042905.tif 000033
になる。ここで、
【００８２】
【数３３】
0005042905.tif 000034
【００８３】
を示す。(35)式と(38)式では、時定数Ｔsが共通であり、図５（ｂ）、（ｃ）に示したように、同様の構成をとっている。よって、(38)式を消去することにすると、ｖx とｒx とは、
【００８４】
【数３４】
0005042905.tif 000035
となる。ここで、δo・≒０と仮定しており、一時遅れ要素を省略しても影響は少ないと仮定している。
よって、横滑り速度ｖは、
【００８５】
【数３５】
0005042905.tif 000036
と表すことができる。(37)式に、(45)式を入力することで、ｖの次数を減らすことができる。
以上の横滑り速度ｖは、斜航角βで置き換えることもでき、その場合、
【００８６】
【数３６】
0005042905.tif 000037
になる。ここで 、Ｋβ ，Ｔβ3 は操縦性指数でそれぞれ斜航角ゲイン、時定数を示す。これより
【００８７】
【数３７】
0005042905.tif 000038
になる。ここで
【００８８】
【数３８】
0005042905.tif 000039
を示す。ｒx・、βx・は時定数Ｔs が共通なので、斜航角の次数を低減すると、（31）式，（46）式より
【００８９】
【数３９】
0005042905.tif 000040
になる。ここで
【００９０】
【数４０】
0005042905.tif 000041
を示す。このとき
【００９１】
【数４１】
0005042905.tif 000042
になる。
【００９２】
２．２  線形近似による誤差
(33)式において、ｕeは０として線形になるように近似している。航路制御系はsway方向の誤差はゼロにする働きを持つが、surge方向については何の作用もしない。しかしながら、潮流成分を推定するためには、この影響を考慮する必要があるので、ｕeの近似による潮流計算値への影響を調べる。
(29)式の定常状態とすると、
【００９３】
【数４２】
0005042905.tif 000043
が成り立つ。ここで
【００９４】
【数４３】
0005042905.tif 000044
である。よって潮流値は
【００９５】
【数４４】
0005042905.tif 000045
【００９６】
になる。しかしながら線形化の近似のために(36)式のｘe・において、ｕe＝Ｕ（cosψe−１）成分による誤差が生じる。誤差を低減するために、(36)式のｘe・にｕe＝Ｕ（cosψe−１）を打ち消す修正操作を後述の推定器１８内に入れる。surge 方向の位置制御は実施せず、また位置のリセットは変針開始時に実施するために、ｕeの修正操作は、推定値ｘe^の増加を防ぐことからも好都合である。潮流の誤差は
【００９７】
【数４５】
0005042905.tif 000046
になる。数値例を示す。Ｕ ＝ 20[knot]，ψe ＝ 14[deg] とすればｕe ＝ −0.6[knot] になり、ｄx ＝ 5[knot] ならば12% の誤差に相当する。
【００９８】
３  フィードバック制御
方位誤差及びsway 方向の航路誤差は、図５と、（46）式，（37）式とから，波浪モデルを省くと
【００９９】
【数４６】
0005042905.tif 000047
になる。ここで
【０１００】
【数４７】
0005042905.tif 000048
を示す。ここで
【０１０１】
【数４８】
0005042905.tif 000049
の関係を用いる。
【０１０２】
フィードバックゲイン器２０を図７に示すように構成する。制御量であるフィードバック舵角δFBは、状態量が既知であるとして
【０１０３】
【数４９】
0005042905.tif 000050
【０１０４】
になる。ここでＫP，ＫD は方位制御系のそれぞれ比例ゲインと微分ゲインとを、Ｋy は航路制御系の航路ゲインを、αは外乱による修正量を示す。角速度はｒでなくｒx を利用するために上式に（39）式を代入すると
【０１０５】
【数５０】
0005042905.tif 000051
になる。これより
【０１０６】
【数５１】
0005042905.tif 000052
を得る。ここでγ は時定数Ｔs3 による修正係数で
【０１０７】
【数５２】
0005042905.tif 000053
を示す。
【０１０８】
上式より特性多項式に相当する項を展開して、方位制御系にのみ関連する項をまとめると、特性多項式λは、
【０１０９】
【数５３】
0005042905.tif 000054
になる。上式の右辺の第１項の（ｓ2＋２ζωnｓ＋ω2n）は、方位制御系における閉ループの特性多項式であり、ωn、ζはその方位制御系固有周波数、減衰係数を示す。
【０１１０】
この減衰係数ζ及び固有周波数ωnは、比例ゲインＫP  及び微分ゲインＫDとの関係を定めることができる。よって、設計パラメータを例えば比例ゲインＫP、減衰係数ζとしたときに、微分ゲインＫDと固有周波数ωnとを設計パラメータＫP、ζと同定器によって同定される船体パラメータＫs,Ｔs,Ｔs3とから一義的に決定することができ（特開２００７−１１８８２８号参照）、
【０１１１】
【数５４】
0005042905.tif 000055
として得られる。
【０１１２】
３．１  閉ループによる検討
基本特性を調べるためＴs3 ＝Ｔv3 ＝０ とし，ＣFB ＝ＫDs ＋ＫP とする。閉ループは
【０１１３】
【数５５】
0005042905.tif 000056
になるから、特性多項式λは
【０１１４】
【数５６】
0005042905.tif 000057
になる。航路ゲインの安定限界は
【０１１５】
【数５７】
0005042905.tif 000058
のようにＫyＵ の１次方程式の解になる。また特性多項式を開ループ伝達関数で表すと
【０１１６】
【数５８】
0005042905.tif 000059
になる。ゼロ点は不安定根ｓ ＝ １／Ｃβ (＞０) に、極は原点と安定な方位制御系の２根とになる。
【０１１７】
ゼロ点は横滑り速度に起因したもので、航路ゲインを極端に大きくすると閉ループが不安定になり発散する（ヨーイング発生）原因になる。
【０１１８】
根軌跡の例を図８に示す。開ループ伝達関数と３次方程式との性質から特性根は
1. Ｋy ＝０ のとき，根は極（記号×）である方位制御系根と原点根とになる
2. Ｋy ＞０ が大きくなると，方位制御系根は安定側から不安定側に原点根は左方向にそれぞれ移動し、この場合、方位制御系根は実軸でBreak-in point になり左右方向に別れて原点根は左方向にそれぞれ移動した後、Breakaway point から安定側から不安定側に移る
3. さらＫy ＞０ が大きくなると、根はゼロ点と左右方向とに向かって移動する
なる特徴をもつ。
【０１１９】
３．２  航路ゲインの求解
以上の航路ゲインＫyの影響を踏まえて、航路ゲインＫyを決定する。
閉ループ系の３次の特性多項式は(61)、(62)式より
【０１２０】
【数５９】
0005042905.tif 000060
になる。Ts3 ＞０ は横滑り角によるＣβ 係数を低減させ安定化に作用する。一方、設計のための特性多項式は
【０１２１】
【数６０】
0005042905.tif 000061
【０１２２】
として定める。ここでａは負の実数を，ζy，ωy はそれぞれ航路制御系減衰係数と固有周波数とを示す。航路制御系の設計パラメータをζy として、ａ，ωy，Ｋy を求める。閉ループと設計とのそれぞれの特性多項式の係数を比較すると
【０１２３】
【数６１】
0005042905.tif 000062
になる。整理すればωy に関する４次方程式
【０１２４】
【数６２】
0005042905.tif 000063
【０１２５】
を得る。ωy の４次と１次との係数は負になり、０次の係数は正になるので，少なくとも上式の根のひとつは負になっている。(78)式から求めたωy から、(75)〜(77)式の関係を用いれば、
【０１２６】
【数６３】
0005042905.tif 000064
が順次求まる。
【０１２７】
（78）式のωyの４次方程式は、代数解法で解ける。Ｔs3，Ｔβ3 に関連する高次項の近似解はＣβ＞Ｔs3 を仮定すると
【０１２８】
【数６４】
0005042905.tif 000065
【０１２９】
になり、式（78）の解の中で最大根に相当し、すぐに減衰する特性を示す。そのため、４次式から近似解を除いた３次式の特性が支配的になる。ωyを前述したＴs3，Ｔβ3を省いた３次の開ループ伝達関数の特性を利用して決め、設計パラメータζy は１／√２≦ζy≦１からBreak-inまたはBreakaway pointsとなりＫy に対する感度が極めて高くパラメータ誤差の影響を受け易いζy ＝１を避けて、ζy ＝ 0.9を選ぶ。Ｔs3，Ｔβ3 を省いた(74)式の３次系の特性多項式の根軌跡は図８に既に示した通りである。同図において大きさの異なる３つの黒丸はζy ＝ 0.9 として(78)式でＴs3，Ｔβ3 を省いた
【０１３０】
【数６５】
0005042905.tif 000066
【０１３１】
を解いたときの３つのωyに対応している。いずれも実根はゼロ点と安定根とになり、共役根は軌跡と直線との交点にある。直線の傾きは安定根ならθ ＝ cos-1ζy，不安定根ならq ＝ cos-1(−ζy) になる。３つの根は大きさが小さい順に航路ゲインＫy1 ＝ ０ ＜ Ｋy2 ＜Ｋy3 になる。方位制御系のζ ＝ 0. 9 よりＫy1 ＝０ のとき根は極と一致している。３次方程式の解と根軌跡との性質から航路ゲインＫy を決定する方法は
1. まず（82）式の解の中から負の解を除いた解を求める
2. 正の解からａ＞ωy になるような根ωy を選ぶ。具体的には正の解で最も小さい値が相当する。ａ＞ωy を用いる理由は閉ループ系の応答速度を上げるために代表根の位置を原点からできるだけ離すことによる
になる。上の方法はそのまま４次の（78）式についても全く同様に成立し、得られた根ωyから（79）式，（80）式を用いてそれぞれ実根ａと航路ゲインＫy とを求めることができる。
【０１３２】
フィードバックゲイン器２０は、同定器で各係数Ｃβ、Ｔβ3、Ｔs3、ωnが与えられる毎に、（78）式の解を解き、上記手順にて決定された航路ゲインＫyにて指令舵角を出力する。
【０１３３】
３．３  船速の変化
船速が変化した場合の航路ゲインＫyの様子を調べる。船体運動のsway 方向速度は対水速度Ｕに比例するため、速度が変化したらＫy はＵ と連動させて適切に更新させる必要がある。
船速変化に対する船体パラメータは
【０１３４】
【数６６】
0005042905.tif 000067
がほぼ成立するとする。Ｋy は式（71）の安定限界を用いると
【０１３５】
【数６７】
0005042905.tif 000068
【０１３６】
になる。数値例を挙げると，ζ ＝ 0. 9，ωn ＝ 0. 05[rad/s]，Ｃβ ＝ 30[s]，Ｕ ＝ 10[m/s]，速度変化Ｕ＊ ＝ 10, 5[m/s] のときそれぞれＫy ＝ 2. 43，3. 83[rad/m] になる．速度が10 から５[m/s]に低下すると航路ゲインは1.6 倍程度大きくする必要がある。
【０１３７】
３．４  外乱項による影響と対策
外乱入力に対する航路誤差の定常値を求め、その誤差を修正する。平均値がゼロであるψw を省き、定常値に影響しないＴs3，Ｔv3 をそれぞれゼロとすれば、舵角オフセットδo と潮流成分オフセットｖd（(28)式参照）との外乱入力に対する方位誤差ψe と航路誤差ｙe との定常解を求める。フィードバック舵角
【０１３８】
【数６８】
0005042905.tif 000069
を（61）式の方位誤差に代入すると
【０１３９】
【数６９】
0005042905.tif 000070
になる。方位誤差と航路誤差とはそれぞれ
【０１４０】
【数７０】
0005042905.tif 000071
になる。これより
【０１４１】
【数７１】
0005042905.tif 000072
を得る。よって一定値の外乱入力δo(s) ＝δo(0)／s，vd(s) ＝ vd(0)／s，修正量
【０１４２】
【数７２】
0005042905.tif 000073
に対する定常値は
【０１４３】
【数７３】
0005042905.tif 000074
を用いて
【０１４４】
【数７４】
0005042905.tif 000075
【０１４５】
になる。ここで(0) は一定値を意味する。ψe は積分効果をもつ航路系のためδo(0) による誤差を生じないが、ｖd(0) による偏差または斜航角を生じる。ｙe はδo(0)，ｖd(0) により誤差を生じる、ｙe をゼロにする修正量の定常値はα(0) ＝−δo(0)−KPvd(0)／U になる。
よって修正量
【０１４６】
【数７５】
0005042905.tif 000076
【０１４７】
として導ける。修正量α をフィードバックし、船首方位ψ が参照方位ψR 基準から見ればvd(0)／U の逆方向に針路を取ることで、船体のsway 方向におけるvd(0)／U成分と偏差ψeとが相殺する関係になる。数値例を挙げる．ＫP ＝１，Ｋy ＝ 0.0007[rad/m]＝0.04[deg/m] とすると，vd(0)／U ＝ 10[deg] のときｙ＝ 250[m] になり，δo ＝ 3[deg] のときｙ＝ 75[m] になる。
【０１４８】
したがって、外乱成分δo(s)，vd(s) を推定器１８によって推定し、上式の修正量をフィードバックすることで、外乱成分に起因した航路誤差を相殺することができる。
【０１４９】
４  推定器
次に、推定器１８の構成について説明する。まず、特許文献１で開示する方位制御系の５次推定器について説明し、次に、その５次推定器に基づき、航路制御系の推定器を構成することを考える。
【０１５０】
４．１  方位制御系の５次推定器
方位制御系の５次推定器の構成例を図９に示す。この推定器は、波浪モデルと舵角オフセットとを組み込んだものであり、(34)式、(35)式及び検出方位誤差ψe-
【０１５１】
【数７６】
0005042905.tif 000077
を用いて、
【０１５２】
【数７７】
0005042905.tif 000078
と表すことができる。ここで、ｋｉ，ｉ＝１〜５は推定ゲインを、ξ^は推定状態量を、ψW^は推定波浪を、δO^は推定舵角オフセットを、それぞれ示す。上式を行列で表すと、
【０１５３】
【数７８】
0005042905.tif 000079
となる。
ここで、
【０１５４】
【数７９】
0005042905.tif 000080
であり、ε は微小項を示す。
この特性多項式を求めるために、
【０１５５】
【数８０】
0005042905.tif 000081
とする。航路制御系の特性多項式は、
【０１５６】
【数８１】
0005042905.tif 000082
【０１５７】
に定める。ここでλhhは船体モデルの状態推定に、λhwは波浪モデルの状態推定に、そしてλhoは舵角オフセットモデルの状態推定にそれぞれ対応し、ζh，ζhwとωh，ωhw，ωho とは方位制御系のそれぞれ減衰係数と固有周波数とで、ζhw，ωhwは、図示しない波浪同定器によって検出方位から同定されるもので、
【０１５８】
【数８２】
0005042905.tif 000083
を示す。ここでωn は方位制御系の操舵系固有周波数（(68)式参照）を、ρh，ρho は推定係数で推定速度やパラメータ不確かさに対する許容度から決定され
【０１５９】
【数８３】
0005042905.tif 000084
とする。
【０１６０】
特性多項式の根を極配置することにより求める推定ゲインＫh＝［ｋ１ｋ２ｋ３ｋ４ｋ５］Ｔは、
【０１６１】
【数８４】
0005042905.tif 000085
になる。ここで、添字(・)-1 は逆行列を意味し、添字(・)T は転置行列を意味する。
【０１６２】
【数８５】
0005042905.tif 000086
になる。
【０１６３】
４．２  航路制御系の制御対象
推定器１８は、極配置を用いるために、確定系オブザーバによって構成される。制御対象は方位制御系と航路誤差系とをまとめて扱うより，おのおの分離するほうが構成する上で便利である。方位制御系は航路誤差系から影響を受けないので、方位制御系推定器１８Ａと、航路誤差系推定器１８Ｂとから構成し、図６の構成は、方位制御系推定器１８Ａとしてそのまま利用し、新たに航路誤差系推定器１８Ｂを追加する。制御対象を方位制御系と航路誤差系とに分割したブロック行列で表現すると、
【０１６４】
【数８６】
0005042905.tif 000087
になる。添字 h,t は、それぞれ方位系，航路系に関することを意味し、航路系は
【０１６５】
【数８７】
0005042905.tif 000088
を示す。ここでε は微小項を、Ａthにおいて
【０１６６】
【数８８】
0005042905.tif 000089
を用いることも可能である。
【０１６７】
４．３  推定ゲイン
航路誤差系を含めた推定器１８の固有値または特性多項式は９次になる。この９つの根を極配置するための推定ゲインを定める。未知数の９×２行列の推定ゲインは特性多項式の９つの係数方程式によって求めることになるが、方程式の数より未知数が多いので不定状態になる。そこで、前述の通り、推定器１８を、図９に示した方位制御系推定器１８Ａと、図１０に示した航路誤差系推定器１８Ｂとに分離して、極配置する推定ゲインを設計する。方位制御系の５次相当部分は、４．１章に示した通り、特許文献１において設計可能であるから、航路誤差系推定器１８Ｂの４次の極配置を設計する。
【０１６８】
推定器の構成は(93)式と、(106)式及び(107)式から
【０１６９】
【数８９】
0005042905.tif 000090
【０１７０】
になる。ここでＫt は航路制御系（ＴＣＳ）の推定ゲインを示す。航路誤差系の推定器を設計し、Ｋt を求める。この推定器の特徴は、
・対地系の北向き（x）と東向き（y）との潮流成分を推定する
・参照系のｘe・とｙe・とに方位角を介して潮流成分が印加する
・変化する方位角のために、非線形系である
である。推定ゲインは対角化されて
【０１７１】
【数９０】
0005042905.tif 000091
のように定める。ここでＫh は方位制御系の推定ゲインを、Ｋt は航路誤差系の状態量
【０１７２】
【数９１】
0005042905.tif 000092
を推定するための推定ゲインを示す。
推定器の特性多項式を求めるために、
【０１７３】
【数９２】
0005042905.tif 000093
を定める。ここでｓはラプラス演算子を、Ｉは単位行列を意味しその括弧内の数字は次数を示す。航路制御系の特性多項式は上式より
【０１７４】
【数９３】
0005042905.tif 000094
になる。よって航路制御系の特性多項式は方位制御系の行列式 det(λh) と航路誤差系の行列式 det(λt ) との積になり、det(λh) は、４．１章で説明した。
一方 det(λt ) はこれから設計されるもので
【０１７５】
【数９４】
0005042905.tif 000095
に定める。ここでζt，ζto とωt，ωto とは航路誤差系のそれぞれ減衰係数と固有周波数とを
【０１７６】
【数９５】
0005042905.tif 000096
を示す。λTCS の固有周波数を
【０１７７】
【数９６】
0005042905.tif 000097
と定める。ここでωy は航路制御系固有周波数（(74)式参照）を、ρt，ρto は推定係数で、ρt＞１，ρto＜１とし、例えば、
【０１７８】
【数９７】
0005042905.tif 000098
と定める（(104)式参照）。
【０１７９】
さてλTCS の極配置を実現する推定ゲインを求める。方程式は上記のａ0，ａ1，ａ2，ａ3 の４つであり、推定ゲインは２×４＝８個である。そのため解析解の導出は難しい。λHCS の場合は５つの方程式から５つの推定ゲインを求める条件であった。本発明による設計では、推定器の構造に着目することによって、極めて簡単な形で推定ゲインを導出することを実現する。
【０１８０】
推定モデルのＡt は対地系から参照系までの座標変換の行列をもっている。そのために方位によって特性多項式の係数が変化する。推定ゲインが対称的な方位の要素をもつことで、係数が一定になることが予想される。これより推定ゲインの要素を
【０１８１】
【数９８】
0005042905.tif 000099
とおいて
【０１８２】
【数９９】
0005042905.tif 000100
を用いて、推定ゲインを用いた特性多項式を求めると
【０１８３】
【数１００】
0005042905.tif 000101
になる。ここで
【０１８４】
【数１０１】
0005042905.tif 000102
を示す。
【０１８５】
さて式(129)式において、係数の要素であるＣ3 ＝０ なる関係を導入することは可能である。それを利用すると
【０１８６】
【数１０２】
0005042905.tif 000103
【０１８７】
になり、残りの変数はＣ1，Ｃ2，ｋ12ｋ21 が残る。Ｃ1 はｓ3 の係数ａ3 に対応させて求める。Ｃ2，ｋ12ｋ21 は３つのｓ2，ｓ1，ｓ0 の係数ａ2，ａ1，ａ0 に対応するために、一意的に定まらない。これはＣ1，Ｃ2，Ｃ3，ｋ12ｋ21 の設定方法では求めることが難しいことを表している。ＨＣＳは任意のζh，ζhw を指定することができるが、ＴＣＳは特定のζt，ζto を指定することによって極配置を実現する。すなわち
【０１８８】
【数１０３】
0005042905.tif 000104
を定める。この理由は
【０１８９】
・ 上記の設定により、設計パラメータがωt、ωto の２つになる。
・ 推定速度の点から見ると上記の設定は適切ではないようだが、応答はωt、ωto を調節すれば改善される。
・ この推定器はサンプル毎に方位角に対応して推定ゲインを計算するので、ゲイン計算量を少なくする必要がある。
である。そのためにζt、ζto の任意指定は必ずしも必要がない。上記の設定によって
ｋ12 ＝ ｋ21 ＝ ０  が得られるので
【０１９０】
【数１０４】
0005042905.tif 000105
になる。また(124)式の係数とＣ3＝０とを用いると
【０１９１】
【数１０５】
0005042905.tif 000106
に得る。これより推定ゲインは
【０１９２】
【数１０６】
0005042905.tif 000107
になる。よって航路誤差系の推定ゲインは極めて単純な計算によって得られる。
【０１９３】
４．４  まとめ
以上のようにして構成される推定器１８の航路誤差系推定器１８Ｂの構成を図１１に示す。ここでの推定ゲインは(135)式で決まり、また、(135)式におけるωt、ωtoは(125)式で決まり、(125)式における航路制御系の固有周波数ωyは、(73)式に示した航路制御の閉ループ系の特性多項式の固有周波数であり、具体的には、(78)式の４次方程式の解のうちの負の解と、最大の正の解を除き、残りの２つの正の解のうちで、特性多項式の根の代表根が原点からより離れる小さい方の解の値となる。
【０１９４】
また、航路誤差系推定器１８Ｂにおいて、２．２章で説明した、船首方向の推定位置誤差ｘeの推定速度誤差ｘ・^に、Ｕ（cosψe^−１）（Ｕは船速）の修正を加えることにより、線形近似による誤差を修正することができる。
【０１９５】
５  まとめ
以上説明したように、推定器１８を、舵角を制御して設定方位に船首方位を追従させる方位制御系を基礎とした方位制御系推定器１８Ａと、航路誤差系の航路誤差系推定器１８Ｂとに分離して構成し、航路誤差系推定器１８Ｂにおいて、航路誤差系の推定位置誤差（ｘe^，ｙe^）、推定潮流ベクトル（ｄx^，ｄy^）を推定し、航路誤差系推定器１８Ｂでは、極配置による設計思想及び座標関係を利用した推定ゲインの関係から簡単な計算によって決まる推定ゲインを設定することで、船体の位置誤差と潮流ベクトルの推定を行うことができる。これらの推定値は、閉ループ安定性確保と、航路誤差の低減に資することができる。
【０１９６】
基本的に方位制御系を基礎に、航路制御系の設計を行っているために、設計工数を削減することができ、また、航路制御系と方位制御系との間での動作切換による不連続減少を防止することができる。
【図面の簡単な説明】
【０１９７】
【図１】本発明による船舶用自動操舵装置の全体構成を表すブロック図である。
【図２】フィードバック制御部の構成を表すブロック図である。
【図３】航路制御系で用いる座標系を表す説明図である。
【図４】舵を取ったときに発生する横流れ速度を表す説明図である。
【図５】制御対象のモデルを表すブロック図である。
【図６】制御対象のモデルを偏差系で表したブロック図である。
【図７】フィードバックゲイン器の構成を表すブロック図である。
【図８】Ｔs3＝Ｔβ3＝０とおいたときの特性多項式の根軌跡を表す図である。
【図９】方位制御系推定器のブロック図である。
【図１０】航路誤差系推定器のブロック図である。
【図１１】航路誤差系推定器のブロック図である。
【符号の説明】
【０１９８】
１０  船舶用自動操舵装置
１２  軌道計画部
１４  軌道航路誤差演算部
１６  フィードバック制御部
１８  推定器
１８Ａ  方位制御系推定器
１８Ｂ  航路誤差系推定器
２０  フィードバックゲイン器
【図１】
0005042905.tif 000112
【図２】
0005042905.tif 000113
【図３】
0005042905.tif 000114
【図４】
0005042905.tif 000115
【図５】
0005042905.tif 000116
【図６】
0005042905.tif 000117
【図７】
0005042905.tif 000118
【図８】
0005042905.tif 000119
【図９】
0005042905.tif 000120
【図１０】
0005042905.tif 000121
【図１１】
0005042905.tif 000122