(57)【特許請求の範囲】
【請求項１】
目標体から放射される音波の到来方位及び周波数を音響センサにより時系列的に観測し、蓄積した時系列観測方位データ及び時系列観測周波数データを用いて、任意の基準時刻における前記目標体の位置、針路及び速力を含む運動ベクトルを推定する目標運動解析方法において、
(a)目標体の変針変速の検出に応じて、前記変針変速を検出した場合には、前記変針変速時刻を仮定し、
(b)前記仮定した変針変速時刻の前後の区間のそれぞれを前記目標体は等速直進運動している２レグ運動と見なして、以前に推定した変針変速時刻の中で最も新しいデータである最新の変針変速時刻から前記仮定した変針変速時刻までの区間と、前記仮定した変針変速時刻から最新の観測時刻までの区間とのそれぞれの区間に等速直進モデルを適用して非線形の最小二乗法により運動パラメータを推定する２レグ解析を行って、前記仮定した変針変速時刻前後の時系列推定方位データ及び時系列推定周波数データを求め、
(c)前記仮定した変針変速時刻前後の前記時系列観測方位データと前記時系列推定方位データとの残差及び前記時系列観測周波数データと前記時系列推定周波数データとの残差の和を評価値として、前記任意の基準時刻における前記目標体の位置、前記変針変速前後の針路及び速力を、導関数を用いた最適化手法により推定する目標運動解析方法。
【請求項２】
前記(a)での変針変速時刻の仮定は、変針変速の検出時刻から過去方向に変化させながら、前記(b)及び(c)の手順を前記評価が改悪方法になるまで繰り返し、前記(a)で仮定した変針変速時刻のうち前記評価値が最小となる時刻を変針変速時刻として推定する請求項１記載の目標運動解析方法。
【請求項３】
目標体から放射される音波の到来方位及び周波数を時系列的に観測する音響センサ、前記音響センサにより観測された前記音波の到来方位及び周波数を、時系列観測方位データ及び時系列観測周波数データとして格納するデータベース、及び
(a)目標体の変針変速の検出に応じて、前記変針変速を検出した場合には、前記変針変速の時刻を仮定する処理、(b)前記仮定した変針変速時刻の前後の区間のそれぞれを前記目標体は等速直進運動している２レグ運動と見なして、以前に推定した変針変速時刻の中で最も新しいデータである最新の変針変速時刻から仮定した変針変速時刻までの区間と、前記仮定した変針変速時刻から最新の観測時刻までの区間とのそれぞれの区間に等速直進モデルを適用して非線形の最小二乗法により運動パラメータを推定する２レグ解析を行って、前記区間毎に、前記時系列観測方位データを基にした時系列推定方位データと前記時系列観測周波数データを基にした時系列推定周波数データを求める処理、および(c) 前記仮定した変針変速時刻前後の前記時系列観測方位データと前記時系列推定方位データとの残差及び前記時系列観測周波数データと前記時系列推定周波数データとの残差の和を評価値として、任意の基準時刻における前記目標体の位置、前記変針変速前後の針路及び速力を、導関数を用いた最適化手法により推定する処理を実行するプロセッサを有する目標運動解析装置。
【請求項４】
前記プロセッサは、前記仮定した変針変速の時刻を変針変速の検出時刻から過去方向に変化させながら、前記(b)及び(c)の処理を前記評価が改悪方法になるまで繰り返した結果、前記評価値が最小となる変針変速時刻を出力する請求項３記載の目標運動解析装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、目標体から放射される音波を時系列データとして観測し、目標の位置、針路、速力等を推定する目標運動解析方法及び装置に関する。
【背景技術】
【０００２】
目標運動解析方法は移動可能な観測体に設けられた音響センサにより、目標体から放射される音波の到来方位及び周波数を時系列観測データとして蓄積し、前記時系列観測データと時系列推定データの残差が最小となる目標体の運動モデルのパラメータを推定するものである。
【０００３】
これまでの目標運動解析方法は目標体の運動モデルを等速直進運動として扱っており、推定するパラメータは基準時刻における目標体の位置座標、針路、速力、固有周波数であった(特許文献１)。
【０００４】
前述した時系列観測データと時系列推定データの残差を最小化し、目標体の運動ベクトルを推定する問題は、非線形最適化問題として知られており、カルマンフィルタによる方法、最小二乗法による方法、遺伝的アルゴリズムによる方法等があった。
【０００５】
また、等速直進運動ではない目標体、すなわち変針変速を行うような目標体については、変針変速を検知し、オペレータに提示する方法があった (特許文献2)。
【０００６】
【特許文献１】特開平10−62511公報
【特許文献２】特開平10−62507公報
【発明の開示】
【発明が解決しようとする課題】
【０００７】
前述した技術では、目標体の運動モデルを等速直進運動として扱っているため、目標体が変針変速を行うと、変針変速前後の目標の運動を正しく解析できないという問題があった。
【０００８】
図１１は前述した従来技術による目標運動解析のシミュレーション結果である。1101は観測体の航跡、1102は目標体の航跡、1103は解析周期30秒として目標運動解析を実施した時の目標体の推定位置座標履歴、1104は観測体と目標体の相対距離、1105は目標運動解析結果による観測体と目標体の推定相対距離、1106は目標体の速力、1107は目標運動解析結果による目標体の推定速力、1108は目標体の針路、1109は目標運動解析結果による目標体の推定針路を表している。シミュレーション開始後900秒で目標体は針路を180°から150°に変針しており、目標体変針以降においては目標運動解析を正しく実施できていないことが分かる。
【０００９】
また、目標体の変針変速を何らかの検定手段により看破し、変針変速後の時系列データを用いて、解析を行う方法においては、変針変速直後は解析が不能となり、解析結果が不安定となるという問題があった。これは、変針変速後は解析に使用するデータ数が減少するためである。
【００１０】
図１２は何らかの検定手段により、変針変速を検出し、変針変速検出後は変針変速時刻以降の観測データを用いて目標の運動を解析する技術による目標運動解析のシミュレーション結果である。1201は観測体の航跡、1202は目標体の航跡、1203は解析周期30秒として目標運動解析を実施した時の目標体の推定位置座標履歴、1204は観測体と目標体の相対距離、1205は目標運動解析結果による観測体と目標体の推定相対距離、1206は目標体の速力、1207は目標運動解析結果による目標体の推定速力、1208は目標体の針路、1209は目標運動解析結果による目標体の推定針路を表している。シミュレーション開始後900秒で目標体は針路を180°から150°に変針しており、目標体の変針以降は一度解析不能状態となり、約600秒間は解析結果が不安定な状態となっていることが分かる。
【００１１】
前述の課題を解決する方法として、目標体の運動モデルとして等速直進運動ではなく、変針変速を考慮して、2区間以上の等速直進区間の組み合わせとして目標体の運動をモデル化する方法が考えられる（以降、このモデルを等速直進モデルに対して、変針変速モデルと呼ぶ）。しかし、変針変速モデルの時系列推定方位関数及び時系列周波数関数は不連続なものとなってしまい、カルマンフィルタあるいは最小二乗法のように評価値の導関数を使用する最適化手法を適用することができない。遺伝的アルゴリズムへの変針変速モデルの適用は比較的容易であるが、遺伝的アルゴリズムは推定パラメータが増加すると、解析を収束させるための計算量は飛躍的に増加することが知られている。
【００１２】
そのため、計算量の大幅な増加を伴わずに変針変速前後の目標体の運動ベクトルを安定して出力する目標運動解析方法及び装置が求められている。
【課題を解決するための手段】
【００１３】
本発明は、目標体から放射される音波の到来方位及び周波数を音響センサにより観測した、時系列観測方位データ及び時系列観測周波数データから、目標体の位置及び針路、速力を推定する目標運動解析方法である。その目標運動解析方法は、目標体の変針変速を検出した上で変針変速時刻を推定し、時系列観測データを用いて、基準時刻における目標体位置座標及び、等速直進区間毎の目標体の針路、速力を推定パラメータとする非線形最適化問題を解くことである。
【００１４】
より具体的な態様は、まず目標体の変針変速の時刻を仮定し、仮定した変針変速時刻の前後の区間のそれぞれを前記目標体は等速直進運動していると見なして、その等速直進運動区間毎に、時系列観測方位データを基にした時系列推定方位データと時系列観測周波数データを基にした時系列推定周波数データとを求め、時系列観測方位データと時系列推定方位データとの残差及び時系列観測周波数データと時系列推定周波数データとの残差の和を評価値として、変針変速の時刻を推定する。この推定した変針変速の時刻の前後の区間のそれぞれを前記目標体は等速直進運動していると見なして、この等速直進運動区間毎に、時系列観測方位データを基にした時系列推定方位データと時系列観測周波数データを基にした時系列推定周波数データとを求め、時系列観測方位データと時系列推定方位データとの残差及び時系列観測周波数データと時系列推定周波数データとの残差の和を評価値として、目標体の運動ベクトルを推定する。
【発明の効果】
【００１５】
本発明の目標運動解析方法及び装置により、変針変速を繰り返すジグザグ運動をする目標体であっても、目標体の運動ベクトルを安定して推定することができる。
【発明を実施するための最良の形態】
【００１６】
以下、本発明による目標運動解析方法及び装置の実施の形態について図面を参照して説明する。
【００１７】
図１は、目標運動解析方法を実現するための機能ブロック図である。変針変速検出部102は、観測情報データベース101と変針変速データベース103を参照し、新たな変針変速があったかを検出すると共に、変針変速時刻を推定する。目標運動解析部104は観測情報データベース101と変針変速データベース103を参照し、基準時刻における目標体位置と等速直進毎の針路、速力を推定する解析を行い、解析履歴データベース105に解析結果を格納する。変針変速データベース103上に変針変速が登録されていない場合には、基準時刻における目標位置と速力を推定する解析を行い、解析結果データベース105に解析結果を格納する。
【００１８】
図２は観測情報データベース部101を説明する図である。図に示すように、観測方位データベースと観測周波数データベースは独立に構成する。[]は時系列にデータを蓄積するために配列構造を示しており、蓄積されているデータ数を格納データ数として格納しておく。なお、観測方位誤差標準偏差及び観測周波数誤差標準偏差は観測した情報の正確さを示すものであり、必ずしも必要な情報ではないが、これを用いて解析することにより、より高精度の解析が可能となる。
【００１９】
図３は変針変速データベース103の構造を説明する図である。変針変速情報として推定変針変速時刻を格納するレコードを持ち、変針変速回数を保持するレコードを持つ構造となっている。
【００２０】
図４は変針変速検出部102の処理を説明するフローチャートある。401の処理では、変針変速データベース103から最新の変針変速時刻を取得し、観測情報データベース101中の最新変速時刻以降のデータ、すなわち最新の区間のデータを取得し（以降、解析用データと記述する）、以降の処理を実施する。ここで、最新の変針変速時刻とは、以前に推定した変針変速時刻の中で、最も新しい時刻データである。したがって、最新の区間のデータとは、最新の変針変速時刻から最も新しい観測時刻までの時系列の観測情報（解析用データ）である。
【００２１】
402の変針変速検出処理は解析用データから変針変速の有無を統計的な手段で検出する処理であり、多数の特許が出願(たとえば、特許文献2)されているため、ここでの説明を省略する。変針変速検出処理の結果を403にて判定し、変針変速が検出された場合には404以降の処理により変針変速時刻を推定する。
【００２２】
考え方は最新の観測データ時刻から最新の変針変速時刻まで、変針変速時刻を仮定値として2レグ解析を実施して、最も評価値の良い仮定値を推定変針変速時刻とするものである。2レグ解析の詳細については後述するが、２つの等速直進区間の時系列観測データを用いて、基準時刻における目標体の位置座標、変針変速前後の速力、固有周波数を推定し、評価値として正規化方位残差と正規化周波数残差の和を出力するものである。
【００２３】
まず、処理404において最良の評価値Φに初期値Φinitを設定する。評価値は正規化方位残差と正規化周波数残差であるため、小さい方が良い評価値である。したがって、Φinitは十分に大きな値とする。次に処理405において変針変速時刻の仮定値tzigに最新の解析用データの観測時刻tnewを設定する。処理406の2レグ解析は変針変速時刻を仮定値tzigとして基準時刻における目標体の位置座標、変針変速前後の速力、固有周波数を推定し、評価値φを出力する。次に分岐407において、最良の評価値Φと２レグ解析にて計算した評価値φを比較し、2レグ解析にて計算した評価値φの方が良ければ、処理408及び処理409にて最良の評価値Φと、最良の評価値となる推定変針変速時刻Tzigを更新する。処理410にて変針変速時刻の仮定値をdtだけ過去に遡り、分岐411にて変針変速最古時刻Tmaxより過去でなければ、処理406以降を繰り返す。最終的に変針変速時刻はTzigとなり、処理412にて変針変速データベース103に格納する。
【００２４】
変針変速最古時刻Tmaxは、最新の区間の中の最も古い時刻の意であり、前述の最新の変針変速時刻である。図４の処理が時間を遡るように繰り返すために、変針変速時刻の仮定値をdtだけ過去に遡るようにしているのであって、図４の処理が変針変速最古時刻Tmaxから最新の観測時刻に向けて処理するならば、時間の推移は逆になることは自明であろう。またdtは、推定する変針変速時刻に要求される分解能より小さな値であればよい。
【００２５】
図５は解析履歴データベース105の構造を説明する図である。解析結果は固有周波数とレグと呼ぶ等速直進区間の運動ベクトルで管理する。解析履歴データベースに格納しているレグ数は格納レグ数に保持する。レグには解析基準時刻、解析基準時刻における目標位置X座標、解析基準時刻におけるY座標、速力X成分、速力Y成分を格納する構造となっている。
【００２６】
図６は目標運動解析部104の処理を説明するフローチャートである。まず、分岐601において、変針変速データベース103を参照して、変針変速回数が０の場合は処理602の１レグ解析を実行し、変針変速回数が１以上の場合は603の２レグ解析を実行する。
【００２７】
ここで、1レグ解析は、前述の最新の変針変速時刻から最新の観測時刻に至る区間の、目標体の運動モデルとして等速直進モデルを適用するものであり、解析基準時刻における目標体の位置座標、速力、固有周波数を推定するものである。２レグ解析は目標体の運動モデルとして、前述の最新の変針変速時刻から推定した変針変速時刻までの区間と推定した変針変速時刻から最新の観測時刻までの区間とのそれぞれの区間に等速直進モデルを適用するものであり、解析基準時刻における目標体の位置座標、固有周波数、等速直進区間毎の針路、速力を推定するものである。
【００２８】
最後に、処理604にて解析結果を解析履歴データベース105に格納する。
【００２９】
図７は１レグ解析において推定する運動パラメータを説明する観測体と目標体の幾何学的関係を示した図である。701は観測体の運動軌跡、702は目標体の運動軌跡を表している。703〜705は目標体から放射される音波を観測した観測体の位置、705が最新の観測データを取得した観測体の位置を表しており、最新観測データを取得した時刻を解析基準時刻とする。１レグ解析において推定するパラメータは解析基準時刻における目標体位置X座標707、解析基準時刻における目標体位置Y座標708、速力X成分710、速力Y成分711及び目標体固有周波数である。
【００３０】
これらのパラメータは最小二乗法により式(1)を最小化することにより推定する。ここで、tは解析基準時刻を0として正規化した時刻を表し、Cは音速、co,moは観測体の針路、速力を表す。また、yi、ziは時刻tiにおける観測方位及び観測周波数であり、式(2)(3)はそれぞれ推定方位、推定周波数である。また、xoi、yoiは時刻tiにおける観測体の位置（xoi，yoi）である。
【００３１】
【数１】
0005203654.tif 000002
【００３２】
【数２】
0005203654.tif 000003
【００３３】
【数３】
0005203654.tif 000004
【００３４】
【数４】
0005203654.tif 000005
【００３５】
【数５】
0005203654.tif 000006
【００３６】
【数６】
0005203654.tif 000007
特に、式(2)(3)から構成される評価値(1)はパラメータx1,x2,x3,x4,x5に対して非線形であるため、反復により解く必要があり、例えば、Levenberg-Marquardt法ではΔXが十分小さくなるまで式(7)(8)の計算繰り返し実施する。なお、λは定数、Ｉは単位行列である。またbは、ヤコビアン行列と残差ベクトルの積であり、ここでは、方位ヤコビアン行列の転置行列と方位残差ベクトルの積、および周波数ヤコビアン行列の転置行列と周波数残差ベクトルの積の二つの積の和である。
【００３７】
【数７】
0005203654.tif 000008
【００３８】
【数８】
0005203654.tif 000009
式(8)のA、Bは方位情報及び周波数情報のヤコビアン行列であり、次式となる。
【００３９】
【数９】
0005203654.tif 000010
【００４０】
【数１０】
0005203654.tif 000011
【００４１】
【数１１】
0005203654.tif 000012
【００４２】
【数１２】
0005203654.tif 000013
【００４３】
【数１３】
0005203654.tif 000014
【００４４】
【数１４】
0005203654.tif 000015
【００４５】
【数１５】
0005203654.tif 000016
【００４６】
【数１６】
0005203654.tif 000017
【００４７】
【数１７】
0005203654.tif 000018
【００４８】
【数１８】
0005203654.tif 000019
【００４９】
【数１９】
0005203654.tif 000020
【００５０】
【数２０】
0005203654.tif 000021
図８はｎレグ解析において等速直進区間が２区間の場合に、推定する運動パラメータを説明する観測体と目標体の幾何学的関係を示した図である。801は観測体の運動軌跡、802は目標体の運動軌跡を表している。803は最古の観測体の位置座標、804は推定変針変速時刻における観測体の位置座標、805は最新の観測体の位置座標を表している。また806は、推定した変針変速時刻における目標体の位置であり、そのX座標を807で、そのY座標を808で表している。
【００５１】
解析基準時刻は最新の観測時刻とし、２レグ解析において推定するパラメータは解析基準時刻における目標体位置812のX座標813、解析基準時刻における目標体位置812のY座標814、変針変速前の速力809のX成分810、変針変速前の速力809のY成分811、変針変速後の速力815のX成分816、変針変速後の速力815のY成分817及び固有周波数である。
【００５２】
これらのパラメータは１レグ解析と同様に式(1)を最小化することにより推定する。なお、図４のフローチャートでは、２レグ解析における式(1)のS(X)をΦで表している。ただし、推定方位及び推定周波数は推定変針変速時刻をtz（図４のフローチャートでは、Tzigと表記している。）として、式(21)〜(24)となる。
【００５３】
【数２１】
0005203654.tif 000022
【００５４】
【数２２】
0005203654.tif 000023
【００５５】
【数２３】
0005203654.tif 000024
【００５６】
【数２４】
0005203654.tif 000025
【００５７】
【数２５】
0005203654.tif 000026
【００５８】
【数２６】
0005203654.tif 000027
【００５９】
【数２７】
0005203654.tif 000028
【００６０】
【数２８】
0005203654.tif 000029
【００６１】
【数２９】
0005203654.tif 000030
Levenberg-Marquardt法等により計算するためのヤコビアンは式(30)〜(56)となる。
【００６２】
【数３０】
0005203654.tif 000031
【００６３】
【数３１】
0005203654.tif 000032
【００６４】
【数３２】
0005203654.tif 000033
【００６５】
【数３３】
0005203654.tif 000034
【００６６】
【数３４】
0005203654.tif 000035
【００６７】
【数３５】
0005203654.tif 000036
【００６８】
【数３６】
0005203654.tif 000037
【００６９】
【数３７】
0005203654.tif 000038
【００７０】
【数３８】
0005203654.tif 000039
【００７１】
【数３９】
0005203654.tif 000040
【００７２】
【数４０】
0005203654.tif 000041
【００７３】
【数４１】
0005203654.tif 000042
【００７４】
【数４２】
0005203654.tif 000043
【００７５】
【数４３】
0005203654.tif 000044
【００７６】
【数４４】
0005203654.tif 000045
【００７７】
【数４５】
0005203654.tif 000046
【００７８】
【数４６】
0005203654.tif 000047
【００７９】
【数４７】
0005203654.tif 000048
【００８０】
【数４８】
0005203654.tif 000049
【００８１】
【数４９】
0005203654.tif 000050
【００８２】
【数５０】
0005203654.tif 000051
【００８３】
【数５１】
0005203654.tif 000052
【００８４】
【数５２】
0005203654.tif 000053
【００８５】
【数５３】
0005203654.tif 000054
【００８６】
【数５４】
0005203654.tif 000055
【００８７】
【数５５】
0005203654.tif 000056
【００８８】
【数５６】
0005203654.tif 000057
【実施例】
【００８９】
実施例として、本発明による目標運動解析装置の構成例を図９に示す。目標運動解析装置901は図９に示すように、CPU904、メモリ905、出力装置906、I/F(インターフェイス)907および908、入力装置909から構成されるような電子計算機で実現可能であり、音響センサ902及び航海センサ903に接続される。目標運動解析装置のサイズや演算速度に対する制約により、専用装置として構成されても良い。
【００９０】
音響センサ902は目標体から放射される音波を受波して、音波の到来方位及び周波数を目標運動解析装置901に入力する。航海センサ903は電磁ログ、ジャイロ、GPS等により、観測体の位置、速力、針路を目標運動解析装置901に入力する。入力された音波の到来方位、周波数、観測体の位置、速力、針路はメモリ905の観測情報データベースに蓄積される。変針変速検出部及び目標運動解析部は前述した方法に基づく処理を行い、その結果を変針変速データベース及び解析履歴データベースに格納する。
【００９１】
メモリ905上で保有する観測情報データベース、変針変速データベース、解析履歴データベースの内容をCRT等の出力装置906に表示し、オペレータに提示する。オペレータは出力装置906に出力された観測情報データベース、変針変速データベース、解析履歴データベースの内容を参照することにより、目標の変針変速を看破した場合には入力装置909により変針変速時刻を変針変速データベースに格納することも可能である。
【００９２】
例えば、図１０のような横軸に時刻、縦軸に観測方位とした図を出力装置906に表示する。図１０は目標の運動パターンを表１のとおりとした場合に音響センサで観測される音波の到来方位であり、図１０からも700秒に観測方位の傾向が変化していることがわかる。
【００９３】
【表１】
0005203654.tif 000058
しかしながら、実際に観測される方位は音波伝搬中の屈折の影響や音響センサにおける計測誤差もあり、図１０のように鮮明な方位の傾向変化を読み取る事は難しいため、本装置は統計的な手段により変針変速を看破する変針変速検出部を備えている。
【００９４】
図１３は本実施例による目標運動解析方法のシミュレーション結果であり、1301は観測体の航跡、1302は目標体の航跡、1303は解析周期30秒として目標運動解析を実施した時の目標体の推定位置座標履歴、1304は観測体と目標体の相対距離、1305は目標運動解析結果による観測体と目標体の推定相対距離、1306は目標体の速力、1307は目標運動解析結果による目標体の推定速力、1308は目標体の針路、1309は目標運動解析結果による目標体の推定針路を表している。シミュレーション開始後900秒で目標体は針路を180°から150°に変針しており、目標体の変針以降も解析結果を安定して出力している。なお、本シミュレーションの条件は図１１、１２と同様であり、本実施例により目標体が変針変速している場合でも、安定した解析を実施することができていることが分かる。
【図面の簡単な説明】
【００９５】
【図１】目標運動解析を実施するためのブロック図である。
【図２】観測情報データベース部の構造を説明する図である。
【図３】変針変速データベース部の構造を説明する図である。
【図４】変針変速検出部のフローチャート図である。
【図５】解析履歴データベース部の構造を説明する図である。
【図６】目標運動解析部のフローチャート図である。
【図７】等速直進モデルの目標運動解析における観測体と目標体の幾何学的関係を示した図である。
【図８】変針変速モデルの目標運動解析における観測体と目標体の幾何学的関係を示した図である。
【図９】目標運動解析装置の構成例である。
【図１０】目標が変針した場合の音波の到来方位を時系列に示したグラフである。
【図１１】変針する目標体に対して、従来技術である等速直進モデルによる目標運動解析を適用した場合のシミュレーション結果である。
【図１２】変針する目標体に対して、従来技術である変針変速検出後の観測データを用いて目標運動解析を適用した場合のシミュレーション結果である。
【図１３】変針する目標体に対して、目標運動解析を適用した場合のシミュレーション結果である。
【符号の説明】
【００９６】
101…観測情報データベース部、102…変針変速検出部、103…変針変速データベース部、104…目標運動解析部、105…解析履歴データベース部、901…目標運動解析装置、902…音響センサ、903…航海センサ、904…ＣＰＵ、905…メモリ、906…出力装置、907〜908…インターフェース、909…入力装置。
【図１】
0005203654.tif 000059
【図２】
0005203654.tif 000060
【図３】
0005203654.tif 000061
【図４】
0005203654.tif 000062
【図５】
0005203654.tif 000063
【図６】
0005203654.tif 000064
【図７】
0005203654.tif 000065
【図８】
0005203654.tif 000066
【図９】
0005203654.tif 000067
【図１０】
0005203654.tif 000068
【図１１】
0005203654.tif 000069
【図１２】
0005203654.tif 000070
【図１３】
0005203654.tif 000071