特許 6767807 震源位置推定方法及び震源位置推定システム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6767807

(24)【登録日】2020年9月24日

(45)【発行日】2020年10月14日

(54)【発明の名称】震源位置推定方法及び震源位置推定システム

(51)【国際特許分類】

   G01S 5/22 20060101AFI20201005BHJP

   G01V 1/00 20060101ALI20201005BHJP

【ＦＩ】

   G01S5/22

   G01V1/00 C

【請求項の数】4

【全頁数】16

(21)【出願番号】特願2016-162296(P2016-162296)

(22)【出願日】2016年8月22日

(65)【公開番号】特開2018-31606(P2018-31606A)

(43)【公開日】2018年3月1日

【審査請求日】2019年3月19日

(73)【特許権者】

【識別番号】000000099

【氏名又は名称】株式会社ＩＨＩ

(73)【特許権者】

【識別番号】899000068

【氏名又は名称】学校法人早稲田大学

(74)【代理人】

【識別番号】100118267

【弁理士】

【氏名又は名称】越前 昌弘

(72)【発明者】

【氏名】小笹 弘晃

(72)【発明者】

【氏名】佐藤 文男

(72)【発明者】

【氏名】田中 浩一郎

(72)【発明者】

【氏名】及川 靖広

(72)【発明者】

【氏名】矢田部 浩平

(72)【発明者】

【氏名】草野 翼

【審査官】 安井 英己

(56)【参考文献】

【文献】 特開昭５５−１２１１７２（ＪＰ，Ａ）

【文献】 特開昭６２−２８８５８６（ＪＰ，Ａ）

【文献】 特開平０４−２３８２８４（ＪＰ，Ａ）

【文献】 特開平０６−３４７５３０（ＪＰ，Ａ）

【文献】 米国特許出願公開第２０１２／００２０１８４（ＵＳ，Ａ１）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０１Ｓ ５／１８− ５／３０，

Ｇ０１Ｓ ７／５２− ７／６４，

Ｇ０１Ｓ １５／００−１５／９６，

Ｇ０１Ｖ １／００− １／４０，

ＣＳＤＢ

(57)【特許請求の範囲】

【請求項1】

水中に沈められた震源から所定の音波を発震する発震ステップと、
水中の音波を略直線状に配列された複数の受振器で観測する観測ステップと、
観測した信号から前記震源の直接波を抽出する抽出ステップと、
抽出した直接波の信号から前記震源の位置を推定する推定ステップと、を備え、
前記推定ステップは、前記震源の位置をパラメータとして、あるパラメータを選択し、そのときの到達時間における各受振器の観測信号からラドン変換の値の最大値又は勾配を計算することによって前記震源の位置を推定する、
ことを特徴とする震源位置推定方法。

【請求項2】

前記発震ステップは、前記震源から既知の周波数帯域の音波を発震する、ことを特徴とする請求項１に記載の震源位置推定方法。

【請求項3】

水中に沈められ所定の音波を発震可能な震源と、
略直線状に配列され水中の音波を観測可能な複数の受振器と、
前記複数の受振器の観測信号から前記震源の直接波を抽出して前記震源の位置を推定する制御装置と、を備え、
前記制御装置は、前記震源の位置をパラメータとして、あるパラメータを選択し、そのときの到達時間における各受振器の観測信号からラドン変換の値の最大値又は勾配を計算することによって前記震源の位置を推定するように構成されている、
ことを特徴とする震源位置推定システム。

【請求項4】

前記震源は、音波の発震波形を任意に調整可能な音波発生装置である、ことを特徴とする請求項３に記載の震源位置推定システム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、震源位置推定方法及び震源位置推定システムに関し、特に、水中に沈めた震源の位置を特定するための震源位置推定方法及び震源位置推定システムに関する。

【背景技術】

【0002】

大陸棚や深海底には豊富な資源（例えば、石油、天然ガス、メタンハイドレート、多金属団塊、マンガン・クラスト、海底熱水鉱床等）が存在しており、近年の資源価格の高騰により海洋資源開発の必要性が高まっている。また、陸上における天然資源は、一定の地域に偏在しており、国内産出量が少ない資源については、外国からの輸入に頼らざるを得ず、地政学的リスクが少なくない。そして、四方を海に囲まれた我が国においては、安定した資源供給のためにも、海洋地域が資源開発の新たなフロンティアとして注目されている。

【0003】

かかる海底資源探査方法としては、例えば、特許文献１に記載されたように、所定の帯域の低周波を発震可能な音源を用い、その反射波をストリーマケーブル（曳航型受振器）で計測する方法が既に提案されている。

【0004】

ところで、震源から発震される音波の反射波を計測する探査方法では、一般に、震源位置と受波位置に基づく幾何学的な考察により地中の物理的特性を予測している。したがって、より高精度な探査を実現するためには、震源位置を正確に把握することが不可欠である。

【0005】

水中測位方法としては、例えば、特許文献２に記載されたように、船側に配置された音響測位装置と水中に沈められる構造物に配置されたトランスポンダとの間で所定の信号を送受信することによって水中の構造物の位置や方位を計測する方法が知られている。

【先行技術文献】

【特許文献】

【0006】

【特許文献1】特開２０１４−１３７３２０号公報

【特許文献2】特開２０１６−６５７５０号公報

【発明の概要】

【発明が解決しようとする課題】

【0007】

近年、より深い海底の資源探査のニーズが高くなっているところ、特許文献２に記載された方法では、船側の音響測位装置と水中側のトランスポンダとの距離が離れてしまい測位精度が低下してしまうという問題があった。また、探査機器とは別の水中測位機器を必要とすることから、コスト及び故障リスクが増加してしまうという懸念もある。

【0008】

本発明は、上述した問題点に鑑み創案されたものであり、震源を水中に深く沈めた場合であっても、別の水中測位機器を用いることなく、震源の位置を正確に推定することができる、震源位置推定方法及び震源位置推定システムを提供することを目的とする。

【課題を解決するための手段】

【0009】

本発明によれば、水中に沈められた震源から所定の音波を発震する発震ステップと、水中の音波を略直線状に配列された複数の受振器で観測する観測ステップと、観測した信号から前記震源の直接波を抽出する抽出ステップと、抽出した直接波の信号から前記震源の位置を推定する推定ステップと、を備え、前記推定ステップは、前記震源の位置をパラメータとして、あるパラメータを選択し、そのときの到達時間における各受振器の観測信号からラドン変換の値の最大値又は勾配を計算することによって前記震源の位置を推定する、ことを特徴とする震源位置推定方法が提供される。

【0012】

また、前記発震ステップは、前記震源から既知の周波数帯域の音波を発震するようにしてもよい。

【0013】

また、本発明によれば、水中に沈められ所定の音波を発震可能な震源と、略直線状に配列され水中の音波を観測可能な複数の受振器と、前記複数の受振器の観測信号から前記震源の直接波を抽出して前記震源の位置を推定する制御装置と、を備え、前記制御装置は、前記震源の位置をパラメータとして、あるパラメータを選択し、そのときの到達時間における各受振器の観測信号からラドン変換の値の最大値又は勾配を計算することによって前記震源の位置を推定するように構成されている、ことを特徴とする震源位置推定システムが提供される。

【0014】

前記震源は、例えば、音波の発震波形を任意に調整可能な音波発生装置である。

【発明の効果】

【0016】

上述した震源位置推定方法及び震源位置推定システムによれば、水中に沈めた震源から所定の音波を発震し、複数の受振器で受信した信号から震源からの直接波を抽出し、抽出した直接波の信号から震源の位置を推定するようにしたことにより、いわゆる海底探査（海底資源及び地形の調査を含む。）で使用する機器以外の水中測位機器を用いることなく、震源の位置を正確に特定することができる。また、直接波から震源の位置を求めることができれば、深い位置まで震源を沈めた場合でも観測信号のみから海底の様子を探索することができる。

【図面の簡単な説明】

【0017】

【図1】本発明の一実施形態に係る震源位置推定システムを示す全体構成図である。

【図2】本発明の第一実施形態に係る震源位置推定方法を示すフロー図である。

【図3】推定ステップを示す概念図である。

【図4】本発明の第二実施形態に係る震源位置推定方法を示すフロー図である。

【図5】指令信号を示す図であり、（ａ）は波形、（ｂ）はスペクトログラム、を示している。

【図6】１回目の発震における観測信号を示す図であり、（ａ）は全ての受振器の信号、（ｂ）はｉ=１の受振器の信号、を示している。

【図7】１回目の発震における観測信号を示すパルス圧縮後の信号を示す図であり、（ａ）は全ての受振器の信号、（ｂ）はｉ=１の受振器の信号、を示している。

【図8】１回目の発震における観測信号を示す直接波抽出後の信号を示す図であり、（ａ）は全ての受振器の信号、（ｂ）はｉ=１の受振器の信号、を示している。

【図9】１回目の発震における観測信号を示す直接波抽出後かつパルス圧縮後の信号を示す図であり、（ａ）は全ての受振器の信号、（ｂ）はｉ=１の受振器の信号、を示している。

【図10】１回目の発震におけるラドン変換した結果を示す図であり、（ａ）はｘｙ平面、（ｂ）はｘｃ平面、（ｃ）はｙｃ平面、を示している。

【発明を実施するための形態】

【0018】

以下、本発明の実施形態及び実施例について図１〜図１０を用いて説明する。ここで、図１は、本発明の一実施形態に係る震源位置推定システムを示す全体構成図である。図２は、本発明の第一実施形態に係る震源位置推定方法を示すフロー図である。

【0019】

本発明の一実施形態に係る震源位置推定システム１は、例えば、図１に示したように、水中に沈められ所定の音波を発震可能な震源２と、略直線状に配列され水中の音波を観測可能な複数の受振器３ｉ（ｉは１以上の整数）と、複数の受振器３ｉの観測信号から震源２の直接波Ｗｄを抽出して震源２の位置を推定する制御装置４と、を備えている。

【0020】

震源２は、例えば、油圧サーボを用いた低周波発生装置である。低周波発生装置としては、例えば、特開平８−２８００８９号公報に記載されたように、振動板の振動を油圧で制御することにより水中に音波を発震させる水中音源が使用可能である。震源２として低周波発生装置を使用することにより、所望の周波数及び波形を有する音波を自在に発生させることができる。また、かかる震源２を使用することにより、既知の周波数帯域の音波を発震することができ、エアガンとは異なり再現性の高い観測を行うことができる。

【0021】

なお、震源２から水中に発震される音波は、例えば、１０〜１００Ｈｚ程度の低周波であることが好ましいが、上述した油圧サーボ式の低周波発生装置であれば、１００〜５００Ｈｚの周波数の音波まで発生させることもできる。また、震源２は、油圧サーボを用いた低周波発生装置に限定されるものではなく、発震する音波の周波数や発震時間等が調整可能である等、発震波形が明確な音波発生装置、すなわち、音波の発震波形を任意に調整可能な音波発生装置であればよい。また、震源２は、震動源として、圧電素子や電動機を用いたものであってもよい。

【0022】

震源２は、例えば、ケーブル２１を介して曳航船５に接続されている。震源２は、例えば、水深数十ｍ〜数百ｍあるいは数千ｍまでの深い水中に沈められる。かかる水深まで震源２を沈めると曳航船５との距離が離れすぎてしまい、ＵＳＢＬ（Ultra Short Base Line）のような水中測位機器では正確な位置を測定することが困難であり、ＧＰＳ（全地球測位システム）を用いて震源２の位置を直接的に測定することもできない。

【0023】

なお、震源２は、曳航船５に曳航されるものに限定されず、例えば、自走式水中航走体に接続したものであってもよいし、海底に固定したものであってもよい。

【0024】

複数の受振器３ｉは、例えば、複数のハイドロフォンをケーブルに連結したストリーマケーブル３（曳航型受振器）により構成される。ストリーマケーブル３は、例えば、一端がケーブルを介して曳航船５に接続され、他端がケーブルを介してブイ６に接続されている。曳航船５及びブイ６を用いることにより、ストリーマケーブル３を略直線状に展張することができる。また、ストリーマケーブル３は、図示したように、水中に沈められた状態で展張される。

【0025】

なお、「略直線状」とは、複数の受振器３ｉが文字通り直線状に配列した場合に限定されるものではなく、左右方向又は上下方向に蛇行している場合や湾曲している場合を含む趣旨であり、全体として一方向に配列した状態を保持することができていればよい。

【0026】

ストリーマケーブル３に配置される受振器３ｉの個数は任意であるが、少なくとも三個以上であることが好ましい。また、受振器３ｉは、震源２に近い側から離れる方向に向かってｉ＝１，２，３，…，Ｎと番号を振るものとし、以下の説明において、各受振器３ｉをチャンネルと呼ぶこともある。

【0027】

制御装置４は、震源２から発震される音波の周波数、波形、発震時刻等の制御を行う他、震源２の位置を推定する処理を行う。なお、制御装置４は、震源２の制御機能と震源２の位置推定機能とを分離して別の制御装置で処理するようにしてもよい。また、震源２の位置推定機能を有する制御装置４は、曳航船５に配置される場合に限定されず、ブイ６に配置されていてもよいし、専用の作業船に配置されていてもよいし、地上基地に配置されていてもよい。また、ストリーマケーブル３から制御装置４へのデータ伝送は有線であってもよいし無線であってもよい。

【0028】

上述した震源位置推定システム１において、震源２から音波を発震すると、音波は水中で球形状に伝播する。各受振器３ｉは、震源２からダイレクトに伝播した音波（直接波Ｗｄ）、海底や地層等に反射して伝播した音波（反射波Ｗｒ）及びノイズ（波の音、船舶のプロペラ音、海洋生物の発する音等）を受信する。各受振器３ｉが受信した信号（観測信号）は制御装置４に伝送される。

【0029】

一般に、海底探査では、反射波Ｗｒの観測信号を使用して海底の物理的特性を予測している。したがって、通常の海底探査では、震源２の直接波Ｗｄはノイズとして扱われる。本実施形態では、このノイズとして捨てられていた直接波Ｗｄの観測信号を用いて震源２の位置推定を行っている。

【0030】

本発明の第一実施形態に係る震源位置推定方法は、図２に示したように、水中に沈められた震源２から所定の音波を発震する発震ステップＳｔｅｐ１と、水中の音波を略直線状に配列された複数の受振器３ｉで観測する観測ステップＳｔｅｐ２と、観測信号から震源２の直接波Ｗｄを抽出する抽出ステップＳｔｅｐ３と、抽出した直接波Ｗｄの観測信号をパルス圧縮するパルス圧縮ステップＳｔｅｐ４と、直接波Ｗｄの観測信号から震源２の位置を推定する推定ステップＳｔｅｐ５と、を備えている。

【0031】

発震ステップＳｔｅｐ１は、図１に示した震源２から既知の周波数帯域の音波を発震するステップである。上述したように、震源２として低周波発生装置を使用することにより、帯域制限された低周波であって予め設定した波形を有する音波（指令信号）を任意に発震させることができる。

【0032】

観測ステップＳｔｅｐ２は、図１に示したストリーマケーブルに内蔵された複数の受振器３ｉにより水中の音波を観測する。各受振器３ｉが受信した音波（観測信号）には、上述したように、震源２の直接波Ｗｄ、反射波Ｗｒ及びノイズが含まれていることから、観測信号から直接波Ｗｄの観測信号のみを抽出する必要がある。

【0033】

抽出ステップＳｔｅｐ３は、受振器３ｉの観測信号から反射波Ｗｒ及びノイズを除去するステップである。例えば、抽出ステップＳｔｅｐ３は、ＦＦＴ（高速フーリエ変換）を行うことにより周波数領域で指令信号に含まれていない帯域を取り除き、窓掛け（重み付けを行う窓関数を掛ける処理）及び逆ＦＦＴを行った後、ＳＴＦＴ（短時間フーリエ変換）を行うことにより時間周波数領域で直接波Ｗｄ以外の成分を取り除き、窓掛け及び逆ＳＴＦＴを行う。

【0034】

本実施形態では、震源２に低周波発生装置を使用していることから、周波数の時間変動が既知であり、時間変動に沿った抽出を行うことができる。また、ＦＦＴによる抽出とＳＴＦＴによる抽出とに二段階に分けることにより、ＳＴＦＴで時間を区切ることによって周波数分解能が低下して低帯域の周波数を取り除くことが難しい場合であっても、ＦＦＴによる抽出によって低帯域の周波数を効果的に取り除くことができる。なお、直接波Ｗｄの抽出方法は上述した方法に限定されるものではない。

【0035】

パルス圧縮ステップＳｔｅｐ４は、直接波Ｗｄの観測信号と震源２の指令信号との相互相関をとることによりピーク値を際立たせるステップである。パルス圧縮することにより、小さな振幅を保ったまま幅の長いパルスを変調させた指令信号を送信した後、観測信号と指令信号との相互相関をとることによって、振幅が大きく幅の短いパルスを送信した場合と同じ効果を得ることができ、高い距離分解能を保持することができる。

【0036】

推定ステップＳｔｅｐ５は、例えば、図示したように、震源２の位置をパラメータａとして、あるパラメータａを選択し、そのときの到達時間ｔにおける各受振器３ｉの観測信号からラドン変換Ｒ（ａ）の値の最大値を計算することによって震源２の位置を推定するステップである。

【0037】

具体的には、図２に示したように、推定ステップＳｔｅｐ５は、初期パラメータを設定する第一ステップＳｔｅｐ５１と、ラドン変換Ｒ（ａ）の値の最大値を計算する第二ステップＳｔｅｐ５２と、所定の停止条件を満たしている否か判定する第三ステップＳｔｅｐ５３と、を備え、所定の停止条件に達するまで第二ステップＳｔｅｐ５２を繰り返す。

【0038】

ここで、図３は、推定ステップを示す概念図である。図３において、横軸はチャンネル（受振器３ｉ）を示し、縦軸は時間（下向きを正とする。）を示している。また、図３は、チャンネルごとに観測信号の時間変動を図示したものである。また、ｔ_ｉ（ｘ，ｙ，ｐ）は、震源２の位置（Ｘ座標及びＹ座標）及びスローネス（音速の逆数）ｐを特定した場合における指令信号の各チャンネルでの到達時間を示している。

【0039】

震源２の座標は、ｉ＝１の受振器３ｉ（１番目のチャンネル）を原点とし、曳航船５に向かう水平方向を正とするＸ座標とし、鉛直方向下向きを正とするＹ座標によって表現される。なお、音速ｃは、水深、温度及び塩分濃度によって変化する変数である。

【0040】

図示したように、指令信号の到達時間（観測時間）は、チャンネルの番号が大きくなるに連れて遅くなり、ｔ_ｉは点線で示した双曲線状に表示することができる。また、双曲線状における各チャンネルにおける観測信号の振幅は、チャンネルの番号が大きくなるに連れて小さくなる。したがって、震源２の位置及び音速を推定して特定し、そのときの各チャンネルにおける指令信号の到達時間における振幅を足し合わせることができれば、その加算値から尤もらしい震源２の位置を推定することが可能となる。

【0041】

例えば、震源２の位置及びスローネスｐを（ｘ_１，ｙ_１，ｐ）、（ｘ₂，ｙ₂，ｐ）、（ｘ₃，ｙ₃，ｐ）と三種類の推定をした場合において、到達時間ｔ_ｉ（ｘ_１，ｙ_１，ｐ）を通る線上の各チャンネルにおける振幅を足し合わせた場合、到達時間ｔ_ｉ（ｘ₂，ｙ₂，ｐ）を通る線上の各チャンネルにおける振幅を足し合わせた場合、到達時間ｔ_ｉ（ｘ₃，ｙ₃，ｐ）を通る線上の各チャンネルにおける振幅を足し合わせた場合を考察する。

【0042】

このとき、到達時間ｔ_ｉ（ｘ_１，ｙ_１，ｐ）の加算値が最も大きく、到達時間ｔ_ｉ（ｘ₂，ｙ₂，ｐ）の加算値、到達時間ｔ_ｉ（ｘ₃，ｙ₃，ｐ）の加算値と順に小さくなることは図から容易に理解することができる。すなわち、到達時間ｔ_ｉ（ｘ_１，ｙ_１，ｐ）の推定値が尤もらしい震源２の位置であると推定することができる。

【0043】

ここで、ラドン変換とは、パラメータで変化する線に沿う線積分をとる変換である。各チャンネルの信号を足し合わせることで、個々のノイズの影響を小さくし、尤もらしいパラメータを推定することができる。震源２と受振器３ｉとの位置関係から、指令信号の到達時間ｔ_ｉ（ｘ，ｙ，ｐ）は、ｐ｛（ｘ−ｈ_ｉ）^２＋ｙ^２｝^1/2と表現することができ、受振器３ｉの位置に対して双曲線状に変化する。なお、ｈ_ｉはｉ番目のチャンネルの原点からの水平距離を示している。

【0044】

いま、パラメータａ＝（ｘ，ｙ，ｐ）と定義すれば、ｔ_ｉ（ａ）に沿って線積分するラドン変換は、直接波Ｗｄを抽出した後の観測信号と指令信号との相互相関をとった信号のｉ番目のチャンネルの時刻ｔにおける振幅をＤ_ｉ（ｔ）とすると、式（１）のように表される。

【0045】

【数1】

【0046】

式（１）は、振幅の最大値を多く通る双曲線において最大となる。つまり、式（１）が最大となるパラメータａが尤もらしいパラメータとなる。なお、ここでの観測信号は、時間に対して離散的であるころから観測信号を一次補間して、時刻ｔ_ｉにおける振幅を求めるようにしてもよい。

【0047】

上述したように、ラドン変換で得られた領域の最大となる点から震源２の位置を推定することが可能であるが、探索範囲の全てを計算する必要がある。したがって、第三ステップＳｔｅｐ５３における停止条件は、例えば、探索範囲の全てについてラドン変換したか否かに設定される。

【0048】

一方で、この方法では、探索範囲の全てについてラドン変換しなければならないことから処理に時間を要することとなる。そこで、ラドン変換Ｒ（ａ）を目的関数とし、目的関数を最大化する解を求める非線形計画問題として、式（２）を勾配法で解くことにより、パラメータ推定を高速化するようにしてもよい。

【0049】

【数2】

【0050】

ここで、図４は、本発明の第二実施形態に係る震源位置推定方法を示すフロー図である。図４に示した震源位置推定方法の推定ステップＳｔｅｐ５は、震源２の位置をパラメータａとして、あるパラメータａを選択し、そのときの到達時間における各受振器３ｉの観測信号からラドン変換Ｒ（ａ）の値の勾配を計算することによって震源２の位置を推定するようにしたものである。

【0051】

具体的には、図４に示したように、推定ステップＳｔｅｐ５は、初期パラメータを設定する第一ステップＳｔｅｐ５１と、ラドン変換Ｒ（ａ）の値の勾配を計算する第四ステップＳｔｅｐ５４と、パラメータａの更新を行う第五ステップＳｔｅｐ５５と、所定の停止条件を満たしている否か判定する第三ステップＳｔｅｐ５３と、を備え、所定の停止条件に達するまで第四ステップＳｔｅｐ５４及び第五ステップＳｔｅｐ５５を繰り返す。

【0052】

第四ステップＳｔｅｐ５４では、例えば、ネステロフの加速勾配法（Nesterov's Accelerated Gradient）を用いることができる。勾配法は、微分を用いて最適解を求める手法である。本実施形態では、観測信号を連続関数で近似することにより、ラドン変換した領域での微分を考えることができる。また、連続関数で近似することで、１次での近似よりも正確にラドン変換を求めることができ、より正確なパラメータを推定することができる。この連続関数の近似には、例えば、ｓｉｎｃ関数を用いることができる。

【0053】

一般に、標本化された信号は標本化定理より、標本化周波数ｆ_ｓの半分までの帯域のみを含む信号であれば、ｓｉｎｃ関数を畳み込むことによって元の連続信号を得ることができる。ここで、ｓｉｎｃ関数とは式（３）に示した式である。

【0054】

【数3】

【0055】

連続信号をｘ（ｔ）、標本化周波数ｆ_ｓで標本化された信号をｘ_ｓ（ｔ）とすると、ｘ_ｓ（ｔ）はディラックのデルタ関数を用いて式（４）のように表せる。ここで、Ｔ_ｓは標本化周期（＝１／ｆ_ｓ）、ｘ_ｊ＝ｘ（ｊＴ_ｓ）である。また、畳み込む関数ｇ（ｔ）は、ｇ（ｔ）＝ｆ_ｓｓｉｎｃ（ｆ_ｓｔ）と表せる。よって、連続信号ｘ（ｔ）は、式（４）から式（５）のように表せる。

【0056】

【数4】

【0057】

したがって、式（１）のラドン変換は、直接波Ｗｄの抽出後の観測信号と指令信号との相互相関をとった信号のｉ番目のチャンネルでのｋサンプル目の振幅をｄ（ｉ，ｊ）とすると式（６）のように表すことができる。

【0058】

【数5】

【0059】

ネステロフの加速勾配法は、最急降下法の勾配とパラメータの間に別のパラメータをいれることで更新による移動方向を反復全体で平均することにより、収束速度を速くするアルゴリズムである。最急降下法とは、勾配法の中で最も基本的なものであり、勾配のみからパラメータを更新するアルゴリズムである。

【0060】

Ｒ（ａ）を最大化するパラメータａを求めるとき、最急降下法は、ａ_ｋ＋１＝ａ_ｋ＋α∇Ｒ（ａ_ｋ）のようにパラメータを更新する。ここで、αは更新の重みを決めるパラメータであり、∇Ｒ（ａ）は、Ｒ（ａ）の勾配を表している。

【0061】

ネステロフの加速勾配法では、ａ_ｋ＋１＝ａ_ｋ＋ｖ_ｋ＋１のようにパラメータを更新する。ここで、ｖ_ｋ＋１＝μｖ_ｋ＋α∇Ｒ（ａ_ｋ＋μｖ_ｋ）であり、μ，αは更新の重みを決めるパラメータである。また、ｖは速度のパラメータであり、勾配∇Ｒ（ａ）によって更新される。速度パラメータｖに勾配∇Ｒ（ａ）が反復ごとに加算され平均化されることにより、移動方向を解の方向へ促すことができる。

【0062】

勾配Ｒ（ａ）は、式（７）のように表される。

【0063】

【数6】

【0064】

ここで、Ｒ_ｉｊ（ａ）＝ｄ（ｉ，ｊ）ｓｉｎｃ（ｆ_ｓｔ_ｉ（ａ）−ｊ＋１）であり、∂Ｒ_ｉｊ／∂ｔ_ｉは式（８）により表され、∇ｔｉ（ａ）は式（９）により表される。

【0065】

【数7】

【0066】

よって、式（７）〜式（９）により∇Ｒ（ａ）は式（10）のように表される。この式（10）を解くことにより、尤もらしい震源２の位置（Ｘ座標及びＹ座標）を推定することができる。

【0067】

【数8】

【0068】

かかる第二実施形態における停止条件は、例えば、微分値（勾配の大きさ）が閾値以下であるか否か、パラメータａの変化量が閾値以下であるか否か、ラドン変換Ｒ（ａ）の値の変化量が閾値以下であるか否か、ループ回数が閾値を超えたか否か等の中から任意に選択することができる。

【0069】

上述した震源位置推定方法及び震源位置推定システム１によれば、水中に沈めた震源２から所定の音波（指令信号）を発震し、複数の受振器３ｉで受信した観測信号から震源２の直接波Ｗｄを抽出し、該直接波Ｗｄの観測信号から震源２の位置を推定するようにしたことにより、いわゆる海底探査（海底資源及び地形の調査を含む。）で使用する機器以外の水中測位機器を用いることなく、震源２の位置を正確に推定することができる。

【0070】

また、ラドン変換Ｒ（ａ）を用いることにより、パラメータの推定処理をより正確かつ円滑に行うことができ、さらに、ラドン変換Ｒ（ａ）の勾配∇Ｒ（ａ）を用いることにより処理の高速化を図ることができる。なお、上述した推定ステップＳｔｅｐ５では、ラドン変換Ｒ（ａ）を用いた場合について説明したが、ラドン変換を用いたものに限定されず、最小二乗法等の他の手法を用いてもよいことはいうまでもない。

【0071】

なお、海底探査を行う際には、通常、ストリーマケーブル３（受振器３ｉ）の絶対位置は、曳航船５及び／又はブイ６に装備されたＧＰＳにより求められていることから、ストリーマケーブル３（受振器３ｉ）に対する震源２の相対位置を正確に推定することにより、震源２の絶対位置を容易に特定することができる。

【実施例1】

【0072】

ここで、上述した実施形態について、新潟県粟島沖で測定された実測データを用いて推定手法の検討を行った。ここでは、４８００回のデータのうち最初の２０回の発震によるデータに対して、直接波を抽出した後に各手法を適用した。図１に示した震源位置推定システム１を用いて発震を定期的に行った。実験条件を表１に示す。なお、以下の説明において、第三実施形態とは推定ステップで最小二乗法を用いたものである。

【0073】

【表1】

【0074】

また、発震する指令信号の波形を図５（ａ）に、指令信号のスペクトログラムを図５（ｂ）に示す。また、１回目の発震によって得られた観測信号を図６（ａ）に、その１チャンネル目の観測信号を図６（ｂ）に示す。また、観測信号と指令信号との相互相関をとった信号（すなわち、パルス圧縮後の信号）を図７（ａ）に、その１チャンネル目の信号を図７（ｂ）に示す。

【0075】

本実施形態に係る推定手法を適用する前処理として、観測信号から直接波の抽出を行った。まず、全体を離散フーリエ変換し、周波数が６〜１５０Ｈｚ以外の成分を取り除いた。さらに、ＳＴＦＴにより、スイープに沿って直接波のみを取り出した。窓関数にはブラックマンウィンドウを用い、窓長を０．２５６ｓ（２５６サンプル）、オーバーラップ長を０．２５０ｓ（２５０サンプル）とした。

【0076】

１回目の発震によって得られた観測信号から直接波を抽出した結果を図８（ａ）に、１チャンネル目の信号を図８（ｂ）に示す。また、１回目の発震の直接波抽出後の信号を指令信号と相互相関をとった信号を図９（ａ）に、１チャンネル目の信号を図９（ｂ）に示す。

【0077】

まず、第一実施形態に係る震源位置推定方法について説明する。かかる第一実施形態では、探索範囲での各点のラドン変換を計算し、最大となる点を求めることによりパラメータ推定を行った。

【0078】

計算する点は、５０≦ｘ≦１５０（ｍ）、２０≦ｙ≦１００（ｍ）の範囲を０．１ｍ刻みに定めた。スローネスｐ（ｓ／ｍ）は、音速ｃ（ｍ／ｓ）を１５０５≦ｃ≦１５１５（ｍ／ｓ）の範囲を０．１ｍ／ｓ刻みとし、ｐ＝１／ｃとしてスローネスｐ（ｓ／ｍ）を定め、合計１００１×８０１×１０１＝８０９８１９０１点とした。

【0079】

１回目の発震による直接波抽出を行った観測信号と指令信号との相互相関をとった信号をラドン変換した結果をｘｙ平面、ｘｃ平面、ｙｃ平面で表示したものをそれぞれ図１０（ａ）、図１０（ｂ）、図１０（ｃ）に示す。

【0080】

次に、第二実施形態に係る震源位置推定方法について説明する。かかる第二実施形態では、まず最急降下法によるパラメータ推定との収束速度の比較を行い、収束速度が速められていることを確かめ、次に２０回の発震データに対してパラメータ推定を行った。

【0081】

まず、三つのパラメータの桁数をそろえるために、スローネスｐの単位を（μｓ／ｍ）とし、振幅の単位を（×１０^１２）とした。更新の重みパラメータはα＝４、μ＝０．５ろした。初期値はａ_０＝（１００，５０，（１／１５００）×１０^６）とし、速度パラメータの初期値はｖ_０＝（０，０，０）とした。

【0082】

次に、ネステロフの加速勾配法と最急降下法とによるパラメータ推定の収束速度の比較を行ったところ、ネステロフの加速勾配法の方が少ない反復回数でパラメータ変化量及び最大値との差が小さくなり、ネステロフの加速勾配法を用いることで少ない反復回数で解を求めることができることが確かめられた。なお、反復回数は５００回とした。

【0083】

次に、第三実施形態に係る震源位置推定方法について説明する。かかる第三実施形態では、直接波の到達時間を求め、最小二乗法によりパラメータ推定を行った。

【0084】

まず、各チャンネルの直接波の到達時間を求める。ここでは、一つ前のチャンネルの到達時間から２０ｍｓ後の間の中で次のチャンネルの差分の最小値を取った。本実施形態では、距離が大きくなるにしたがって振幅が小さくなりノイズとの区別がつきにくく、到達時間を求めることが難しくなることから、１〜７０チャンネルまでのデータを用いて推定を行った。ニュートン法は、初期値ａ_０＝（１００，５０，１／１５００）、反復回数は１５回とした。

【0085】

上述した第一実施形態〜第三実施形態の推定手法を用いて、２０回の発震データに対する震源位置の推定を行った結果の平均値を表２に示す。

【0086】

【表2】

【0087】

ここで、震源の中心点とストリーマケーブルの中心点とのオフセット距離を他の機器で実測したところ、ｘ＝９６．１５５（ｍ）であった。また、震源の深さは、表１に記載したように、ｙ＝５０（ｍ）に設定した。これらのデータと表２のデータとを比較すれば、いずれの実施形態においても、震源の位置をほぼ正確に推定できていることがわかる。

【0088】

計算時間については、最小二乗法を用いた第三実施形態が最も速いが、各チャンネルにおける到達時間を正しく求めることができないような場合には誤差が大きくなってしまうことが予想される。また、ラドン変換の勾配を用いた第二実施形態では、ラドン変換の最大値を用いた第一実施形態と比較して、計算時間を約１／１００に縮めることができ、処理時間の高速化を実現することができる。

【0089】

本発明は上述した実施形態及び実施例に限定されず、本発明の趣旨を逸脱しない範囲で種々変更が可能であることは勿論である。

【符号の説明】

【0090】

１ 震源位置推定システム
２ 震源
３ ストリーマケーブル
３ｉ 受振器
４ 制御装置
５ 曳航船
６ ブイ
２１ ケーブル
Ｓｔｅｐ１ 発震ステップ
Ｓｔｅｐ２ 観測ステップ
Ｓｔｅｐ３ 抽出ステップ
Ｓｔｅｐ４ パルス圧縮ステップ
Ｓｔｅｐ５ 推定ステップ
Ｓｔｅｐ５１ 第一ステップ
Ｓｔｅｐ５２ 第二ステップ
Ｓｔｅｐ５３ 第三ステップ
Ｓｔｅｐ５４ 第四ステップ
Ｓｔｅｐ５５ 第五ステップ

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】