JP6637592B2 - 信号処理システム、信号処理方法、信号処理装置、及びプログラム - Google Patents

Description

本発明は、信号処理システム、信号処理方法、信号処理装置、及びプログラムに関し、特に地震探査等の振動に基づく信号を処理する技術に関する。

地盤調査、及び既存構造物内部の損傷調査等においては、調査対象となる媒質を破壊しないで媒質の特性を計測する非破壊計測が行われている。特許文献１には、複数の発震器を媒質に接するように設置し、これらの発震器から送出された信号を受信器で受信して解析することにより、媒質の特性を計測する方法が開示されている。この方法においては、夫々の発震器が異なるパターンの疑似ランダムコード信号に基づく震動波を送出し、受信器において、複数の発震器から送出された震動波が合波された合成震動波を受信する。発震器が送出した震動波と受信器が受信した合成震動波との相関性を用いることで、特定の震動波に関する情報が分離される。

このような手法は、複数の震源から同時に送出した信号を受信後に弁別する技術である。特に、連続する複数の送信信号を受信後に送信周期で加算処理する、いわゆるスタッキングと呼ばれる技術は、理論的に優れた相関性を得ることができる。この技術は長時間の測定時間を要する。

特開２００４−１６３３２２号公報

しかしながら、スタッキングと呼ばれる技術においては、スタック回数が少ない場合、発震器が送出した震動波と受信器が受信した合成震動波との相関性が悪化する。これを防止するため、あらかじめダミーの送信信号を送信する方法も提案されているが、例えば送信元源が頻繁に移動する場合などは弁別精度を上げられない問題があった。

本発明はこれらの点に鑑みてなされたものであり、複数の発震器から同時に送出した信号を受信後に弁別する際の弁別精度を高めることを目的とする。

本発明の第１の態様は、信号処理システムである。このシステムは、媒質に震動波を送出するＮ個（Ｎは１以上の整数）の発震器と、前記Ｎ個の発震器のそれぞれが送出して前記媒質を伝搬した震動波の合波を受信する受信器と、前記受信器が受信した合波の時系列データの末尾から所定の長さＤのデータを切り出して、前記時系列データの先頭から前記所定の長さＤまでの時系列データに加算した修正時系列データを生成する修正データ生成部と、前記Ｎ個の発震器がそれぞれ送出したＮ個の震動波の時系列データのそれぞれと、前記修正時系列データとの相関を示すデータを算出する相関算出部と、を備える。

前記Ｎ個の発震器のそれぞれは、前記所定の長さＤの範囲内において互いに無相関な疑似ランダム信号に基づく震動波を送出してもよい。

前記Ｎ個の発震器のそれぞれは互いに同期して前記震動波を送出してもよく、前記受信器は、前記Ｎ個の発震器のそれぞれによる前記震動波の送出と同期して前記時系列データの受信を開始してもよい。

前記Ｎ個の発震器のそれぞれは、前記所定の長さＤにＮ＋１を乗じた長さ以上の周期を有する疑似ランダム信号を、前記所定の長さのＭ倍（Ｍは０以上Ｎ−１以下の整数であって前記Ｎ個の発震器のそれぞれで異なる）の長さだけ循環的にシフトして得られた信号に基づく震動波を送出してもよい。

前記相関算出部は、前記修正時系列データを０から前記所定の長さＤにＮを乗じた長さに至るまで循環的にシフトしながら、順番に前記疑似ランダム信号との自己相関係数を算出してもよい。

前記相関算出部が算出したデータを解析して、前記Ｎ個の発震器がそれぞれ前記震動波を送出してから前記受信器に到達するまでに要した時間であるそれぞれの震動波の遅延時間を算出する遅延時間算出部をさらに備えてもよい。

前記Ｎ個の発震器はそれぞれ、前記震動波それぞれの遅延時間の最大値の時間が経過するまでの間に送出する震動波の長さが、前記所定の長さＤ以下となるように、前記震動波を送出してもよい。

本発明の第２の態様は、信号処理方法である。この方法は、媒質にＮ個（Ｎは１以上の整数）の震動波を送出するステップと、前記媒質を伝搬した前記Ｎ個の震動波の合波を受信して受信波の時系列データを生成するステップと、前記時系列データの末尾から所定の長さＤのデータを切り出して、前記時系列データの先頭から前記所定の長さＤまでの時系列データに加算した修正時系列データを生成するステップと、送出した前記Ｎ個の震動波の時系列データのそれぞれと、前記修正時系列データとの相関を示すデータを算出するステップと、をプロセッサが実行する。

本発明の第３の態様は、信号処理装置である。この装置は、媒質に震動波を送出するＮ個（Ｎは１以上の整数）の発震器のそれぞれが送出して前記媒質を伝搬した震動波の合波を受信して得られた受信波の時系列データを生成する時系列データ生成部と、前記時系列データの末尾から所定の長さＤのデータを切り出して、前記時系列データの先頭から前記所定の長さＤまでの時系列データに加算した修正時系列データを生成する修正データ生成部と、前記Ｎ個の発震器がそれぞれ送出したＮ個の震動波の時系列データのそれぞれと、前記修正時系列データとの相関を示すデータを算出する相関算出部と、を備える。

本発明の第４の態様は、プログラムである。このプログラムは、コンピュータに、媒質に震動波を送出するＮ個（Ｎは１以上の整数）の発震器のそれぞれが送出して前記媒質を伝搬した震動波の合波を受信して得られた受信波の時系列データを生成する機能と、前記時系列データの末尾から所定の長さＤのデータを切り出して、前記時系列データの先頭から前記所定の長さＤまでの時系列データに加算した修正時系列データを生成する機能と、前記Ｎ個の発震器がそれぞれ送出したＮ個の震動波の時系列データのそれぞれと、前記修正時系列データとの相関を示すデータを算出する機能と、を実現させる。

本発明によれば、複数の発震器から同時に送出した信号を受信後に弁別する際の弁別精度を高めることができる。

実施の形態に係る信号処理システムの概要を説明するための図である。 実施の形態に係る発震器及び受信器の配置例を示す図である。 実施の形態に係る発震器及び受信器の機能構成を模式的に示す図である。 実施の形態に係る疑似乱数生成部及び乱数シフト部が生成する疑似ランダム信号の一例を示す図である。 実施の形態に係る震動波発生部が生成する震動波を説明するための図である。 実施の形態に係る発震器が送出し、媒質を伝搬して受信器が受信した震動波の波形を模式的に示す図である。 実施の形態に係る修正データ生成部による修正時系列データの生成処理を説明するための図である。 第１発震器が送出した震動波の受信波形に関する修正時系列データを説明するための模式図である。 実施の形態に係る相関算出部が算出した自己相関係数の時間変化のグラフを例示する図である。 比較例に係る自己相関係数の算出方法を説明するための図である。 比較例に係る自己相関係数の算出方法を用いて算出した自己相関係数の時間変化のグラフを例示する図である。 地震探査における自己相関関数の実測データの例を示すグラフである。

図１は、実施の形態に係る信号処理システム１００の概要を説明するための図である。信号処理システム１００は、Ｎ個（Ｎは１以上の整数）の発震器１と、少なくとも一つの受信器２とを有する。図１に示す例ではＮの値は４であり、第１発震器１ａから第４発震器１ｄまでの４つの発震器１が図示されている。

複数の発震器１のそれぞれは、所定の長さＤの範囲内において互いに無相関な疑似ランダム信号を生成し、生成した信号に基づいて生成した震動波を、媒質中に送出する。なお、発震器１は互いに同期しており、震動波を同時に送出する。

「所定の長さＤ」とは、「所定の時間Ｔ」が経過するまでの間に各発震器１が送出する震動波の長さである。詳細は後述するが、「所定の時間Ｔ」とは、Ｎ個の発震器１のそれぞれが送出した震動波が、測定対象の媒質を伝搬して受信器２に至るまでに要する時間として見込まれる時間のうち最大の時間である。所定の時間Ｔは、いわば発震器１が震動波を送出してから受信器２に至るまでに要する時間として定義される「遅延時間」の最大値の見込みである。

Ｎ個の発震器１はそれぞれ、所定の時間Ｔが経過するまでの間に送出する震動波の長さが所定の長さＤ以下となるように、震動波を送出する。なお、所定の時間Ｔは遅延時間における最大値以上の長さであれば、最大値と同じ長さでなくてもよい。

受信器２は、媒質中を伝搬した震動波を受信する。受信器２は、Ｎ個の発震器１のそれぞれによる震動波の送出と同期して、発震器１が震動波を送出するのと同時に震動波の受信を開始する。ここで受信器２は、各発震器１が震動波の送出を終了してから少なくとも所定の長さＤの震動波を受信するまでの間は受信を継続する。

各発震器１は同じ媒質に対して同時に震動波を送出するため、受信器２は、複数の発震器１が発震することにより生成される複数の震動波が合成された合成震動波を受信する。受信器２が受信した合成振動は、信号処理装置３（例えば、コンピュータ）において合成震動波に対応する時系列データに変換され、解析される。

信号処理装置３は、受信器２が受信した合成震動波を時系列データに変換する。信号処理装置３はさらに、時系列データの末尾から所定の長さＤに対応する長さのデータを切り出す。信号処理装置３は、切り出したデータを、時系列データの先頭から所定の長さＤに対応する長さまでの時系列データに加算した修正時系列データを生成する。

信号処理装置３は、修正時系列データと、Ｎ個の受信器２がそれぞれ送出したＮ個の震動波それぞれとの相関を算出する。これにより、信号処理装置３は、受信器２が受信した合成震動波に基づく受信信号から、特定の発震器１が送出した震動波に基づく成分を分離することができ、Ｎ個の発震器１の夫々との間の媒質の特性を特定することができる。

より具体的には、信号処理装置３は、算出した相関値がピークとなるタイミングと、ピークにおける相関係数の大きさとを特定する。相関係数がピークになるタイミングは上述した遅延時間に対応し、この遅延時間から媒質における震動波の伝達速度が推定できる。また、相関係数のピークの大きさから、媒質を伝搬したことによる震動波の減衰量が推定できる。これらの結果から、信号処理装置３は、特定の発震器１と受信器２との間の媒質の特性を特定できる。

以上が、実施の形態に係る信号処理システム１００における信号処理装置３が実行する信号処理の概要である。信号処理装置３は、受信器２が受信した合成震動波の末尾から所定の長さＤの信号を切り出し、その信号を合成震動波の先頭に加算して生成される信号を相関算出の対象とすることが、実施の形態に係る信号処理システム１００の一つの特徴である。
以下、信号処理システム１００についてより詳細に説明する。

図２は、実施の形態に係る発震器１及び受信器２の配置例を示す図である。図２においては、４つの発震器１（第１発震器１ａ〜第４発震器１ｄ）、及び受信器２が示されている。夫々の発震器１が送出した震動波は、媒質Ｍの一例である地中を伝搬し、地中の物質の特性に応じて反射及び散乱をした後に受信器２に到達する。

衛星Ｓは、例えばＧＰＳ（Global Positioning System）等の航法衛星である。各発震器１及び受信器２は衛星Ｓが送信する電波を受信し、受信電波に含まれる時刻情報を利用することで互いに同期することができる。

受信器２はケーブル又は無線回線を介して信号処理装置３に接続されており、受信器２は受信した震動波に基づく信号を信号処理装置３に送る。信号処理装置３は、受信器２から受信した信号を解析することにより、媒質Ｍの特性を特定する。

図３は、実施の形態に係る発震器１及び受信器２の機能構成を模式的に示す図である。
発震器１は、疑似乱数生成部１１、乱数シフト部１２、震動波発生部１３、クロック生成部１４、及び震動波記憶部１５を備える。発震器１は、例えばバイブロサイスである。

クロック生成部１４は、衛星Ｓから受信した電波に基づいてクロックを生成する。クロック生成部１４は、生成したクロックを疑似乱数生成部１１に出力する。クロック生成部１４が生成するクロックは、衛星Ｓから受信した電波に含まれる時刻情報に同期しているので、信号処理システム１００が有するすべての発震器１は同期して動作することができる。

疑似乱数生成部１１は、クロック生成部１４から入力されたクロックに基づいて、所定の周期Ｐの長さを有する疑似ランダム信号を生成する。より具体的には、疑似乱数生成部１１は、所定の周期Ｐの長さを有するＭ系列信号より選択された、自己相関性疑似乱数パターンを生成する。このような、自己相関性疑似乱数パターンを適切な変調方式で表現した疑似ランダム信号は、ラグタイム０において１であり、ラグタイムが０以外の場合に−１／Ｐとなる。

ここで所定の周期Ｐは、上述した所定の長さＤに、信号処理システム１００が備える発震器１の数Ｎ＋１を乗じた長さ以上の値である。すなわち、
Ｐ≧（Ｎ＋１）Ｄ （１）
を満たす長さである。

乱数シフト部１２は、疑似乱数生成部１１が生成した疑似ランダム信号を、所定の長さＤのＭ倍の長さだけ循環的にシフトした信号を生成する。ここでＭは０以上Ｎ−１以下の整数であって、Ｎ個の発震器１のそれぞれで異なる。実施の形態に係る信号処理システム１００の場合、第１発震器１ａ〜第４発震器１ｄに、それぞれＭ＝０〜３が割り当てられる。

図４（ａ）−（ｂ）は、実施の形態に係る疑似乱数生成部１１及び乱数シフト部１２が生成する疑似ランダム信号の一例を示す図である。より具体的には、図４（ａ）は、各発震器１の疑似乱数生成部１１が共通に生成する疑似ランダム信号の一例を示しており、図４（ｂ）は、第１発震器１ａ、第２発震器１ｂ、第３発震器１ｃ、及び第４発震器１ｄのそれぞれの乱数シフト部１２が出力する疑似ランダム信号Ｓ１〜Ｓ４を示している。

詳細は後述するが、震動波発生部１３は、疑似ランダム信号の値が１のとき正の振幅、０のとき負の振幅となるような震動波を、１ビット当たり０．１［ｓ］＝１００［ｍｓ］の間隔で生成する。

実施の形態に係る信号処理システム１００において、発震器１が震動波を送出してから受信器２に至るまでに要する時間として定義される「遅延時間」における最大値は、１０［ｓ］と見込まれている。このため、所定の時間Ｔも１０［ｓ］に設定されており、結果として所定の長さＤは１００ビット（１０［ｓ］／０．１［ｓ］＝１００）となる。

実施の形態に係る信号処理システム１００が備える発震器１の数は４つであるから、式（１）より所定の周期Ｐは、
Ｐ≧（４＋１）×１００＝５００ （２）
を満たす値として設定される。

図４（ａ）において、記号Ａ〜Ｅで示される矩形内には、それぞれ１００ビットの長さの信号、すなわち所定の長さＤと同じ長さの信号が含まれる。また記号Ｆで示される矩形内には、１１ビットの信号が含まれる。発震器１が生成する信号は、記号Ａ〜Ｆで示される矩形内に含まれる信号をつなげた５１１ビットの長さを持つ信号である。震動波発生部１３は１ビット当たり０．１［ｓ］の間隔で震動波を生成するため、発震器１が生成する信号は５１．１［ｓ］の長さの信号となる。

ここで、疑似乱数生成部１１が発生するＭ系列の周期は５１１ビット以上に設定されている。このため、第１発震器１ａが生成する信号の相関は、ラグタイム０において１であり、それ以外はほぼ０となる。したがって、第１発震器１ａが生成する信号は式（２）を満たす。

以下説明の便宜のため、記号Ａ〜Ｆで示される矩形内に含まれる信号を、それぞれ信号Ａ〜Ｆと記載することがある。また、疑似乱数生成部１１が生成し、乱数シフト部１２が出力する疑似ランダム信号の長さである５１１ビットを、信号長Ｌと記載することがある。

各発震器１における乱数シフト部１２の図示しない記憶部には、上述したＭの値があらかじめ設定されている。例えば第１発震器１ａの乱数シフト部１２に設定されたＭの値は０であるため、第１発震器１ａの乱数シフト部１２は疑似乱数生成部１１が生成した疑似ランダム信号をシフトしない。結果として第１発震器１ａの乱数シフト部１２が出力する疑似ランダム信号Ｓ１は、疑似乱数生成部１１が生成した疑似ランダム信号と同一となり、図４（ｂ）に示すように信号Ａ〜Ｆがこの順でつながった信号となる。

第２発震器１ｂの乱数シフト部１２には、Ｍ＝１が設定されている。このため第２発震器１ｂの乱数シフト部１２は、疑似乱数生成部１１が生成した疑似ランダム信号を、所定の長さＤだけ循環的にシフトして得られた信号を出力する。このため図４（ｂ）に示すように、第２発震器１ｂの乱数シフト部１２が出力する疑似ランダム信号Ｓ２は、信号Ｂ〜Ｆ、及び信号Ａがこの順でつながった信号となる。

同様に、第３発震器１ｃの乱数シフト部１２及び第４発震器１ｄの乱数シフト部１２には、それぞれＭ＝２及び３が設定されている。このため、第３発震器１ｃの乱数シフト部１２が出力する疑似ランダム信号Ｓ２は、信号Ｃ〜Ｆ、及び信号Ａ、Ｂがこの順でつながった信号となり、第４発震器１ｄの乱数シフト部１２が出力する疑似ランダム信号Ｓ２は、信号Ｄ〜Ｆ、及び信号Ａ〜Ｃがこの順でつながった信号となる。

このように、疑似ランダム信号Ｓ１〜Ｓ４は、信号長Ｌ以上の周期を持つ同一の疑似ランダム信号を少なくとも所定の長さＤだけ循環的にシフトして生成された信号である。したがって、疑似ランダム信号Ｓ１〜Ｓ４は、少なくとも所定の長さＤの範囲内において互いに無相関な疑似ランダム信号となっている。

図５（ａ）−（ｃ）は、実施の形態に係る震動波発生部１３が生成する震動波を説明するための図である。震動波発生部１３は、乱数シフト部１２が出力した疑似ランダム信号に基づいて、図示しない油圧駆動のアクチュエータを駆動し、震動波を生成する。

図５（ａ）は、信号の値が「１」の場合における、アクチュエータの送油量とアクチュエータの変位との関係を示す図である。図５（ａ）において横軸は時間の経過を示し、縦軸はアクチュエータの変位又は送油量を示す。また図５（ａ）において破線はアクチュエータを駆動するための送油量の時間変動を示し、実線はアクチュエータの変位の時間変動を示す。

震動波発生部１３は、時刻０［ｍｓ］の時点でアクチュエータに送油を開始し、時刻２５［ｍｓ］に至るまでは送油量を増加する。震動波発生部１３は時刻２５［ｍｓ］を過ぎると送油量を減少させ、時刻５０［ｍｓ］の時点で送油量を０とする。その後、震動波発生部１３は時刻７５［ｍｓ］に至るまで送油量の減少を継続し、時刻７５［ｍｓ］を過ぎた後は再び送油量を増加させる。時刻１００［ｍｓ］において、震動波発生部１３は送油量を０とする。以上が信号の値が「１」の場合において震動波発生部１３が実行するアクチュエータの送油量とアクチュエータの変位の制御である。震動波発生部１３は、信号１ビットにつき１００［ｍｓ］の時間でアクチュエータを駆動する。

ここで、アクチュエータの変位は、アクチュエータを駆動するための送油量を時間積分した値となる。このため図５（ａ）に示すように、震動波発生部１３は、時刻５０［ｍｓ］で変位が最大となり、時刻０［ｍｓ］及び１００［ｍｓ］では変位の速度が０（すなわち、図５（ａ）における実線の傾きが０）となる。このように、震動波発生部１３がアクチュエータの変位の開始時と終了時とにおいて変位の速度を０とすることにより、アクチュエータにかかる負荷を低減し、バイブロサイスの故障率を低減することができる。

図５（ｂ）は、第１発震器１ａの乱数シフト部１２が生成した疑似ランダム信号Ｓ１の最初の５ビットに基づいて震動波発生部１３が生成した震動波の概形を模式的に示す図である。疑似ランダム信号Ｓ１の最初の５ビットは「１０００１」であるため、図５（ｂ）に示すように、震動波発生部１３は正の振幅の震動波を発生させた後に、３回続けて負の振幅の震動波を発生させ、その後再び正の振幅の震動波を発生させている。このように、震動波発生部１３は疑似ランダム信号に基づいてアクチュエータを駆動して震動波を発生させることにより、離散信号である疑似ランダム信号を連続的な信号である震動波に変換する。

図５（ｃ）は、震動波発生部１３が疑似ランダム信号に基づいて発生させた震動波の自己相関関数の概形を示す図である。図５（ｃ）に示すように、震動波発生部１３が発生させる震動波は良好な自己相関性を有している。

震動波記憶部１５は、震動波発生部１３が発生させた震動波をサンプリングして震動波の時系列データを生成して記憶する。震動波記憶部１５は、既知のＡＤコンバータ及び記憶装置によって実現される。

より具体的には、震動波記憶部１５は、震動波発生部１３が発生させた震動波を６．２５［ｍｓ］のサンプリングインターバルでサンプリングし、６４ビットの浮動小数点形式で時系列的に蓄積する。震動波発生部１３は乱数シフト部１２が出力する疑似ランダム信号を１ビットあたり１００［ｍｓ］の長さで震動波を生成する。このため、結果として震動波記憶部１５が記憶するデータは５１１×（１００／６．２５）＝５１１×１６＝８１７６個の要素を持つデータとなる。以下本明細書において、震動波発生部１３が記憶した時系列データを、送信信号配列Ｓ（ｊ）、（ｊ＝０，・・・，８１７５）と記載することがある。

図６（ａ）−（ｅ）は、実施の形態に係る発震器１が送出し、媒質Ｍを伝搬して受信器２が受信した震動波の波形を模式的に示す図である。発震器１の震動波発生部１３が生成した震動波は、発震器１と受信器２との間に存在する媒質Ｍを伝搬する距離及び媒質Ｍにおける震動波の速度に応じて遅延及び減衰して受信器２に到達する。なお図６（ａ）−（ｅ）は、シミュレーションにより生成した振動波の波形を示す図である。以下、図７及び図９に示す図も同様に、シミュレーションによる計算結果を示す。

図６（ａ）−（ｄ）は、それぞれ第１発震器１ａ〜第４発震器１ｄが送出した震動波における受信波形Ｓｒ１〜Ｓｒ４を示している。図６（ａ）−（ｄ）において、横軸は時間、縦軸は受信器２に到達したときの震動波の振幅を示している。図６（ａ）−（ｄ）において、τ１〜τ４はそれぞれ第１発震器１ａ〜第４発震器１ｄが媒質Ｍに送出した震動波が受信器２に到達するまでの時間、すなわち上述した遅延時間を示している。

図２に示すように、各発震器１が送出した震動波が受信器２に到達するまでに媒質Ｍを伝搬する経路は互いに異なる。図２に示す例では、各発震器１が送出した震動波が受信器２に到達するまでに媒質Ｍを伝搬する経路の長さは、第１発震器１ａ、第２発震器１ｂ、第３発震器１ｃ、及び第４発震器１ｄの順に短い。このため、遅延時間もτ１、τ２、τ３、及びτ４の順に短くなっている。また震動波は媒質Ｍを伝搬するうちに減衰する。このため、受信器２に到達したときの震動波の受信波形Ｓｒ１〜Ｓｒ４の振幅はこの順に小さくなっている。

各発震器１は同期しており、各発震器１は媒質Ｍ中に震動波を同時に送出する。このため受信器２は、図６（ａ）−（ｄ）に示す震動波を個別に受信するのではなく、各発震器１のそれぞれが送出して媒質Ｍを伝搬した震動波の合波である合成震動波を受信する。

図６（ｅ）は、受信器２が受信した合成震動波の波形を示す図である。図６（ｅ）に示すように、受信器２が受信する合成震動波は図６（ａ）−（ｄ）に示す受信波形Ｓｒ１〜Ｓｒ４の重ね合わせであり、一見しただけでは図６（ａ）−（ｄ）に示す震動波を分別することはできない。信号処理装置３は、疑似ランダム信号Ｓ１〜Ｓ４が備える統計的な性質を利用して、合成震動波から遅延時間もτ１〜τ４と各震動波の振幅を示す情報を抽出する。

図３の説明に戻る。信号処理装置３は、受信部３１、時系列データ生成部３２、震動波再生部３３、クロック生成部３４、修正データ生成部３５、相関算出部３６、及び遅延時間算出部３７を備える。
クロック生成部３４は、クロック生成部１４と同等の機能を有しており、衛星Ｓから受信した電波に基づいてクロックを生成する。クロック生成部３４は、生成したクロックを受信部３１に出力する。受信部３１は、クロック生成部３４から取得したクロックに同期して受信器２が生成した合成震動波の受信信号を受信する。

ここで、受信部３１が受信器２から受信信号を受信することを含めて信号の受信と考える。この場合、各発震器１のクロック生成部１４と受信器２のクロック生成部３４とは、いずれも衛星Ｓの時刻情報に同期しているため、共通の時間軸でクロック信号を出力できる。結果として、各発震器１が震動波を送出するのと同時に、信号処理装置３は信号の受信を開始することができる。

発震器１が震動波の送出を開始してから受信器２に到達するまではいくらかの時間を要するため、受信器２が受信を開始した直後にはまだ震動波は計測できない。そのため、クロック生成部３４が受信器２からの信号の受信を開始し、震動波の送出時間である５１．１［ｓ］受信を継続した時点では、発震器１が送出した震動波をすべて受信しきることができない。

ここで最大遅延時間は１０［ｓ］と見込まれているため、クロック生成部３４は、受信器２からの信号の受信を開始してから５１．１［ｓ］＋１０［ｓ］＝６１．１［ｓ］＝６１１００［ｍｓ］受信を継続すれば、発震器１が送出した震動波をすべて受信しきることができる。

時系列データ生成部３２は、受信部３１が受信した合成震動波の時系列データを生成する。時系列データ生成部３２は既知のＡＤコンバータであり、２５［ｍｓ］のサンプリングインターバルで合成震動波をデジタルデータに変換する。結果として、時系列データ生成部３２が生成するデータは６１１００［ｍｓ］／２５［ｍｓ］＝２４４４個の要素を持つデータとなる。このデータは最終的に６４ビットの浮動小数点形式のデータに変換され、６４ビット×２４４４レコードのデータとなる。以下本明細書において、時系列データ生成部３２が生成した時系列データを、受信信号配列Ｒ（ｊ）、（ｊ＝０，・・・，２４４３）と記載することがある。

修正データ生成部３５は、時系列データ生成部３２が生成した受信信号配列Ｒ（ｊ）の末尾から所定の長さＤに相当する長さの信号を切り出す。修正データ生成部３５は、受信信号配列Ｒ（ｊ）の先頭から所定の長さＤに相当する長さまでの時系列データに、切り出した信号を加算した修正時系列データＣｄ（ｊ）、（ｊ＝０，・・・，２０４３）を生成する。

ここで、疑似乱数生成部１１が生成する疑似ランダム信号において所定の長さＤは１００ビットに相当する。時系列データ生成部３２は震動波発生部１３が１ビット分の震動波を生成する間に４回サンプリングするため、修正時系列データＣｄ（ｊ）において所定の長さＤは４００個のデータ数に相当することになる。以上より、修正時系列データＣｄ（ｊ）は、以下の式（３）で表現できる。

図７（ａ）−（ｂ）は、実施の形態に係る修正データ生成部３５による修正時系列データＣｄ（ｊ）の生成処理を説明するための図である。図７（ａ）は、修正データ生成部３５による受信信号配列Ｒ（ｊ）の修正処理を説明するための図である。また図７（ｂ）は、修正データ生成部３５が修正処理の結果生成した修正時系列データＣｄ（ｊ）を示す図である。図７に示すように、修正データ生成部３５は受信信号配列Ｒ（ｊ）の末尾の４００個分のデータを受信信号配列Ｒ（ｊ）の先頭に加算した修正時系列データＣｄ（ｊ）を生成する。最終的に、修正時系列データＣｄ（ｊ）の要素数は２０４４個となる。

図８（ａ）―（ｂ）は、第１発震器１ａが送出した震動波の受信波形Ｓｒ１に関する修正時系列データＣｄ（ｊ）を説明するための模式図である。図８（ａ）に示すように、受信波形Ｓｒ１は、クロック生成部３４が受信信号配列Ｒ（ｊ）の生成を開始してから遅延時間τ１経過した後に、受信器２に到達する。このため、クロック生成部３４が受信信号配列Ｒ（ｊ）の生成を開始してから遅延時間τ１が経過するまでは、受信波形Ｓｒ１には第１発震器１ａが送出した震動波に由来する信号は含まれない。修正データ生成部３５が受信信号配列Ｒ（ｊ）の末尾の４００個分のデータを受信信号配列Ｒ（ｊ）の先頭に加算することにより、修正時系列データＣｄ（ｊ）は、全域にわたって第１発震器１ａが送出した震動波に由来する信号を含むことになる。

発震器１が送出した震動波は、媒質Ｍを伝搬しても基本的には震動波発生部１３が発生させた震動波の形状を保っている。このため、図８（ｂ）に示すように、修正時系列データＣｄ（ｊ）と震動波の波形を表す送信信号配列Ｓ（ｊ）とを相対的に循環シフトしながら相関を算出することにより、遅延時間τを求めることができる。

図３の説明に戻る。震動波再生部３３は、震動波記憶部１５が記憶した送信信号配列Ｓ（ｊ）を取得する。震動波再生部３３は、例えば震動波記憶部１５から有線通信又は無線通信によって送信信号配列Ｓ（ｊ）を取得してもよいし、持ち運び可能な記録媒体にコピーされた送信信号配列Ｓ（ｊ）を記録媒体から読み出すことによって取得してもよい。

相関算出部３６は、Ｎ個の発震器１がそれぞれ送出したＮ個の震動波の送信信号配列Ｓ（ｊ）それぞれと、修正時系列データＣｄ（ｊ）との相関を示すデータを算出する。より具体的には、相関算出部３６は、修正時系列データＣｄ（ｊ）を０から所定の長さＤにＮを乗じた長さに至るまで、循環的にシフトしながら、順番に送信信号配列Ｓ（ｊ）との自己相関係数を算出する。すなわち、相関算出部３６は修正時系列データＣｄ（ｊ）送信信号配列Ｓ（ｊ）との自己相関係数の時間変化を算出する。

より具体的には、相関算出部３６は、遅延時間が０のときにおける自己相関係数ρ（０）を、以下の式（４）にしたがって算出する。

式（４）において送信信号配列Ｓ（ｊ）のインデックスに４が乗じられているのは、送信信号配列Ｓ（ｊ）は修正時系列データＣｄ（ｊ）の４倍のサンプリングレート（すなわち、サンプリング間隔は４分の１）でサンプリングされているからである。また、Ｃａｖは修正時系列データＣｄ（ｊ）の時間平均を表し、Ｓａｖは送信信号配列Ｓ（ｊ）の時間平均を表す。すなわち、Ｃａｖ及びＳａｖは、それぞれ以下の式（５）及び式（６）にしたがって算出される。

続いて相関算出部３６は、送信信号配列Ｓ（ｊ）を以下の式（７−１）及び式（７−２）にしたがって１データ分だけ循環的にシフトする。
Ｓ（０）←Ｓ（８１７５） （７−１）
Ｓ（ｊ）←Ｓ（ｊ−１）、（ｊ＝１，・・・，８１７５） （７−２）

相関算出部３６は、式（７−１）及び式（７−２）にしたがって循環的にシフトされた送信信号配列Ｓ（ｊ）を式（４）に代入することにより、送信信号配列Ｓ（ｊ）の１サンプリング分である６．２５［ｍｓ］の遅延時間における自己相関係数ρ（１）を算出できる。

実施の形態に係る信号処理システム１００において、発震器１が震動波を送出してから受信器２に至るまでに要する時間として定義される遅延時間における最大値は１０［ｓ］＝１００００［ｍｓ］と見込まれている。このため、送信信号配列Ｓ（ｊ）においては、１００００［ｍｓ］／６．２５［ｍｓ］＝１６００個分のデータが所定の長さＤに相当する。

相関算出部３６は式（７−１）及び式（７−２）にしたがって送信信号配列Ｓ（ｊ）を１データずつ循環的にシフトしながら式（４）の自己相関係数ρを１６００回算出し、ρ（０）〜ρ（１５９９）を得る。ρ（０）〜ρ（１５９９）が、発震器１の送信開始から所定の長さＤに対応する時間までの自己相関係数ρの時間変化である。相関算出部３６が算出した自己相関係数ρ（０）〜ρ（１５９９）のうち、第１発震器１ａが送出した震動波が受信器２に到達した時間である遅延時間τ１に対応する自己相関係数ρにおいて、自己相関係数ρは急峻なピークを発生する。

図４を参照して説明したように、第１発震器１ａの乱数シフト部１２が出力した疑似ランダム信号Ｓ１を所定の長さＤだけ循環的にシフトすると、第２発震器１ｂの乱数シフト部１２が出力した疑似ランダム信号Ｓ２と一致する。同様に、第２発震器１ｂの乱数シフト部１２が出力した疑似ランダム信号Ｓ２を所定の長さＤだけ循環的にシフトすると、第３発震器１ｃの乱数シフト部１２が出力した疑似ランダム信号Ｓ３と一致する。また第３発震器１ｃの乱数シフト部１２が出力した疑似ランダム信号Ｓ３を所定の長さＤだけ循環的にシフトすると、第４発震器１ｄの乱数シフト部１２が出力した疑似ランダム信号Ｓ４と一致する。

このことは、相関算出部３６が式（７−１）及び式（７−２）にしたがって送信信号配列Ｓ（ｊ）を１６００回循環的にシフトして得られる送信信号配列Ｓ（ｊ）は、第２発震器１ｂの震動波発生部１３が送出した震動波の送信信号配列Ｓ（ｊ）となることを意味する。このため相関算出部３６が式（７−１）、式（７−２）、及び式（４）にしたがって自己相関係数ρ（１６００）〜ρ（３１９９）を算出することにより、第２発震器１ｂが送出した震動波における送信信号配列Ｓ（ｊ）と、修正時系列データＣｄ（ｊ）との自己相関係数ρの時間変化を算出することができる。

同様に、相関算出部３６は自己相関係数ρ（３２００）〜ρ（４７９９）を算出することにより、第３発震器１ｃが送出した震動波における送信信号配列Ｓ（ｊ）と、修正時系列データＣｄ（ｊ）との自己相関係数ρの時間変化を算出することができる。また、相関算出部３６は自己相関係数ρ（４８００）〜ρ（６３９９）を算出することにより、第４発震器１ｄが送出した震動波における送信信号配列Ｓ（ｊ）と、修正時系列データＣｄ（ｊ）との自己相関係数ρの時間変化を算出することができる。結果として、相関算出部３６は、式（７−１）、式（７−２）、及び式（４）にしたがって自己相関係数ρ（０）〜ρ（６３９９）を計算することにより、第１発震器１ａ〜第４発震器１ｄのそれぞれが送出した震動波における送信信号配列Ｓ（ｊ）と、修正時系列データＣｄ（ｊ）との自己相関係数ρの時間変化を一度に算出することができる。

図９は、実施の形態に係る相関算出部３６が算出した自己相関係数ρの時間変化のグラフを例示する図である。図９に示すように、自己相関係数ρの時間変化は、横軸の値がＴ１、Ｔ２、Ｔ３及びＴ４の位置において急峻なピークが認められる。これは、第１発震器１ａ〜第４発震器１ｄがそれぞれ送出した震動波の遅延時間であるτ１〜τ４に対応する位置である。

図９において、Ｔ１は、第１発震器１ａが送出した震動波の遅延時間τ１そのものを示す。Ｔ２は、第２発震器１ｂが送出した震動波の遅延時間τ２に、送信信号配列Ｓ（ｊ）における１６００個分のデータに対応する時間、すなわち所定の長さＤに対応する所定の時間Ｔ＝１０［ｓ］のバイアスが乗っている。このため、Ｔ２＝τ２＋Ｔである。また、Ｔ３は、第３発震器１ｃが送出した震動波の遅延時間τ３に所定の長さＤの２倍に対応する２０［ｓ］のバイアスが乗っており、Ｔ３＝τ３＋２Ｔとなっている。同様にＴ４＝τ４＋３Ｔである。

遅延時間算出部３７は、相関算出部３６が算出したデータである自己相関係数ρを解析してピークを特定し、各発震器１がそれぞれ震動波を送出してから受信器２に到達するまでに要した時間であるそれぞれの震動波の遅延時間を算出する。これにより、各発震器１が送出した震動波の伝搬経路となった媒質Ｍにおけるそれぞれの震動波の速度が推定できる。信号処理システム１００のユーザは、推定された震動波の速度に基づいて、震動波が伝搬した媒質Ｍの特性を推定することができる。

＜比較例＞
図１０は、比較例に係る自己相関係数の算出方法を説明するための図である。比較例に係る自己相関係数の算出方法においては、受信器２が受信した合成震動波の末尾から所定の長さＤの信号を切り出し、その信号を合成震動波の先頭に加算して生成される信号を相関算出の対象とすることはしない。その代わり、受信器２が受信した合成震動波を修正処理することなく、自己相関係数を算出する。

図１０に示すように、第１発震器１ａが送出した震動波は、遅延時間τ１において受信器２に到達したとする。この場合、時系列データ生成部３２が受信波形Ｓｒ１の生成を開始してから遅延時間τ１経過した後に、第１発震器１ａが送出した震動波に由来する信号が到達することになる。

比較例に係る自己相関係数の算出方法においては、送信信号配列Ｓ（ｊ）を単にシフトさせながら受信波形Ｓｒ１との相関係数を算出する。この場合、図１０に示すように、送信信号配列Ｓ（ｊ）を遅延時間τ１に相当する量だけシフトした場合には、受信波形Ｓｒ１と送信信号配列Ｓ（ｊ）とをそれぞれ構成する信号が一致する。しかしながら、それ以外の場合には、送信信号配列Ｓ（ｊ）が相関係数を算出する対象となる信号には、Ｓ（ｊ）の情報が一部含まれていないため、相関係数が零となるとは限らない。このため、算出した相関係数の時間変化には、遅延時間を示す各ピークのほか、微細なダミーのピークが発生しうる。

図１１は、比較例に係る自己相関係数の算出方法を用いて算出した自己相関係数の時間変化のグラフを例示する図である。図１１は、図９に示す実施の形態に係る相関算出部３６が算出した自己相関係数ρの時間変化のグラフと比較すると、遅延時間を示す各ピークのほか微細なダミーのピークが発生していることを示している。

図１２は、地震探査における自己相関関数の実測データの例を示すグラフである。地中を伝搬する信号は、途中の地下構造を反映し、位相の反転も含む複数の遅延時間を示すピークが含まれる。図１２において、破線で示す楕円Ｅ１に出現している微小なピークは地震探査で特に重視されるいわゆる初動で、受信器に最初に到達する振動を示すピークである。破線で示す楕円Ｅ２に出現している各ピークは、直接到達する振動のうち、地中のインピーダンス界面で屈折または反射により遅延した振動、さらに、楕円Ｅ３に出現している反射を複数回繰り返して到達した信号を示す各ピークが続く。

このように、石油等の探査における地質構造調査、あるいは土木分野における地盤調査や既存構造物内部の損傷調査等において震動波の初動のような小さなピークを検出することは重要である。また、地震探査分野で活用が期待されるフルウェーブインバージョンなどの手法では微細なピークの検知がより重視される。実施の形態に係る自己相関係数の算出方法はこの初動も含め、全てのピークが出現する領域においても遅延時間を示すピーク以外のノイズが抑制されるため、比較例に係る自己相関係数の算出方法と比べてピークの検出精度を高めることもできる。

＜実施の形態の効果＞
以上説明したように、実施の形態に係る信号処理システム１００によれば、複数の発震器から同時に送出した信号を受信後に弁別する際の弁別精度を高めることができる。特に地震探査で重視される初動（送信器から受信機に最初に到達する振動）はそれに続く振動より、極めて小さいことが多いが、本発明は、サイドローブが極めて小さいため、はこのような初動も含め、微少な信号を容易に弁別することが可能になる。

以上、本発明を実施の形態を用いて説明したが、本発明の技術的範囲は上記実施の形態に記載の範囲には限定されない。上記実施の形態に、例えば、自己相関係数の算出に、タイムドメインによる処理に替え、周波数ドメインによる解析を使用するなど、多様な変更又は改良を加えることが可能であることが当業者に明らかである。そのような変更又は改良を加えた形態も本発明の技術的範囲に含まれ得ることが、特許請求の範囲の記載から明らかである。

＜変形例＞
上記の説明では、信号処理システム１００が１つの受信器２を含む場合について説明したが、受信器２は複数あってもよい。地震探査等においては、複数地点において測定を実施するのが一般的である。このため、異なる地点に複数の受信器２を展開して測定を実施する。仮に受信器２の数が足りなければ、受信器２を移動して測定を繰り返す。複数の受信器２を用いることにより、測定時間を短縮することができる。

上記では、信号長Ｌやサンプリング間隔１００［ｍｓ］、２５［ｍｓ］、６．２５［ｍｓ］等、具体的な数値を用いて説明したが、これらに限定されるものではない。これらの具体的な値は、信号処理システム１００が測定の対象とする媒質Ｍの事前情報や実際の測定結果等を考慮して、実験により適宜定めればよい。

１・・・発震器
２・・・受信器
３・・・信号処理装置
１１・・・疑似乱数生成部
１２・・・乱数シフト部
１３・・・震動波発生部
１４・・・クロック生成部
１５・・・震動波記憶部
３１・・・受信部
３２・・・時系列データ生成部
３３・・・震動波再生部
３４・・・クロック生成部
３５・・・修正データ生成部
３６・・・相関算出部
３７・・・遅延時間算出部
１００・・・信号処理システム

Claims (10)

1. 媒質に震動波を送出するＮ個（Ｎは１以上の整数）の発震器と、
   前記Ｎ個の発震器のそれぞれが送出して前記媒質を伝搬した震動波の合波を受信する受信器と、
   前記受信器が受信した合波の時系列データの末尾から所定の長さＤのデータを切り出して、前記時系列データの先頭から前記所定の長さＤまでの時系列データに加算した修正時系列データを生成する修正データ生成部と、
   前記Ｎ個の発震器がそれぞれ送出したＮ個の震動波の時系列データのそれぞれと、前記修正時系列データとの相関を示すデータを算出する相関算出部と、
   を備え、
   前記所定の長さＤは、前記Ｎ個の発震器がそれぞれ前記震動波を送出してから前記受信器に到達するまでに要した時間である遅延時間の最大値の時間が経過するまでの間に前記発震器が送出する震動波の長さであることを特徴とする信号処理システム。
2. 前記Ｎ個の発震器のそれぞれは、前記所定の長さＤの範囲内において互いに無相関な疑似ランダム信号に基づく震動波を送出する、
   ことを特徴とする請求項１に記載の信号処理システム。
3. 前記Ｎ個の発震器のそれぞれは互いに同期して前記震動波を送出し、
   前記受信器は、前記Ｎ個の発震器のそれぞれによる前記震動波の送出と同期して前記時系列データの受信を開始する、
   ことを特徴とする請求項１又は２に記載の信号処理システム。
4. 前記Ｎ個の発震器のそれぞれは、前記所定の長さＤにＮ＋１を乗じた長さ以上の周期を有する疑似ランダム信号を、前記所定の長さのＭ倍（Ｍは０以上Ｎ−１以下の整数であって前記Ｎ個の発震器のそれぞれで異なる）の長さだけ循環的にシフトして得られた信号に基づく震動波を送出する、
   ことを特徴とする請求項１から３のいずれか一項に記載の信号処理システム。
5. 前記相関算出部は、前記修正時系列データを０から前記所定の長さＤにＮを乗じた長さに至るまで循環的にシフトしながら、順番に前記疑似ランダム信号との自己相関係数を算出する、
   ことを特徴とする請求項４に記載の信号処理システム。
6. 前記相関算出部が算出したデータを解析して、前記Ｎ個の発震器それぞれの前記遅延時間を算出する遅延時間算出部をさらに備える、
   ことを特徴とする請求項１から５のいずれか一項に記載の信号処理システム。
7. 前記Ｎ個の発震器はそれぞれ、前記震動波それぞれの遅延時間の最大値の時間が経過するまでの間に送出する震動波の長さが、前記所定の長さＤ以下となるように、前記震動波を送出する、
   ことを特徴とする請求項６に記載の信号処理システム。
8. プロセッサが実行する信号処理方法であって、
   Ｎ個（Ｎは１以上の整数）の発震器のそれぞれで媒質に震動波を送出するステップと、
   前記媒質を伝搬した前記震動波の合波を受信器で受信して受信波の時系列データを生成するステップと、
   前記時系列データの末尾から所定の長さＤのデータを切り出して、前記時系列データの先頭から前記所定の長さＤまでの時系列データに加算した修正時系列データを生成するステップと、
   送出した前記Ｎ個の震動波の時系列データのそれぞれと、前記修正時系列データとの相関を示すデータを算出するステップと、
   を含み、
   前記所定の長さＤは、前記Ｎ個の震動波それぞれが送出されてから前記受信器に到達するまでに要した時間である遅延時間の最大値の時間が経過するまでの間に前記発震器が送出する震動波の長さであることを特徴とする信号処理方法。
9. 媒質に震動波を送出するＮ個（Ｎは１以上の整数）の発震器のそれぞれが送出して前記媒質を伝搬した震動波の合波を受信器で受信して得られた受信波の時系列データを生成する時系列データ生成部と、
   前記時系列データの末尾から所定の長さＤのデータを切り出して、前記時系列データの先頭から前記所定の長さＤまでの時系列データに加算した修正時系列データを生成する修正データ生成部と、
   前記Ｎ個の発震器がそれぞれ送出したＮ個の震動波の時系列データのそれぞれと、前記修正時系列データとの相関を示すデータを算出する相関算出部と、
   を備え、
   前記所定の長さＤは、前記Ｎ個の発震器がそれぞれ前記震動波を送出してから前記受信器に到達するまでに要した時間である遅延時間の最大値の時間が経過するまでの間に前記発震器が送出する震動波の長さであることを特徴とする信号処理装置。
10. コンピュータに、
    媒質に震動波を送出するＮ個（Ｎは１以上の整数）の発震器のそれぞれが送出して前記媒質を伝搬した震動波の合波を受信器で受信して得られた受信波の時系列データを生成する機能と、
    前記時系列データの末尾から所定の長さＤのデータを切り出して、前記時系列データの先頭から前記所定の長さＤまでの時系列データに加算した修正時系列データを生成する機能と、
    前記Ｎ個の発震器がそれぞれ送出したＮ個の震動波の時系列データのそれぞれと、前記修正時系列データとの相関を示すデータを算出する機能と、
    を実現させ、
    前記所定の長さＤは、前記Ｎ個の発震器がそれぞれ前記震動波を送出してから前記受信器に到達するまでに要した時間である遅延時間の最大値の時間が経過するまでの間に前記発震器が送出する震動波の長さであることを特徴とするプログラム。