(57)【要約】【課題】１つの信号に対する複数のピークの発生を抑制し、信号の到来方位を容易に特定できる整相器および整相処理方法を提供する。【解決手段】整相器は、複数のセンサに到来する信号のセンサ出力に基づく共分散行列から、空間を信号部分空間および雑音部分空間に分離し、各空間のそれぞれに属する固有ベクトルを取得し、雑音部分空間の固有ベクトルから推定したセンサ出力の到来方位と、信号部分空間の固有ベクトルとから、固有ベクトルに対応付ける方位を決定し、決定された方位から、この方位に対応する固有ベクトルの成分を除外する除外重みを算出し、信号部分空間の固有ベクトルと、算出した除外重みとから、各方位に対応する適応重みを算出し、信号部分空間の固有ベクトルに対応付ける方位を決定する際に、この方位の周囲の複数の方位を、信号部分空間の固有ベクトルに対応付け、複数の方位について、除外重みを設定する。【選択図】図４
2018189602.tif 000002
【特許請求の範囲】
【請求項１】
アレイを構成する複数のセンサに到来する信号を示すセンサ出力信号に基づき、前記アレイの周囲に存在する信号源に関する情報を取得する整相器であって、
前記センサ出力信号に基づき、共分散行列を算出する共分散行列算出部と、
前記共分散行列に基づき、空間を信号部分空間および雑音部分空間に分離するとともに、前記信号部分空間および前記雑音部分空間に属する固有ベクトルを取得する固有値分解処理部と、
前記雑音部分空間に属する固有ベクトルに基づき、前記センサ出力信号の到来方位を推定する方位推定処理部と、
推定された前記到来方位と、前記信号部分空間に属する固有ベクトルとに基づき、該固有ベクトルに対応付ける到来方位を示す方位インデックスを決定する方位インデックス決定部と、
決定された前記方位インデックスに基づき、該方位インデックスに対応付けられた前記固有ベクトルの成分を除外する除外重みを算出する除外重み算出部と、
前記信号部分空間に属する固有ベクトルと、算出された前記除外重みとに基づき、各方位に対応する適応重みを算出する重み算出部と
を備え、
前記方位インデックス決定部は、
前記信号部分空間に属する固有ベクトルに対応付ける前記方位インデックスを決定する際に、決定された前記方位インデックスの周囲の複数の方位を示す方位インデックスを、前記信号部分空間に属する前記固有ベクトルに対応付け、
前記除外重み算出部は、
決定された前記方位インデックスに対応付けられた前記複数の方位について、除外重みを設定する
ことを特徴とする整相器。
【請求項２】
前記除外重み算出部は、
前記複数の方位について設定された前記除外重みを、前記方位インデックスの関数を用いて算出する
ことを特徴とする請求項１に記載の整相器。
【請求項３】
決定された前記方位インデックスの時間的な平均および分散を推定する推定平均分散算出部をさらに備え、
前記除外重み算出部は、
推定された前記方位インデックスの平均および分散に基づき、前記方位インデックスの関数を変化させる
ことを特徴とする請求項２に記載の整相器。
【請求項４】
前記方位推定処理部は、
ＭＵＳＩＣ処理によって前記センサ出力信号の到来方位を推定し、
前記方位インデックス決定部は、
前記ＭＵＳＩＣ処理によって得られた推定結果と、前記信号部分空間に属する固有ベクトルおよび方位を示すステアリングベクトルの内積の結果とを乗算し、
前記乗算の結果を示す乗算値の大きさが最大となる方位を示す方位インデックスを、前記信号部分空間に属する前記固有ベクトルに対応付ける
ことを特徴とする請求項１〜３のいずれか一項に記載の整相器。
【請求項５】
アレイを構成する複数のセンサに到来する信号を示すセンサ出力信号に基づき、前記アレイの周囲に存在する信号源に関する情報を取得する整相処理方法であって、
前記センサ出力信号に基づき、共分散行列を算出する共分散行列算出ステップと、
前記共分散行列に基づき、空間を信号部分空間および雑音部分空間に分離するとともに、前記信号部分空間および前記雑音部分空間に属する固有ベクトルを取得する固有値分解ステップと、
前記雑音部分空間に属する固有ベクトルに基づき、前記センサ出力信号の到来方位を推定する方位推定ステップと、
推定された前記到来方位と、前記信号部分空間に属する固有ベクトルとに基づき、該固有ベクトルに対応付ける到来方位を示す方位インデックスを決定する方位インデックス決定ステップと、
決定された前記方位インデックスに基づき、該方位インデックスに対応付けられた前記固有ベクトルの成分を除外する除外重みを算出する除外重み算出ステップと、
前記信号部分空間に属する固有ベクトルと、算出された前記除外重みとに基づき、各方位に対応する適応重みを算出する重み算出ステップと
を有し、
前記方位インデックス決定ステップは、
前記信号部分空間に属する固有ベクトルに対応付ける前記方位インデックスを決定する際に、決定された前記方位インデックスの周囲の複数の方位を示す方位インデックスを、前記信号部分空間に属する前記固有ベクトルに対応付け、
前記除外重み算出ステップは、
決定された前記方位インデックスに対応付けられた前記複数の方位について、除外重みを設定する
ことを特徴とする整相処理方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、複数のセンサで構成されたアレイによって受信した信号を処理する整相器および整相処理方法に関する。
【背景技術】
【０００２】
従来、アレイを構成する複数のセンサからのセンサ出力を用いたアレイ信号処理は、ソーナー、レーダー、および地震源の推定等の種々の分野で使用されている。このような分野においては、信号源の位置を推定するために、空間的に異なる位置に複数のセンサを設置し、センサ出力間に生じる時間差情報または位相差情報に基づき、周囲に存在する信号源に関する情報を取得している。
【０００３】
アレイ信号処理の代表的な例としては、例えば整相処理が挙げられる。整相処理では、センサ出力に対して遅延時間を設定し、その結果に対して重み付き加算を行うことにより、特定の方位から到来した信号を抽出する。整相処理は、整相処理の重み付き加算における重みの与え方に応じて、例えば、重みが固定的な「従来整相処理」と、周囲の環境に応じて重みを動的に変化させる「適応整相処理」との２つの処理に大別して分類することができる。
【０００４】
ここで、ＤＩＦＡＲ（ＤＩｒｅｃｔｉｏｎａｌ  Ｆｒｅｑｕｅｎｃｙ  Ａｎａｌｙｓｉｓ  ａｎｄ  Ｒｅｃｏｒｄｉｎｇ）を対象とした適応整相処理の１つに、ＭＶＤＲ（Ｍｉｎｉｍｕｍ  Ｖａｒｉａｎｃｅ  Ｄｉｓｔｏｒｔｉｏｎｌｅｓｓ  Ｒｅｓｐｏｎｓｅ）方式と称する処理方式が提案されている（例えば、非特許文献１参照）。
【０００５】
しかしながら、適応整相処理としてＭＶＤＲ方式を用いた場合には、整相方位と信号の到来方位とが一致しないと、ビーム応答のピークレベルが劣化してしまうという問題があった。そこで、最近では、このような問題を解決するために、ＥＢＡＥ（Ｅｉｇｅｎｖｅｃｔｏｒ／Ｂｅａｍ  Ａｓｓｏｃｉａｔｉｏｎ  ａｎｄ  Ｅｘｃｉｓｉｏｎ）方式と称する処理方式が提案されている。
【０００６】
ＥＢＡＥ方式は、センサ出力の共分散行列の固有値分解によって得られる固有ベクトルと、方位を示すとともに整相方位を決定するステアリングベクトルとを対応付ける。そして、このステアリングベクトルを用いて整相処理を行う際に、対応する固有ベクトルの情報を除外して適応的な重み計算を行うものである。
【０００７】
このようなＥＢＡＥ方式は、固有ベクトルとステアリングベクトルとを１対１で対応付けるとともに、適応重みを計算する際に固有ベクトルの情報を除外することにより、高い方位分解能を有する。また、様々な雑音の影響を受けるような環境であっても、センサ出力に対してロバストな整相結果を得ることができる。
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】D.Desrochers and R.F.Marsden,“High Resolution Beamforming Techniques Applied To a DIFAR Sonobuoy”, Canadian Acoustics, 1999
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかしながら、上述したＥＢＡＥ方式をＤＩＦＡＲに適用した場合には、信号の到来方位付近に複数のピークが発生してしまうことがある。そのため、信号の到来方位を特定することが困難であるという課題があった。
【００１０】
そこで、適応整相処理としてＥＢＡＥ方式を適用した場合でも、１つの信号に対する複数のピークの発生を抑制し、信号の到来方位を容易に特定することができる整相器および整相処理方法が望まれている。
【課題を解決するための手段】
【００１１】
本発明に係る整相器は、アレイを構成する複数のセンサに到来する信号を示すセンサ出力信号に基づき、アレイの周囲に存在する信号源に関する情報を取得する整相器であって、センサ出力信号に基づき、共分散行列を算出する共分散行列算出部と、共分散行列に基づき、空間を信号部分空間および雑音部分空間に分離するとともに、信号部分空間および雑音部分空間に属する固有ベクトルを取得する固有値分解処理部と、雑音部分空間に属する固有ベクトルに基づき、センサ出力信号の到来方位を推定する方位推定処理部と、推定された到来方位と、信号部分空間に属する固有ベクトルとに基づき、固有ベクトルに対応付ける到来方位を示す方位インデックスを決定する方位インデックス決定部と、決定された方位インデックスに基づき、方位インデックスに対応付けられた固有ベクトルの成分を除外する除外重みを算出する除外重み算出部と、信号部分空間に属する固有ベクトルと、算出された除外重みとに基づき、各方位に対応する適応重みを算出する重み算出部とを備え、方位インデックス決定部は、信号部分空間に属する固有ベクトルに対応付ける方位インデックスを決定する際に、決定された方位インデックスの周囲の複数の方位を示す方位インデックスを、信号部分空間に属する固有ベクトルに対応付け、除外重み算出部は、決定された方位インデックスに対応付けられた複数の方位について、除外重みを設定するものである。
【００１２】
また、本発明に係る整相処理方法は、アレイを構成する複数のセンサに到来する信号を示すセンサ出力信号に基づき、アレイの周囲に存在する信号源に関する情報を取得する整相処理方法であって、センサ出力信号に基づき、共分散行列を算出する共分散行列算出ステップと、共分散行列に基づき、空間を信号部分空間および雑音部分空間に分離するとともに、信号部分空間および雑音部分空間に属する固有ベクトルを取得する固有値分解ステップと、雑音部分空間に属する固有ベクトルに基づき、センサ出力信号の到来方位を推定する方位推定ステップと、推定された到来方位と、信号部分空間に属する固有ベクトルとに基づき、固有ベクトルに対応付ける到来方位を示す方位インデックスを決定する方位インデックス決定ステップと、決定された方位インデックスに基づき、方位インデックスに対応付けられた固有ベクトルの成分を除外する除外重みを算出する除外重み算出ステップと、信号部分空間に属する固有ベクトルと、算出された除外重みとに基づき、各方位に対応する適応重みを算出する重み算出ステップとを有し、方位インデックス決定ステップは、信号部分空間に属する固有ベクトルに対応付ける方位インデックスを決定する際に、決定された方位インデックスの周囲の複数の方位を示す方位インデックスを、信号部分空間に属する固有ベクトルに対応付け、除外重み算出ステップは、決定された方位インデックスに対応付けられた複数の方位について、除外重みを設定するものである。
【発明の効果】
【００１３】
以上のように、本発明によれば、信号部分空間の固有ベクトルに対応付ける方位インデックスを決定する際に、決定された方位インデックスの周囲の複数の方位を示す方位インデックスを、信号部分空間の固有ベクトルに対応付け、複数の方位について除外重みを設定することにより、１つの信号に対する複数のピークの発生を抑制し、信号の到来方位を容易に特定することができる。
【図面の簡単な説明】
【００１４】
【図１】ＤＩＦＡＲの指向性と信号の到来方位との関係について説明するための概略図である。
【図２】従来のＥＢＡＥ方式を用いた整相器の構成の一例を示すブロック図である。
【図３】従来型ＥＢＡＥ方式による問題点について説明するための概略図である。
【図４】実施の形態１に係る改善型ＥＢＡＥ方式を用いた整相器の構成の一例を示すブロック図である。
【図５】実施の形態１における固有ベクトル除外重みの形状の一例を示す概略図である。
【図６】実施の形態１に係る整相器による整相処理結果について説明するための概略図である。
【図７】実施の形態２に係る改善型ＥＢＡＥ方式を用いた整相器の構成の一例を示すブロック図である。
【図８】実施の形態２における信号到来方位の推定分散と定数との関係の一例を示すグラフである。
【図９】実施の形態２における固有ベクトル除外重みの形状の一例を示す概略図である。
【図１０】実施の形態２に係る整相器による整相処理結果について説明するための概略図である。
【発明を実施するための形態】
【００１５】
実施の形態１．
以下、本発明の実施の形態１に係る整相器について説明する。本実施の形態１に係る整相器では、ＤＩＦＡＲに対して、従来のＥＢＡＥ方式を改善した改善型ＥＢＡＥ方式による適応整相処理を適用する。以下では、まず、本実施の形態１に係る整相器について説明する前に、本発明の前提となる「ＤＩＦＡＲに対する適応整相処理の適用」と、「従来型ＥＢＡＥ方式を用いた適応整相処理」とについて説明する。
【００１６】
（ＤＩＦＡＲについて）
図１は、ＤＩＦＡＲの指向性と信号の到来方位との関係について説明するための概略図である。ＤＩＦＡＲは、無指向性センサ、正弦指向性センサ、および余弦指向性センサの３つの音響センサによって構成されている。図１に示すように、正弦指向性センサおよび余弦指向性センサは、指向性の最大応答軸が互いに直交となるように配置されている。
【００１７】
このように正弦指向性センサおよび余弦指向性センサが配置された場合において、北の方角を示す真方位軸と余弦指向性センサの最大応答軸とのなす角を角度αとし、余弦指向性センサの最大応答軸と信号の到来方位とのなす角を角度φとする。また、真方位軸を基準とした場合の信号の到来方位を角度θとする。
【００１８】
このとき、無指向性センサ出力の時系列ｓ０［ｎ］、正弦指向性センサ出力の時系列ｓｓ［ｎ］、および余弦指向性センサ出力の時系列ｓｃ［ｎ］は、式（１）で表される。式（１）において、記号「ｎ」は、時刻を表すインデックスを示す。
【００１９】
【数１】
2018189602.tif 000003
【００２０】
また、式（１）で示される各センサ出力の時系列ｓ０［ｎ］、ｓｓ［ｎ］およびｓｃ［ｎ］に対して離散フーリエ変換を適用した場合の、ｋ番目の周波数の各センサ出力ｘ０［ｋ］、ｘｓ［ｋ］およびｘｃ［ｋ］は、式（２）および（３）で表される。
【００２１】
【数２】
2018189602.tif 000004
【００２２】
【数３】
2018189602.tif 000005
【００２３】
（ＤＩＦＡＲへの適応整相処理の適用）
次に、ＤＩＦＡＲに対する適応整相処理の適用について説明する。ここでは、適応整相処理の１つであるＭＶＤＲ方式をＤＩＦＡＲに対して適用する場合について説明する。
【００２４】
ＭＶＤＲ方式において、整相方位を決定するステアリングベクトルｖ［ｐ，ｋ］は、式（４）で表される。式（４）において、「ｐ」は、全周において整相する方位を表すインデックス（以下、「方位インデックス」と適宜称する）を示す。「φｐ」は、方位インデックスｐに対応した整相方位を示す。
【００２５】
【数４】
2018189602.tif 000006
【００２６】
また、ＭＶＤＲ方式のビーム応答Ｅ＾ＭＶＤＲは、ステアリングベクトルｖ［ｐ，ｋ］を用いて、式（５）によって表される。なお、以下の説明においては、数式中で「Ｅ」などの文字の上に「＾」がつけられた文字を「Ｅ＾」と記載するものとする。式（５）において、「［・］Ｔ」は、転置の操作を示し、「Ｒ［ｋ］−１」は、３つのセンサ出力の共分散行列Ｒ［ｋ］の逆行列を示す。
【００２７】
【数５】
2018189602.tif 000007
【００２８】
ここで、センサ出力の共分散行列Ｒ［ｋ］は、ｋ番目の周波数のセンサ出力の母平均または母分散によって定義されるため、一般的には、直接的に求めることが困難である。そのため、例えば、時刻「ｊ−Ｂ＋１」から時刻「ｊ」までのＢ個の瞬時の共分散行列を平均することによって、時刻「ｊ」における共分散行列の推定値Ｒ＾ｊ［ｋ］を求める。そして、求めた共分散行列の推定値Ｒ＾ｊ［ｋ］を、式（５）において、共分散行列Ｒ［ｋ］の代わりに使用する。
【００２９】
このときの共分散行列の推定値Ｒ＾ｊ［ｋ］は、式（６）により算出することができる。なお、式（６）における「［・］Ｈ」は、共役転置の操作を示す。また、以下の説明では、共分散行列の推定値Ｒ＾ｊ［ｋ］を、簡単のため「Ｒ＾［ｋ］」と表記する。
【００３０】
【数６】
2018189602.tif 000008
【００３１】
以上により、ＭＶＤＲ方式による適応整相処理をＤＩＦＡＲに適用することができる。なお、適応整相処理としてＭＶＤＲ方式を用いた場合には、整相方位と信号の到来方位とが一致しないなどの問題が発生すると、ビーム応答のピークのレベルが劣化する。
【００３２】
（ＥＢＡＥ方式について）
次に、従来のＥＢＡＥ方式による適応整相処理（以下、「従来型ＥＢＡＥ方式」と適宜称する）について説明する。従来型ＥＢＡＥ方式は、上述したＭＶＤＲ方式において生じるビーム応答のピークのレベルが劣化するという問題を解決するために提案されたものである。
【００３３】
従来型ＥＢＡＥ方式は、概略的には、センサ出力の共分散行列の固有値分解によって得られる固有ベクトルと、整相方位を決定するステアリングベクトルとを対応付ける。そして、このステアリングベクトルを用いて整相処理を行う際に、対応する固有ベクトルの情報を除外して重みを算出する適応整相処理である。
【００３４】
図２は、従来のＥＢＡＥ方式を用いた整相器１００の構成の一例を示すブロック図である。図２に示すように、整相器１００は、共分散行列算出部１０１、固有値分解処理部１０２、ＭＵＳＩＣ（ＭＵｌｔｉｐｌｅ  ＳＩｇｎａｌ  Ｃｌａｓｓｉｆｉｃａｔｉｏｎ）処理部１０３、方位インデックス決定部１０４、固有ベクトル除外重み算出部１０５、適応重み算出部１０６、整相処理部１０７、および積分処理部１０８で構成されている。整相器１００には、アレイを構成する複数のセンサによって得られる、周波数ｋの出力を並べたベクトルであるセンサ出力ｘ［ｋ］が入力される。
【００３５】
共分散行列算出部１０１は、整相器１００に入力されたセンサ出力ｘ［ｋ］に基づき、式（６）を用いて共分散行列の推定値Ｒ＾［ｋ］を算出する。
【００３６】
固有値分解処理部１０２は、共分散行列算出部１０１で算出された共分散行列の推定値Ｒ＾［ｋ］に対して固有値分解処理を行う。固有値分解処理では、算出された共分散行列の推定値Ｒ＾［ｋ］を、式（７）に示すように、行列の加算形式に変換する。式（７）において、「Ｍ」は、アレイを構成するセンサ数を示す。また、「λｍ［ｋ］」は、周波数ｋにおけるｍ番目の固有値を示し、「ｑｍ［ｋ］」は、周波数ｋにおけるｍ番目の固有値に対応する固有ベクトルを示す。
【００３７】
【数７】
2018189602.tif 000009
【００３８】
このように変換された共分散行列の推定値Ｒ＾［ｋ］は、信号部分空間および雑音部分空間と称する２つの空間に分離することができる。信号部分空間は、Ｍ個の固有値のうち、Ｄ個の大きな固有値に対応する固有ベクトルによって形成される空間である。雑音部分空間は、Ｍ個の固有値のうち、残りの固有値に対応する固有ベクトルによって形成される空間である。
【００３９】
したがって、式（７）は、分離された信号部分空間および雑音部分空間の２つの空間を考慮することにより、式（８）のように表すことができる。式（８）における右辺の第１項は、信号部分空間の固有値および固有ベクトルの集まりを示し、第２項は、雑音部分空間の固有値および固有ベクトルの集まりを示す。固有値λｍ［ｋ］は、信号または雑音のパワーに関する情報を含んでいる。また、信号部分空間の固有ベクトルｑｍ［ｋ］は、信号の到来方位に関する情報を含んでいる。
【００４０】
【数８】
2018189602.tif 000010
【００４１】
このように、固有値分解処理部１０２は、共分散行列算出部１０１で算出された共分散行列の推定値Ｒ＾［ｋ］を、２つの部分空間に分離し、信号部分空間および雑音部分空間それぞれの固有値λｍ［ｋ］および固有ベクトルｑｍ［ｋ］を取得する。
【００４２】
なお、信号部分空間に属する固有ベクトルの数Ｄを推定する方法としては、例えば、ＡＩＣ（Ａｋａｉｋｅ  ｉｎｆｏｒｍａｔｉｏｎ  ｃｒｉｔｅｒｉｏｎ）、ＭＤＬ（Ｍｉｎｉｍｕｍ  ｄｅｓｃｒｉｐｔｉｏｎ  ｌｅｎｇｔｈ）等の方法がある。これらの方法は、信号部分空間に属する固有値の数Ｄを推定する方法として周知の方法であり、本発明における特徴を示すものではないため、ここでは説明を省略する。
【００４３】
ＭＵＳＩＣ処理部１０３は、固有値分解処理部１０２で得られた雑音部分空間の固有ベクトルに基づき、信号の到来方位を推定する処理の一例であるＭＵＳＩＣ処理を行い、ＭＵＳＩＣスペクトルを算出する。ＭＵＳＩＣスペクトルＭ［ｐ，ｋ］は、式（９）に基づき算出される。このようなＭＵＳＩＣスペクトルＭ［ｐ，ｋ］は、信号の到来方位に関する情報を含んでおり、高い方位分解能を有している。
【００４４】
【数９】
2018189602.tif 000011
【００４５】
方位インデックス決定部１０４は、信号部分空間に属する固有ベクトルｑｍ［ｋ］に対応付ける方位インデックスｐｍを決定する。具体的には、方位インデックス決定部１０４は、信号部分空間に属する固有ベクトルｑｍ［ｋ］と、ＭＵＳＩＣスペクトルＭ［ｐ，ｋ］とに基づき、信号部分空間に属する固有ベクトルｑｍ［ｋ］と方位インデックスｐとの対応付けを行う。
【００４６】
この場合には、まず、信号部分空間に属するｍ番目の固有ベクトルｑｍ［ｋ］とすべての方位のステアリングベクトルｖ［ｐ，ｋ］との内積ρｍ［ｐ，ｋ］を、式（１０）に基づき算出する。
【００４７】
【数１０】
2018189602.tif 000012
【００４８】
次に、式（１０）によって算出された内積ρｍ［ｐ，ｋ］とＭＵＳＩＣスペクトルＭ［ｐ，ｋ］との乗算結果が最大となる方位インデックスｐを、式（１１）に基づき、固有ベクトルｑｍ［ｋ］に対応付ける方位インデックスｐｍとして算出する。
【００４９】
【数１１】
2018189602.tif 000013
【００５０】
固有ベクトル除外重み算出部１０５は、方位インデックス決定部１０４で決定された方位インデックスｐｍに基づき、「ｍ」および「ｐ」の関数である固有ベクトル除外重みδ［ｍ，ｐ］を導出する。固有ベクトル除外重みδ［ｍ，ｐ］は、式（１２）によって定義される。
【００５１】
【数１２】
2018189602.tif 000014
【００５２】
適応重み算出部１０６は、固有ベクトル除外重みδ［ｍ，ｐ］に基づき、各方位に対応する適応重みｗ［ｐ，ｋ］を算出する。信号部分空間に属する固有ベクトルｑｍ［ｋ］に対応付けられた方位インデックスｐｍを用いて適応重みを算出する場合には、対応する固有ベクトルｑｍ［ｋ］の情報を除外して適応重みを算出する。
【００５３】
すなわち、この場合の適応重みｗ［ｐ，ｋ］は、式（１３）に基づき算出される。式（１３）において、βｍは、ｍ番目の固有値λｍ［ｋ］に対応した信号のＳＮ（Ｓｉｇｎａｌ／Ｎｏｉｓｅ）比を示す。
【００５４】
【数１３】
2018189602.tif 000015
【００５５】
整相処理部１０７は、センサ出力ｘ［ｋ］および適応重みｗ［ｐ，ｋ］に基づき、ビーム応答ｙ［ｐ，ｋ］を算出する。ビーム応答ｙ［ｐ，ｋ］は、式（１４）に基づき算出される。
【００５６】
【数１４】
2018189602.tif 000016
【００５７】
積分処理部１０８は、ビーム応答ｙ［ｐ，ｋ］の絶対値を２乗した後に時間積分を行う。
【００５８】
（ＤＩＦＡＲを対象とした従来型ＥＢＡＥ方式の適用の問題点）
ところで、従来型ＥＢＡＥ方式においてＭＵＳＩＣスペクトルを算出する場合には、上述した式（９）の分母に示すように、雑音部分空間に属する固有ベクトルとすべてのステアリングベクトルとの内積の結果を加算する。これにより、従来型ＥＢＡＥ方式では、実際に信号が到来する方位に対してピークの方位が安定するようにしている。
【００５９】
しかしながら、ＤＩＦＡＲにおいては、アレイを構成するセンサ数Ｍが「３」であるため、信号部分空間に属する固有ベクトル数Ｄを「１」もしくは「２」とした場合、雑音部分空間に属する固有ベクトル数Ｄは、「２」もしくは「１」である。したがって、式（９）の分母における固有ベクトルを用いた結果の加算数「Ｍ−（Ｄ＋１）」は、多くとも「２」であり、ピークの方位を安定させるには、加算数が極端に不足してしまう。そのため、ＤＩＦＡＲにおいては、ＭＵＳＩＣスペクトルを算出する場合に、信号のピークの方位が安定しない。
【００６０】
図３は、従来型ＥＢＡＥ方式による問題点について説明するための概略図である。図３（ａ）は、時間変化に伴ってＭＵＳＩＣスペクトルのピークの方位が変化する様子を示す。図３（ｂ）は、時間変化に伴ってビーム応答ピークの方位が変化する様子を示す。図３（ｃ）は、図３（ｂ）のビーム応答を時間積分した結果を示す。ここでは、２つの信号＃１および信号＃２の到来方位を特定する場合についての一例を示す。この例において、点線で示す方位が実際の信号＃１および信号＃２の到来方位であり、ピークの方位が特定した信号＃１および信号＃２の到来方位である。
【００６１】
図３（ａ）に示すように、信号のピークの方位が安定しないことにより、ＭＵＳＩＣスペクトルのピークの方位は、時間変化に伴ってばらつくことになる。その結果、方位インデックス決定部１０４において、信号部分空間に属する固有ベクトルｑｍ［ｋ］に対応付けられる方位インデックスｐｍは、時間変化に伴ってばらつくこととなってしまう。
【００６２】
また、これにより、図３（ｂ）に示すように、整相処理によって得られるビーム応答ｙ［ｐ，ｋ］のピークもばらついてしまう。そして、各時刻におけるビーム応答ｙ［ｐ，ｋ］のピーク方位がばらついているため、このようにばらついたビーム応答ｙ［ｐ，ｋ］を時間積分した場合には、図３（ｃ）に示すように、信号の到来方位付近に複数のピークが発生する。
【００６３】
このように、ＤＩＦＡＲに対して従来型ＥＢＡＥ方式を用いた場合には、信号の到来方位付近に複数のビーム応答のピークが発生するため、信号の到来方位を特定することが困難となる。
【００６４】
そこで、本実施の形態１では、従来型ＥＢＡＥ方式を改善した改善型ＥＢＡＥ方式を適応整相処理として用い、信号の到来方位を特定できるようにした。具体的には、本実施の形態１では、信号部分空間に属する固有ベクトルと方位インデックスとの対応付けを行う際に、信号部分空間に属する固有ベクトルｑｍに対応した方位インデックスｐｍを中心に前後複数の方位インデックスを対応付ける。これにより、ビーム応答の幅が広がり、ビーム応答を時間積分した際に、ビーム応答になだらかな１つのピークを生じさせることができる。
【００６５】
［整相器の構成］
図４は、本実施の形態１に係る改善型ＥＢＡＥ方式を用いた整相器１０の構成の一例を示すブロック図である。図４に示すように、整相器１０は、共分散行列算出部１、固有値分解処理部２、ＭＵＳＩＣ処理部３、方位インデックス決定部４、固有ベクトル除外重み算出部５、適応重み算出部６、整相処理部７、および積分処理部８で構成されている。
【００６６】
共分散行列算出部１は、アレイを構成する複数のセンサによって得られる、周波数ｋのセンサ出力ｘ［ｋ］が入力される。共分散行列算出部１は、入力されたセンサ出力ｘ［ｋ］に基づき、式（６）を用いて共分散行列の推定値Ｒ＾［ｋ］を算出する。
【００６７】
固有値分解処理部２は、共分散行列算出部１で算出された共分散行列の推定値Ｒ＾［ｋ］が入力される。固有値分解処理部２は、入力された共分散行列の推定値Ｒ＾［ｋ］に対して固有値分解処理を行い、式（７）および式（８）を用いて、信号部分空間および雑音部分空間それぞれの固有値λｍ［ｋ］および固有ベクトルｑｍ［ｋ］を取得する。
【００６８】
ＭＵＳＩＣ処理部３は、センサ出力ｘ［ｋ］の到来方位を推定する方位推定処理部であり、固有値分解処理部２で得られた雑音部分空間の固有ベクトルｑｍ［ｋ］が入力される。ＭＵＳＩＣ処理部３は、入力された雑音部分空間の固有ベクトルｑｍ［ｋ］に基づきＭＵＳＩＣ処理を行い、式（９）を用いてＭＵＳＩＣスペクトルＭ［ｐ，ｋ］を算出する。
【００６９】
方位インデックス決定部４は、固有値分解処理部２で得られた信号部分空間の固有ベクトルｑｍ［ｋ］と、ＭＵＳＩＣ処理部３で算出されたＭＵＳＩＣスペクトルＭ［ｐ，ｋ］とが入力される。方位インデックス決定部４は、入力された信号部分空間の固有ベクトルｑｍ［ｋ］およびＭＵＳＩＣスペクトルＭ［ｐ，ｋ］に基づき、信号部分空間に属する固有ベクトルｑｍ［ｋ］に対応する方位インデックスｐｍを決定する。
【００７０】
また、方位インデックス決定部４は、決定した方位インデックスｐｍを中心として周囲の方位インデックスを、信号部分空間に属する固有ベクトルｑｍ［ｋ］に同時に対応付ける。具体的には、対応付けられる方位インデックスは、前後複数の方位インデックス「ｐｍ−（Ｕ−１）／２，・・・，ｐｍ−（Ｕ−１）／２」である。ここで、記号「Ｕ」は、１つの固有ベクトルに対して対応付ける方位インデックスの数であり、奇数とする。
【００７１】
固有ベクトル除外重み算出部５は、方位インデックス決定部４で決定された方位インデックスｐｍが入力される。固有ベクトル除外重み算出部５は、入力された方位インデックスｐｍに基づき、「ｍ」および「ｐ」の関数である固有ベクトル除外重みδ［ｍ，ｐ］を導出する。このときの固有ベクトル除外重みδ［ｍ，ｐ］は、式（１５）によって定義される。
【００７２】
【数１５】
2018189602.tif 000017
【００７３】
ここで、改善型ＥＢＡＥ方式におけるビーム応答のピーク付近の形状は、式（１５）における値「ξ」の形状によって決定される。本実施の形態１では、固定的な固有ベクトル除外重みの形状が、例えば下に凸の２次関数の形状となるように、ξの値を決定する。これは、一般的に、整相処理のビーム応答のピーク付近の形状が、ＥＢＡＥ方式に限らず上に凸の２次関数で近似できるような形状であるためである。そして、このビーム応答のピーク付近の形状を決定するためには、整相処理のビーム応答の逆の特性となる固有ベクトル除外重みの形状が必要となるからである。したがって、固有ベクトル除外重みの形状が下に凸の２次関数となるようにξの値を決定することにより、ビーム応答のピーク付近の形状を尤もらしい形状に保つことができる。
【００７４】
図５は、本実施の形態１における固有ベクトル除外重みの形状の一例を示す概略図である。ここでは、図５に示すような固有ベクトル除外重みδ［ｍ，ｐ］を与える。また、この場合のξの値は、２次関数である「ξ＝ａｍ（ｐ−ｐｍ）２」を満たす値とする。値「ａｍ」は、固有ベクトル除外重みの２次関数としての形状を決定する定数である。
【００７５】
説明は図４に戻り、適応重み算出部６は、固有ベクトル除外重み算出部５で得られた固有ベクトル除外重みδ［ｍ，ｐ］と、固有値分解処理部２で得られた信号部分空間に属する固有値および固有ベクトル、ならびに雑音部分空間に属する固有値と、方位インデックス決定部４で決定された方位インデックスｐｍとが入力される。適応重み算出部６は、入力されたこれらの情報に基づき、式（１３）を用いて適応重みｗ［ｐ，ｋ］を算出する。
【００７６】
整相処理部７は、センサ出力ｘ［ｋ］と、適応重み算出部６で算出された適応重みｗ［ｐ，ｋ］が入力される。整相処理部７は、入力されたセンサ出力ｘ［ｋ］および適応重みｗ［ｐ，ｋ］に基づき、式（１４）を用いてビーム応答ｙ［ｐ，ｋ］を算出する。
【００７７】
積分処理部８は、整相処理部７で算出されたビーム応答ｙ［ｐ，ｋ］が時系列的に入力される。積分処理部８は、入力されたビーム応答ｙ［ｐ，ｋ］の絶対値を２乗した後に時間積分を行う。
【００７８】
［整相器の動作］
次に、上記構成を有する整相器１０の動作について、図４を参照して説明する。まず、アレイを構成する複数のセンサから整相器１０に対してセンサ出力ｘ［ｋ］が入力されると、入力されたセンサ出力ｘ［ｋ］は、共分散行列算出部１および整相処理部７に供給される。
【００７９】
共分散行列算出部１は、供給されたセンサ出力ｘ［ｋ］に基づき、共分散行列の推定値Ｒ＾［ｋ］を算出する。算出された共分散行列の推定値Ｒ＾［ｋ］は、固有値分解処理部２に供給される。
【００８０】
固有値分解処理部２は、供給された共分散行列の推定値Ｒ＾［ｋ］に基づき、固有値分解処理を行う。そして、信号部分空間および雑音部分空間の分離を行うとともに、信号部分空間の固有値λｍ［ｋ］および固有ベクトルｑｍ［ｋ］、ならびに雑音部分空間の固有値λｍ［ｋ］および固有ベクトルｑｍ［ｋ］を取得する。
【００８１】
取得した雑音部分空間に属する固有ベクトルｑｍ［ｋ］は、ＭＵＳＩＣ処理部３に供給される。また、取得した信号部分空間に属する固有ベクトルｑｍ［ｋ］は、方位インデックス決定部４に供給される。さらに、取得した信号部分空間に属する固有値λｍ［ｋ］および固有ベクトルｑｍ［ｋ］、ならびに雑音部分空間に属する固有値λｍ［ｋ］は、適応重み算出部６に供給される。
【００８２】
ＭＵＳＩＣ処理部３は、供給された雑音部分空間に属する固有ベクトルｑｍ［ｋ］に基づきＭＵＳＩＣ処理を行い、ＭＵＳＩＣスペクトルＭ［ｐ，ｋ］を算出する。算出されたＭＵＳＩＣスペクトルＭ［ｐ，ｋ］は、方位インデックス決定部４に供給される。
【００８３】
方位インデックス決定部４は、供給されたＭＵＳＩＣスペクトルＭ［ｐ，ｋ］および信号部分空間に属する固有ベクトルｑｍ［ｋ］に基づき、信号部分空間に属する固有ベクトルｑｍ［ｋ］に対応する方位インデックスｐｍを決定する。また、方位インデックス決定部４は、方位インデックスｐｍを中心とした前後複数の方位インデックスｐｍ−（Ｕ−１）／２，・・・，ｐｍ−（Ｕ−１）／２を、信号部分空間に属する固有ベクトルｑｍ［ｋ］に同時に対応付ける。決定された方位インデックスｐｍは、固有ベクトル除外重み算出部５および適応重み算出部６に供給される。
【００８４】
固有ベクトル除外重み算出部５は、供給された方位インデックスｐｍに基づき、固有ベクトル除外重みδ［ｍ，ｐ］を導出する。導出された固有ベクトル除外重みδ［ｍ，ｐ］は、適応重み算出部６に供給される。
【００８５】
適応重み算出部６は、供給された信号部分空間の固有値λｍ［ｋ］および固有ベクトルｑｍ［ｋ］、ならびに雑音部分空間の固有値λｍ［ｋ］と、方位インデックスｐｍと、固有ベクトル除外重みδ［ｍ，ｐ］とに基づき、適応重みｗ［ｐ，ｋ］を算出する。算出された適応重みｗ［ｐ，ｋ］は、整相処理部７に供給される。
【００８６】
整相処理部７は、供給されたセンサ出力ｘ［ｋ］と、適応重みｗ［ｐ，ｋ］とに基づき、ビーム応答ｙ［ｐ，ｋ］を算出する。算出されたビーム応答ｙ［ｐ，ｋ］は、積分処理部８に供給される。積分処理部８は、供給されたビーム応答ｙ［ｐ，ｋ］の絶対値を２乗した後に時間積分を行い、ビーム応答ｙ［ｐ，ｋ］の積分値を出力する。
【００８７】
［整相処理結果］
図６は、本実施の形態１に係る整相器１０による整相処理結果について説明するための概略図である。図６は、２つの信号＃１および信号＃２が到来する場合の整相処理の結果を示す。図６（ａ）は、従来型ＥＢＡＥ方式を用いた場合の整相処理の結果を示す。図６（ｂ）は、改善型ＥＢＡＥ方式を用いた場合の整相処理の結果を示す。図６（ａ）および図６（ｂ）において、縦軸はビーム応答の出力レベルを示し、横軸は整相方位を示す。
【００８８】
また、図６（ｃ）は、ＤＩＦＡＲに対して到来する２つの信号＃１および信号＃２の方位およびレベルを示す。図６（ｃ）に示すように、この例では、信号＃１は、真方位１６０［ｄｅｇ］からレベルが０［ｄＢ］で入射する。また、信号＃２は、真方位２５０［ｄｅｇ］からレベルが０［ｄＢ］で入射する。
【００８９】
図６（ａ）に示すように、従来型ＥＢＡＥ方式を用いた場合には、信号＃１および信号＃２の各信号の到来方位に対して複数のピークが出現する。そのため、この結果からは、各信号の到来方位を正確に判断することができない。
【００９０】
一方、図６（ｂ）に示すように、改善型ＥＢＡＥ方式を用いた場合には、信号＃１および信号＃２の各信号の到来方位に対して１つのピークが出現する。また、このピークのレベルは、各信号のレベルと同様の０［ｄＢ］となっている。これにより、信号＃１および信号＃２の到来方位およびレベルが正確に検出できていることがわかる。
【００９１】
このように、本実施の形態１に係る整相器１０では、改善型ＥＢＡＥ方式を用いることにより、従来のＥＢＡＥ方式を用いた場合に、１つの信号に対して複数のピークが発生するという問題を解決することができる。
【００９２】
以上のように、本実施の形態１に係る整相器１０は、アレイを構成する複数のセンサに到来する信号を示すセンサ出力信号に基づき、アレイの周囲に存在する信号源に関する情報を取得するものである。この整相器１０は、センサ出力信号に基づき、共分散行列を算出する共分散行列算出部１と、共分散行列に基づき、空間を信号部分空間および雑音部分空間に分離するとともに、信号部分空間および雑音部分空間に属する固有ベクトルを取得する固有値分解処理部２と、雑音部分空間に属する固有ベクトルに基づき、センサ出力信号の到来方位を推定するＭＵＳＩＣ処理部３と、推定された到来方位と、信号部分空間に属する固有ベクトルとに基づき、固有ベクトルに対応付ける到来方位を示す方位インデックスを決定する方位インデックス決定部４と、決定された方位インデックスに基づき、方位インデックスに対応付けられた固有ベクトルの成分を除外する除外重みを算出する固有ベクトル除外重み算出部５と、信号部分空間に属する固有ベクトルと、算出された除外重みとに基づき、各方位に対応する適応重みを算出する適応重み算出部６とを備えている。
【００９３】
方位インデックス決定部４は、信号部分空間に属する固有ベクトルに対応付ける方位インデックスを決定する際に、決定された方位インデックスの周囲の複数の方位を示す方位インデックスを、信号部分空間に属する固有ベクトルに対応付ける。また、固有ベクトル除外重み算出部５は、決定された方位インデックスに対応付けられた複数の方位について、除外重みを設定する。
【００９４】
このように、本実施の形態１では、信号部分空間に属する固有ベクトルに対応付ける方位インデックスを決定する際に、決定された方位インデックスの周囲の複数の方位を示す方位インデックスを、信号部分空間に属する固有ベクトルに対応付け、決定された方位インデックスに対応付けられた複数の方位について、除外重みを設定する。これにより、１つの信号に対する複数のピークの発生を抑制し、信号の到来方位を容易に特定することができる。
【００９５】
実施の形態２．
次に、本実施の形態２に係る整相器について説明する。本実施の形態２では、ＭＵＳＩＣ処理部３で算出されたＭＵＳＩＣスペクトルＭ［ｐ，ｋ］のばらつきの大きさに応じて固有ベクトル除外重みδ［ｍ，ｐ］を変化させる点で、上述した実施の形態１と相違する。
【００９６】
上述した実施の形態１では、固定的な固有ベクトル除外重みδ［ｍ，ｐ］を用いている。そのため、例えば、ＭＵＳＩＣスペクトルのピークのばらつきが小さい場合であっても、ビーム応答の幅が必要以上に広がってしまう可能性がある。そこで、本実施の形態２では、このようなビーム応答幅の広がりを抑制するようにした。
【００９７】
［整相器の構成］
図７は、本実施の形態２に係る改善型ＥＢＡＥ方式を用いた整相器１０の構成の一例を示すブロック図である。なお、以下の説明において、上述した実施の形態１と共通する部分には、同一の符号を付し、詳細な説明を省略する。
【００９８】
図７に示すように、整相器１０は、共分散行列算出部１、固有値分解処理部２、ＭＵＳＩＣ処理部３、方位インデックス決定部４、推定平均分散算出部２１、固有ベクトル除外重み算出部２５、適応重み算出部６、整相処理部７、および積分処理部８で構成されている。なお、共分散行列算出部１、固有値分解処理部２、ＭＵＳＩＣ処理部３、方位インデックス決定部４、適応重み算出部６、整相処理部７、および積分処理部８については、実施の形態１と同様であるため、その構成および動作についての説明を省略する。
【００９９】
推定平均分散算出部２１は、時刻ｊにおける方位インデックス決定部４の出力である信号部分空間に属する固有ベクトルｑｍ，ｊに対応付ける方位インデックスの中心ｐｍ，ｊが入力される。なお、この方位インデックスｐｍ，ｊは、ＭＵＳＩＣスペクトルによって推定された、信号到来方位に最も近い方位インデックスであるため、以下では、「信号到来方位」と適宜称して説明する。
【０１００】
推定平均分散算出部２１は、時刻ｊ−Ｈ＋１から時刻ｊまでの各時刻におけるＨ個の信号到来方位から、時刻ｊにおける推定平均μ＾ｍ，ｊおよび推定分散σ＾２ｍ，ｊを算出する。推定平均μ＾ｍ，ｊは、式（１６）に基づき算出することができる。また、推定分散σ＾２ｍ，ｊは、式（１７）に基づき算出することができる。
【０１０１】
【数１６】
2018189602.tif 000018
【０１０２】
【数１７】
2018189602.tif 000019
【０１０３】
固有ベクトル除外重み算出部２５は、方位インデックス決定部４で決定された信号到来方位ｐｍ，ｊと、推定平均分散算出部２１で算出された信号到来方位ｐｍ，ｊの推定平均μ＾ｍ，ｊおよび推定分散σ＾２ｍ，ｊとが入力される。固有ベクトル除外重み算出部２５は、入力された信号到来方位ｐｍ，ｊと、信号到来方位ｐｍ，ｊの推定平均μ＾ｍ，ｊおよび推定分散σ＾２ｍ，ｊとに基づき、固有ベクトル除外重みδｊ［ｍ，ｐ］を算出する。方位インデックスｐが固有ベクトルｑｍに対応付けられている場合、固有ベクトル除外重みδｊ［ｍ，ｐ］は、式（１８）に基づき算出することができる。
【０１０４】
【数１８】
2018189602.tif 000020
【０１０５】
式（１８）において、２次関数を決定する定数ａｍ，ｊは、信号到来方位ｐｍ，ｊの推定分散σ＾２ｍ，ｊの関数「ａｍ，ｊ＝ｆ（σ＾２ｍ，ｊ）」とする。ここで、図８は、本実施の形態２における信号到来方位ｐｍ，ｊの推定分散σ＾２ｍ，ｊと定数ａｍ，ｊとの関係の一例を示すグラフである。関数「ａｍ，ｊ＝ｆ（σ＾２ｍ，ｊ）」によって表される推定分散σ＾２ｍ，ｊと定数ａｍ，ｊとの関係は、この図８に示すように設定される。
【０１０６】
図９は、本実施の形態２における固有ベクトル除外重みの形状の一例を示す概略図である。図９（ａ）および図９（ｂ）は、それぞれ、信号到来方位ｐｍ，ｊの推定分散σ＾２ｍ，ｊが大きい場合および小さい場合の固有ベクトル除外重みδｊ［ｍ，ｐ］の形状の一例を示す。図９に示すように、固有ベクトル除外重みδｊ［ｍ，ｐ］は、信号到来方位ｐｍ，ｊの推定分散σ＾２ｍ，ｊの大きさによって形状が変化する。
【０１０７】
［整相処理結果］
図１０は、本実施の形態２に係る整相器１０による整相処理結果について説明するための概略図である。図１０は、図６と同様に、２つの信号＃１および信号＃２が到来する場合の整相処理の結果を示す。到来する信号＃１および信号＃２の状態は、実施の形態１と同様である（図６（ｃ）参照）。
【０１０８】
図１０（ａ）は、実施の形態１における改善型ＥＢＡＥ方式を用いた場合の整相処理の結果を示す。図１０（ｂ）は、本実施の形態２における改善型ＥＢＡＥ方式を用いた場合の整相処理の結果を示す。図１０（ａ）および図１０（ｂ）において、縦軸はビーム応答の出力レベルを示し、横軸は整相方位を示す。
【０１０９】
図１０に示すように、実施の形態１における改善型ＥＢＡＥ方式を用いた場合には、ビーム応答の幅がなだらかである。これに対して、本実施の形態２における改善型ＥＢＡＥ方式を用いた場合には、ビーム応答の幅が細くなっていることがわかる。
【０１１０】
このように、本実施の形態２に係る整相器１０では、ＭＵＳＩＣスペクトルのばらつきの度合いに応じて固有ベクトル除外重みδｊ［ｍ，ｐ］を変化させることにより、より適切に適応整相処理を行うことができる。
【０１１１】
以上のように、本実施の形態２に係る整相器１０は、上述した実施の形態１に係る構成に加えて、決定された方位インデックスの時間的な平均および分散を推定する推定平均分散算出部２１をさらに備え、固有ベクトル除外重み算出部２５は、複数の方位について設定された除外重みを算出するための方位インデックスの関数を、推定された方位インデックスの平均および分散に基づき変化させる。これにより、ＭＵＳＩＣスペクトルのばらつきが小さい場合に、ビーム応答の幅が必要以上に広がることを抑制し、信号の到来方位をより正確に特定することができる。
【０１１２】
以上、実施の形態１および２について説明したが、本発明は、上述した実施の形態１および２に限定されるものではなく、本発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。例えば、上述した例では、３つのセンサで構成されるＤＩＦＡＲを対象として、ＥＢＡＥ方式の適応整相処理を適用したが、これはこの例に限られない。例えば、ＤＩＦＡＲ以外の少数のセンサで構成されるアレイに対してＥＢＡＥ方式の適応整相処理を適用した場合でも、同様に１つの信号に対して複数のピークが発生することが考えられる。その場合には、本実施の形態１および２と同様に、１つの信号に対する複数のピークの発生を抑制することができる。
【符号の説明】
【０１１３】
１、１０１  共分散行列算出部、２、１０２  固有値分解処理部、３、１０３  ＭＵＳＩＣ処理部、４、１０４  方位インデックス決定部、５、１０５  固有ベクトル除外重み算出部、６、１０６  適応重み算出部、７、１０７  整相処理部、８、１０８  積分処理部、１０、１００  整相器、２１  推定平均分散算出部、２５  固有ベクトル除外重み算出部。
【図１】
2018189602.tif 000021
【図２】
2018189602.tif 000022
【図３】
2018189602.tif 000023
【図４】
2018189602.tif 000024
【図５】
2018189602.tif 000025
【図６】
2018189602.tif 000026
【図７】
2018189602.tif 000027
【図８】
2018189602.tif 000028
【図９】
2018189602.tif 000029
【図１０】
2018189602.tif 000030