地中レーダ画像処理法
【課題】 管路等の人工物を地中レーダなどにより探査した場合の画像を処理して測定現場で目的の人工物を容易に判定できる効率的な処理法を提供する。
【解決手段】 各測定位置における受信信号をデコンボリューション等により処理して反射信号の遅延時間を推定し、推定された遅延時間と測定位置の差が所定の範囲内にあるものを同じグループに帰属させることにより推定遅延時間をグループ化し、グループごとに同一グループの画素の中で最も遅延時間が短い位置の近傍に反射点があると仮定して、同グループの画素による反射点の推定を行い、推定された反射点に同グループに属する遅延時間の評価値と推定誤差によって決まる評価値を与え、近距離にある推定位置を統合することにより、同一対象物による可能性が高い推定位置を統合して、推定位置の確度を高める。
【解決手段】 各測定位置における受信信号をデコンボリューション等により処理して反射信号の遅延時間を推定し、推定された遅延時間と測定位置の差が所定の範囲内にあるものを同じグループに帰属させることにより推定遅延時間をグループ化し、グループごとに同一グループの画素の中で最も遅延時間が短い位置の近傍に反射点があると仮定して、同グループの画素による反射点の推定を行い、推定された反射点に同グループに属する遅延時間の評価値と推定誤差によって決まる評価値を与え、近距離にある推定位置を統合することにより、同一対象物による可能性が高い推定位置を統合して、推定位置の確度を高める。
【発明の詳細な説明】
【技術分野】
【０００１】
管路等の人工物を地中レーダなどにより探査した場合の画像を処理して測定現場で目的の人工物を容易に判定できる効率的な処理法を提供する。
【背景技術】
【０００２】
地中レーダ等の近距離レーダの画像を処理する方法に関する。特に、人工物など周囲との境界が明確な対象からの反射信号を特定して対象の判定を容易にする方法に関する。本方法は地中レーダだけでなく、対象物の像を得るのに開口合成的な手法を必要とする音波探査および地震探査にも適用可能である。
【０００３】
【特許文献１】特開２００４−１９８１９５号公報
【特許文献２】特開２００１−０３３５６４号公報
【特許文献３】特開２００２−１０７４４９号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
地中レーダは指向性が低いので対象からの反射を分離するための処理が必要である。
現在このために使用されている手法は凡そ以下のように分類される。
・マイグレーション
開口合成的な方法である。一般に、複数の対象物の反射が重なり合っている場合には対象物の無い場所に虚像を生じる欠点がある。計算方法には遅延時間を計算する方法（ディフラクションスタックマイグレーション）とFFTを用いる方法がある。前者は格子点ごとに計算処理を行なう必要があるために多大な量の計算を必要とし、調査結果を現場で処理するような用途には適していない。一方、後者は計算処理に要する時間は少ないが、虚像を発生しやすいという欠点がある。
特に、地中レーダの場合は測定対象とする空間の大きさに比べてレーダ信号の持続時間が長いという特殊な事情により、これらの方法には対象物の明瞭な位置を得るのが難しいという欠点がある。
・トモグラフィ
計算量が極めて多大であるために、現場での処理にはまったく適していない。
以上のように、現在使われている処理法は、多大な計算時間を要する、また、特に複数の対象物が近接して存在する場合に虚像を生じるという欠点がある。本発明はこれら２つの課題を解決することを目的とする。
【課題を解決するための手段】
【０００５】
図５は本発明の処理フローの一例を示す図である。
ステップ１０１：地中レーダの送信アンテナおよび受信アンテナの少なくとも一方を移動しながら反射信号を取得する。
ステップ１０２：デコンボリューション、MUSIC法等により処理して反射信号の遅延時間を推定する。
ステップ１０３：受信信号を孤立反射体からの反射信号の重なりとみなし、ステップ１０２で得られた遅延時間を最も適切に説明可能な孤立物体の位置を推定する。
ステップ１０４：推定された遅延時間を確からしさを基に統合して、遅延時間の推定誤差を軽減するとともに、確からしさの低い遅延時間を除外して偽の対象物が発生するのを抑制する。
ステップ１０５：推定された遅延時間と測定位置の差が所定の範囲内にあるものを同じグループに帰属させることにより推定遅延時間をグループ化し、グループごとに同一グループの画素の中で最も遅延時間が短い位置の近傍に反射点があると仮定して、同グループの画素による反射点の推定を行い、推定された反射点に同グループに属する遅延時間の評価値と推定誤差によって決まる評価値を与える。
ステップ１０６：評価値の低い推定位置を除外することにより偽の対象物が発生するのを抑制し、 近距離にある推定位置を統合することにより、同一対象物による可能性が高い推定位置を統合して、推定位置の確度を高める。
以下で各処理の詳細を述べる。
【０００６】
［測定方法］
解決法を説明するために、本発明が扱うデータの取得方法について地中レーダを例に説明する。
地中レーダでは、通常受信アンテナおよび送信アンテナの少なくとも一方を測定対象領域上の地表面で移動しながら送信アンテナから電磁波が発射された時刻に同期して受信アンテナで信号を受信する。地震波探査では受波器および送波器のいずれか一方を移動して測定する方法が多く行なわれるが、地中レーダでは送信アンテナと受信アンテナを対にして同時に移動しながら測定する場合が大部分である。
【０００７】
本発明による画像処理は
(1) 送信アンテナおよび受信アンテナの少なくとも一方を移動しながら反射信号を取得する。
(2) 各測定位置における受信信号をデコンボリューション等により処理して、反射信号の遅延時間を推定する。
(3) 受信信号を孤立反射体からの反射信号の重なりとみなし、(2)で得られた遅延時間を最も適切に説明可能な孤立物体の位置を推定する。
という手順で行う。
【０００８】
手順(2)で推定された遅延時間は一般に誤差および偽の推定時間を含む。そこで本発明では、
(A) 推定された遅延時間を確からしさを基に統合して、遅延時間の推定誤差を軽減するとともに、確からしさの低い遅延時間を除外して偽の対象物が発生するのを抑制する。
また、本発明の目的は孤立物体とみなせる対象物探索のための現場における迅速な処理であるので、手順(3)においては、
(B) 手順(2)で推定された遅延時間と測定位置の差が所定の範囲内にあるものを同じグループに帰属させることにより推定遅延時間をグループ化する。グループごとに、同一グループの画素の中で最も遅延時間が短い位置の近傍に反射点があると仮定して、同グループの画素による反射点の推定を行うことにより、処理時間を短縮する。このとき、推定された反射点に同グループに属する遅延時間の評価値と推定誤差によって決まる評価値を与える。
【０００９】
(C) 適切さの低い推定位置除外することにより偽の対象物が発生するのを抑制する。
(D) 近距離にある推定位置を統合することにより、同一対象物による可能性が高い推定位置を統合して、推定位置の確度を高める。
【発明の効果】
【００１０】
(1) 処理速度が速い。
(2) 対象物の位置と反射強度を特定できるために、連続物体の同定など高度な処理を行なうことが容易になる。
(3) 対象物の推定位置が明瞭になる。そのため、測定結果の解釈が容易になる。
【実施例】
【００１１】
［受信データ］
説明を簡単にするために、また、通常の地中レーダによる測定を考慮して、送信機XT（送信アンテナ）と受信機XR（受信アンテナ）を所定の相互位置関係に保持しながら移動して受信信号を受信する。ここで、送信機XTとは測定対象領域に電磁波、地震波等のセンシングを行なう信号を送出する手段を指す。また、受信機XRとは電磁波、地震波等センシングを行なう信号を検知する手段を指す。
説明のために座標を定義する。信号が送出される方向をｘ軸、送信機XT、受信機XRが移動する方向をｙ軸、これら2軸に直交して右手系をなす方向をz軸とする。
ｘ = 0 の面に沿って送信機XT、受信機XRを移動して測定を行なう場合を例にとって説明する。送信機XT、受信機XRを移動する面が平坦である限りこの仮定により一般性を失うことは無い。説明の都合上以後この面をセンシング面と呼ぶことにする。測定点はセンシング面上に並んでいるものとする。送信機XT、受信機XRの移動は直線に沿ってでも、曲線に沿ってでも構わない。
【００１２】
以下の説明では、アンテナを直線的に移動して2次元のデータを取得するものとする。また、送信アンテナのみを移動する場合、受信アンテナのみを移動する場合もフィルタに変更を加えるだけで全く同様に扱うことができる。
測定点がｙ座標方向にＮy個並んでいる。測定点の位置ＰMをｙ座標方向に番号付けして
【数１】
とする。ここで、z座標は省略した。また、上付き添え字Tは転置を表す。
送信アンテナと受信アンテナはｙ軸方向に並んでおり、一定の間隔を保って移動するから、それぞれの基準点（送信アンテナと受信アンテナの中点）から位置をｄT、ｄRとすると、測定点の位置ベクトルの組ＰMを用いて送信アンテナの位置ベクトルの組ＰTと受信アンテナの位置ベクトルの組ＰRは
【数２】
と表すことができる。
【００１３】
受信信号の長さはMであり、等間隔tsでサンプリングされているものとする。簡単のために送信アンテナから信号が送出された時刻を0とする。すると、受信信号を受信する受信時刻の組ｔMは信号を発射した時刻を始点とする相対値であり、
【数３】
となる。
測定点ＰMjおける受信信号sjを
【数４】
と表す。全受信データの組sは
【数５】
と表すことができる。
【００１４】
［画像空間］
測定対象領域内の対象物の位置および性状の推定値を表示するために画像空間を定義する。画像空間は測定対象領域内に設けられた矩形格子点から構成される。ｙ座標軸方向のj 番目の格子点のｙ座標ｙI jは受信データを取得した測定点のｙ座標である。つまり、
【数６ａ】
ｘ座標は、信号受信時刻tMの間に信号が往復する距離である。したがって、ｘ 軸方向のm番目の格子点のｘ座標ｘImは、時刻tMmを用いると
【数６ｂ】
となる。ここで、csignalは信号の伝搬速度である。
【００１５】
［反射の遅延時間と振幅の推定］
地中レーダが送・受信する信号は、送信アンテナと受信アンテナの過渡応答により拡がりを持っているパルス状信号であるために、個々の対象物からの反射信号は、レーダ画像上で無視できない拡がりを持っている。同一の測定点で周波数を所定の間隔で変えながら測定した周波数ごとの受信信号の振幅と位相である場合も、逆フーリエ変換を施せば同様なパルス状信号になるので、状況は同じである。
この拡がりは測定対象である空間の大きさの数10％に達するので、このような画像から対象物の位置を明瞭に推定することは難しい。そこで、本発明においては、反射の遅延時間と振幅から対象物の位置を推定する。一般に逆フィルタによるデコンボリューションは雑音が発生しやすい欠点がある。送・受信アンテナの過渡応答を正確に知ることが難しいこと、反射信号のスペクトルが対象物（反射体）の深さに依存することから、地中レーダの場合は良好なデコンボリューションを行うことは特に困難である。そこで、事後確率を最大にする方法(Jensen,J. A.,Leeman,S (1994): Non-parametric Estimation of Ultrasound Pulses,IEEE Trans. Biomed. Eng. 41,No. 10,929-936.)や、いわゆるMUSIC法など(S. Shrestha and I. Arai: Signal ＰRocessing of Ground Penetrating Radar Using Spectral Estimation Techniques to Estimate the Position of Buried Targets)を用いて反射信号の遅延時間と振幅を推定する。後者は前者に比べて処理時間が少ないものの、偽の反射を発生しやすい欠点がある。
そこで、本発明では処理時間が短い後者の方法により遅延時間を推定し、推定した遅延時間に付随する振幅を基に遅延時間の取捨選択と統合して偽の反射を低減する手法を考案している。
【００１６】
［MUSIC法による遅延時間推定］
以下処理方法について説明するが、その前に上記文献に基づいてMUSIC法による遅延時間の推定について説明する。
MUSIC法による遅延時間の推定は、受信信号Ｓj= [Sj,0 ，Sj,1 ，… Sj,M-1]Tの周波数スペクトルＳj= [Sj,0 ，Sj,1 ,… Sj,M-1]TをFFTなどにより計算する。次に、例えば振幅の絶対値が所定の値より大きいなど適当な周波数帯域中のスペクトル成分Ｓj'= [Sj,ms' ，Sj,ms+1' ，… Sj,me']Tを用いて遅延時間の推定を行う。
MUSIC法では、受信信号Ｓj'が真の信号
と雑音Ｎjから構成されていると考えて、
【数７】
とモデル化する。ここで、Ａは混合行列、つまり遅延を与える行列である。またLは反射信号（反射点）の個数である。受信信号Ｓj'の共分散行列ＲSSを求めて固有値分解する。
【数８】
行列Ｄは共分散行列ＲSSの固有値の２乗を対角要素とする行列である。ここで、絶対値が小さな固有値が雑音Ｎjによるもの、絶対値が大きな固有値を真の信号Ｓjによるものと仮定し、(8)式をさらに
【数９】
と分解する。ここで、ＤSは絶対値が大きな固有値、つまり、真の信号
に対応する固有値、ＵSはこれらの固有値に対応する固有ベクトルである。ＤNは絶対値が小さな固有値、つまり、雑音Ｎjに対応する固有値、ＵNはこれらの固有値に対応する固有ベクトルである。
【００１７】
受信信号Ｓj'がなす空間を、信号空間（真の信号に対応する固有ベクトルなす空間）と雑音空間（雑音に対応する固有ベクトルがなす空間）に分割する。信号空間のベクトルは雑音空間に直交するから、混合行列Ａの列ベクトルａj=[ａj,0 ，ａj,1 ，… ａj,L'-1]Tで、雑音空間に直交するものを求めると、遅延時間を得ることができる。列ベクトルａj,l (l = 0，… L'-1)を信号空間を張る固有ベクトルで展開して、展開係数とＤSから列ベクトルａj,l (l = 0，… L'-1)に付随する電力を得ることができる。
【００１８】
［偽反射の抑制法］
次に、本発明で用いる偽の反射を抑制する方法について説明する。
MUSIC法では、雑音に対応する固有値の個数、あるいは真の信号に対応する固有値の個数を指定する必要がある。地中レーダの受信信号では、反射点の個数は未知である。適当に余裕を持った個数の反射点があるもとして処理を行う。すると、偽の反射が生ずる。予め反射点の個数が分かっている場合でも、付随する電力が大きなものから与えれた個数だけ反射を選択した場合に、検知すべき反射が全てその中に含まれるとは限らず、個数に余裕を持って遅延時間を推定する必要がある場合が多い。
図６（ａ）は測定された信号を表し、図６（ｂ）および図６（ｃ）は図６（ａ）に示される測定信号を処理して得られた信号を表している。この場合、図５（ｂ）の黒線で表される真の反射点以外に、破線で表される多くの偽の反射点が得られる。図から分かるようにこれらの偽の反射点は真の反射点の周辺に生ずることが多い。真の反射点に近い位置に多くの反射点が推定されることに変わりは無いが、実際の位置は雑音の状態等により変化する。これらの反射点に付随する電力は顕著な差が無く、電力の差だけで選択することは難しい。
【００１９】
そこで、本発明ではこれらの偽の反射点を電力を重みとして反射点の平均位置を求めることにより統合する。つまり、
(1) 互い位置が所定の範囲内の反射点をグループ化する。
(2) グループgに属する反射点を以下の方法で統合する
【数１０】
により、統合した反射点の位置（遅延時間）<ｔg>と、それに付随する電力Egを求める。ここで、ｔgiはグループgに属する反射点の遅延時間、Ｅgiは同反射点に付随する電力である。また、Ｌgは同グループに含まれる反射点の個数である。
(3) 統合後の反射点を付随する電力が大きなものから所定の個数選択する。
前述の統合前の反射点にこのような処理を施した結果が図６（ｃ）であり、黒線で表される真の反射点とともに、処理された結果が破線で表されている。真の反射点からずれた位置に統合される場合もあるが、概ね真の反射点に近い位置に反射点が得られている。
【００２０】
［対象物の位置の推定］
統合前の反射点の位置が雑音の状態によって変化することから、統合後の反射点の位置も雑音の状況に応じて変化する。そこで、位置推定の信頼性を高めるために、全受信データの組ｓ = [ｓ0，ｓ1，… ｓNx-1]を用いて対象物の位置を推定する。
【００２１】
［真の遅延時間の計算法］
位置推定法の説明の前提として、対象物の位置が与えられたときに測定点で得られる反射信号の遅延時間を計算する方法を説明する。
点ｐTから送信された信号が位置ｑ(ｘq,ｙq) にある点状対象物Ｑによって反射され、点ｐRで受信されるときの反射信号の遅延時間T (ｐT,ｐR,ｑ) は、送信機ＸＴを発射した信号が点状対象物Ｑに達し、さらに対象物Ｑで反射して受信機ＸＲに達するのに要する時間であるから
【数１１】
で与えられる。ここで|| ｖ || はベクトルｖのノルムを表す。
【００２２】
［遅延時間のグループ化］
全受信データの組ｓ = [ｓ0，ｓ1，… ｓNx-1]に対して、［反射の遅延時間と振幅の推定］の項で述べられている方法により遅延時間の推定を行う。すると、各受信信号ｓj，j=0，…Nx-1に対して遅延時間の組
【数１２ａ】
が得られる。適当なしきい値δt、δｙを用意して
【数１２ｂ】
である遅延時間を同一グループに帰属させ、遅延時間をグループ化する。この処理によりそれぞれの要素数がＬGm個(m=0，…NG-1)のNG個のグループが作成されたとする。グループｍ (m=0，…NG-1)は推定遅延時間の組
【数１２ｃ】
によって構成される。ここで、Ｔmj(ｐMkj)はｔm = [tM0，tM1，… tMM-1]T のいずれかであり、ｐMkjは測定点の座標ＰM = [ｐM0，ｐM1，… ｐMNy-1]の何れかである。
【００２３】
［推定法］
測定位置ｐMj ，j = 1，… Ｎy-1での測定における反射点Ｑからの反射信号の真の遅延時間を遅延時間ｔq(ＰMj)とする。同一グループjに属する推定遅延時間
に対して、遅延時間の推定誤差をｅ(ｐMj) とすると
【数１３】
である。推定誤差e(ｐMj)が平均値0、分散σの正規分布をなすとすると、真の遅延時間が
であるときに推定遅延時間が
である確率は
【数１４】
であるから、
【数１５ａ】
を最大にする
を求めれば、与えられた遅延時間の推定値の下で最も事後確率が高くなる対象物の位置が求められる。つまり、
【数１５ｂ】
真の遅延時間は未知であるから、分散σも未知である。そこで、本発明では分散 の適当な推定値を用いる。例えば、遅延時間の推定値の変動の分散を用いる。さらに、遅延時間の推定値に重みをつけることも可能である。そこで、本発明では遅延時間の推定値
に付随する電力
を重みとして、
【数１６ａ】
を最大とする方法も採用する。つまり、
【数１６ｂ】
本発明では遅延時間をグループ化し、グループごとに式(15)あるいは(16)を適用する。このとき、以下のような手順で処理を行う。
【００２４】
(1) グループに含まれる遅延時間推定値
の中で最小のものを見つける。この遅延時間推定値の最小値をｔgMINとし、この最小値が測定された測定位置をｐgMIN = (ｘgMIN,ｙgMIN)とする。
(2) 測定位置ｐgMINの直下（同じｙ座標）に対象物があると仮定して、つまり、 ｙObj = ｙgMINと仮定して式(11)により対象物のｘ座標ｘObjを求める。
(3) 以下の(a)または(b)の方法で対象物の位置を推定する。
(a) (2)で得た座標(ｘObj,ｙObj)を初期値として、ニュートン法、共役勾配法、最急降下法などで(15b)または(16b)を解いて対象物の位置ｑ(ｘq,ｙq) と評価値ｆ(q)を求める。評価値が所定の値より大きい場合は、これらを対象物の推定位置ｑE(ｘq,ｙq) と評価値ｆE (q)とする。
(b) 座標(ｘObj,ｙObj)の周辺に適当な間隔で格子点を与え、各格子点で式(15a)あるいは(16a)を評価し、評価値が所定の値より大きかった場合、最大の格子点を対象物の推定位置ｑE (ｘq,ｙq)とする。また、この格子点での評価値を対象物の推定位置での評価値ｆE (ｑ)とする。
(4) 各グループを処理して得られた対象物の推定位置に対して、推定位置間の距離が所定の距離より小さなものをグループ化する。各グループに属する指定位置ｑE (ｘq,ｙq)に対して評価値ｆE (ｑ)とする重み付け平均を求め、平均位置を対象物の位置、評価値の和を対象物の反射振幅とする。
【産業上の利用可能性】
【００２５】
本発明は、地中レーダ等の近距離レーダの画像を処理する方法であり、各測定位置における受信信号をデコンボリューション等により処理して反射信号の遅延時間を推定し、推定された遅延時間と測定位置の差が所定の範囲内にあるものを同じグループに帰属させることにより推定遅延時間をグループ化する。次に、グループごとに同一グループの画素の中で最も遅延時間が短い位置の近傍に反射点があると仮定して、同グループの画素による反射点の推定を行い、推定された反射点に同グループに属する遅延時間の評価値と推定誤差によって決まる評価値を与え、近距離にある推定位置を統合することにより、同一対象物による可能性が高い推定位置を統合して、推定位置の確度を高める。本発明の画像処理方法により、管路等の人工物を地中レーダなどにより探査する場合、測定現場で目的の人工物を容易に判定できる効率的な処理が可能となる。
【図面の簡単な説明】
【００２６】
【図１】受信アンテナのみを移動して測定する場合の測定方法
【図２】図１の測定した場合の受信信号の様子
【図３】送信／受信アンテナを一体で移動して測定する場合の測定方法
【図４】送信／受信アンテナを一体で移動して測定した場合の受信信号の様子
【図５】本発明に係る実施例の処理フローを示す図
【図６】本発明に係る実施例に基づいて処理された測定信号を表す図
【技術分野】
【０００１】
管路等の人工物を地中レーダなどにより探査した場合の画像を処理して測定現場で目的の人工物を容易に判定できる効率的な処理法を提供する。
【背景技術】
【０００２】
地中レーダ等の近距離レーダの画像を処理する方法に関する。特に、人工物など周囲との境界が明確な対象からの反射信号を特定して対象の判定を容易にする方法に関する。本方法は地中レーダだけでなく、対象物の像を得るのに開口合成的な手法を必要とする音波探査および地震探査にも適用可能である。
【０００３】
【特許文献１】特開２００４−１９８１９５号公報
【特許文献２】特開２００１−０３３５６４号公報
【特許文献３】特開２００２−１０７４４９号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
地中レーダは指向性が低いので対象からの反射を分離するための処理が必要である。
現在このために使用されている手法は凡そ以下のように分類される。
・マイグレーション
開口合成的な方法である。一般に、複数の対象物の反射が重なり合っている場合には対象物の無い場所に虚像を生じる欠点がある。計算方法には遅延時間を計算する方法（ディフラクションスタックマイグレーション）とFFTを用いる方法がある。前者は格子点ごとに計算処理を行なう必要があるために多大な量の計算を必要とし、調査結果を現場で処理するような用途には適していない。一方、後者は計算処理に要する時間は少ないが、虚像を発生しやすいという欠点がある。
特に、地中レーダの場合は測定対象とする空間の大きさに比べてレーダ信号の持続時間が長いという特殊な事情により、これらの方法には対象物の明瞭な位置を得るのが難しいという欠点がある。
・トモグラフィ
計算量が極めて多大であるために、現場での処理にはまったく適していない。
以上のように、現在使われている処理法は、多大な計算時間を要する、また、特に複数の対象物が近接して存在する場合に虚像を生じるという欠点がある。本発明はこれら２つの課題を解決することを目的とする。
【課題を解決するための手段】
【０００５】
図５は本発明の処理フローの一例を示す図である。
ステップ１０１：地中レーダの送信アンテナおよび受信アンテナの少なくとも一方を移動しながら反射信号を取得する。
ステップ１０２：デコンボリューション、MUSIC法等により処理して反射信号の遅延時間を推定する。
ステップ１０３：受信信号を孤立反射体からの反射信号の重なりとみなし、ステップ１０２で得られた遅延時間を最も適切に説明可能な孤立物体の位置を推定する。
ステップ１０４：推定された遅延時間を確からしさを基に統合して、遅延時間の推定誤差を軽減するとともに、確からしさの低い遅延時間を除外して偽の対象物が発生するのを抑制する。
ステップ１０５：推定された遅延時間と測定位置の差が所定の範囲内にあるものを同じグループに帰属させることにより推定遅延時間をグループ化し、グループごとに同一グループの画素の中で最も遅延時間が短い位置の近傍に反射点があると仮定して、同グループの画素による反射点の推定を行い、推定された反射点に同グループに属する遅延時間の評価値と推定誤差によって決まる評価値を与える。
ステップ１０６：評価値の低い推定位置を除外することにより偽の対象物が発生するのを抑制し、 近距離にある推定位置を統合することにより、同一対象物による可能性が高い推定位置を統合して、推定位置の確度を高める。
以下で各処理の詳細を述べる。
【０００６】
［測定方法］
解決法を説明するために、本発明が扱うデータの取得方法について地中レーダを例に説明する。
地中レーダでは、通常受信アンテナおよび送信アンテナの少なくとも一方を測定対象領域上の地表面で移動しながら送信アンテナから電磁波が発射された時刻に同期して受信アンテナで信号を受信する。地震波探査では受波器および送波器のいずれか一方を移動して測定する方法が多く行なわれるが、地中レーダでは送信アンテナと受信アンテナを対にして同時に移動しながら測定する場合が大部分である。
【０００７】
本発明による画像処理は
(1) 送信アンテナおよび受信アンテナの少なくとも一方を移動しながら反射信号を取得する。
(2) 各測定位置における受信信号をデコンボリューション等により処理して、反射信号の遅延時間を推定する。
(3) 受信信号を孤立反射体からの反射信号の重なりとみなし、(2)で得られた遅延時間を最も適切に説明可能な孤立物体の位置を推定する。
という手順で行う。
【０００８】
手順(2)で推定された遅延時間は一般に誤差および偽の推定時間を含む。そこで本発明では、
(A) 推定された遅延時間を確からしさを基に統合して、遅延時間の推定誤差を軽減するとともに、確からしさの低い遅延時間を除外して偽の対象物が発生するのを抑制する。
また、本発明の目的は孤立物体とみなせる対象物探索のための現場における迅速な処理であるので、手順(3)においては、
(B) 手順(2)で推定された遅延時間と測定位置の差が所定の範囲内にあるものを同じグループに帰属させることにより推定遅延時間をグループ化する。グループごとに、同一グループの画素の中で最も遅延時間が短い位置の近傍に反射点があると仮定して、同グループの画素による反射点の推定を行うことにより、処理時間を短縮する。このとき、推定された反射点に同グループに属する遅延時間の評価値と推定誤差によって決まる評価値を与える。
【０００９】
(C) 適切さの低い推定位置除外することにより偽の対象物が発生するのを抑制する。
(D) 近距離にある推定位置を統合することにより、同一対象物による可能性が高い推定位置を統合して、推定位置の確度を高める。
【発明の効果】
【００１０】
(1) 処理速度が速い。
(2) 対象物の位置と反射強度を特定できるために、連続物体の同定など高度な処理を行なうことが容易になる。
(3) 対象物の推定位置が明瞭になる。そのため、測定結果の解釈が容易になる。
【実施例】
【００１１】
［受信データ］
説明を簡単にするために、また、通常の地中レーダによる測定を考慮して、送信機XT（送信アンテナ）と受信機XR（受信アンテナ）を所定の相互位置関係に保持しながら移動して受信信号を受信する。ここで、送信機XTとは測定対象領域に電磁波、地震波等のセンシングを行なう信号を送出する手段を指す。また、受信機XRとは電磁波、地震波等センシングを行なう信号を検知する手段を指す。
説明のために座標を定義する。信号が送出される方向をｘ軸、送信機XT、受信機XRが移動する方向をｙ軸、これら2軸に直交して右手系をなす方向をz軸とする。
ｘ = 0 の面に沿って送信機XT、受信機XRを移動して測定を行なう場合を例にとって説明する。送信機XT、受信機XRを移動する面が平坦である限りこの仮定により一般性を失うことは無い。説明の都合上以後この面をセンシング面と呼ぶことにする。測定点はセンシング面上に並んでいるものとする。送信機XT、受信機XRの移動は直線に沿ってでも、曲線に沿ってでも構わない。
【００１２】
以下の説明では、アンテナを直線的に移動して2次元のデータを取得するものとする。また、送信アンテナのみを移動する場合、受信アンテナのみを移動する場合もフィルタに変更を加えるだけで全く同様に扱うことができる。
測定点がｙ座標方向にＮy個並んでいる。測定点の位置ＰMをｙ座標方向に番号付けして
【数１】
とする。ここで、z座標は省略した。また、上付き添え字Tは転置を表す。
送信アンテナと受信アンテナはｙ軸方向に並んでおり、一定の間隔を保って移動するから、それぞれの基準点（送信アンテナと受信アンテナの中点）から位置をｄT、ｄRとすると、測定点の位置ベクトルの組ＰMを用いて送信アンテナの位置ベクトルの組ＰTと受信アンテナの位置ベクトルの組ＰRは
【数２】
と表すことができる。
【００１３】
受信信号の長さはMであり、等間隔tsでサンプリングされているものとする。簡単のために送信アンテナから信号が送出された時刻を0とする。すると、受信信号を受信する受信時刻の組ｔMは信号を発射した時刻を始点とする相対値であり、
【数３】
となる。
測定点ＰMjおける受信信号sjを
【数４】
と表す。全受信データの組sは
【数５】
と表すことができる。
【００１４】
［画像空間］
測定対象領域内の対象物の位置および性状の推定値を表示するために画像空間を定義する。画像空間は測定対象領域内に設けられた矩形格子点から構成される。ｙ座標軸方向のj 番目の格子点のｙ座標ｙI jは受信データを取得した測定点のｙ座標である。つまり、
【数６ａ】
ｘ座標は、信号受信時刻tMの間に信号が往復する距離である。したがって、ｘ 軸方向のm番目の格子点のｘ座標ｘImは、時刻tMmを用いると
【数６ｂ】
となる。ここで、csignalは信号の伝搬速度である。
【００１５】
［反射の遅延時間と振幅の推定］
地中レーダが送・受信する信号は、送信アンテナと受信アンテナの過渡応答により拡がりを持っているパルス状信号であるために、個々の対象物からの反射信号は、レーダ画像上で無視できない拡がりを持っている。同一の測定点で周波数を所定の間隔で変えながら測定した周波数ごとの受信信号の振幅と位相である場合も、逆フーリエ変換を施せば同様なパルス状信号になるので、状況は同じである。
この拡がりは測定対象である空間の大きさの数10％に達するので、このような画像から対象物の位置を明瞭に推定することは難しい。そこで、本発明においては、反射の遅延時間と振幅から対象物の位置を推定する。一般に逆フィルタによるデコンボリューションは雑音が発生しやすい欠点がある。送・受信アンテナの過渡応答を正確に知ることが難しいこと、反射信号のスペクトルが対象物（反射体）の深さに依存することから、地中レーダの場合は良好なデコンボリューションを行うことは特に困難である。そこで、事後確率を最大にする方法(Jensen,J. A.,Leeman,S (1994): Non-parametric Estimation of Ultrasound Pulses,IEEE Trans. Biomed. Eng. 41,No. 10,929-936.)や、いわゆるMUSIC法など(S. Shrestha and I. Arai: Signal ＰRocessing of Ground Penetrating Radar Using Spectral Estimation Techniques to Estimate the Position of Buried Targets)を用いて反射信号の遅延時間と振幅を推定する。後者は前者に比べて処理時間が少ないものの、偽の反射を発生しやすい欠点がある。
そこで、本発明では処理時間が短い後者の方法により遅延時間を推定し、推定した遅延時間に付随する振幅を基に遅延時間の取捨選択と統合して偽の反射を低減する手法を考案している。
【００１６】
［MUSIC法による遅延時間推定］
以下処理方法について説明するが、その前に上記文献に基づいてMUSIC法による遅延時間の推定について説明する。
MUSIC法による遅延時間の推定は、受信信号Ｓj= [Sj,0 ，Sj,1 ，… Sj,M-1]Tの周波数スペクトルＳj= [Sj,0 ，Sj,1 ,… Sj,M-1]TをFFTなどにより計算する。次に、例えば振幅の絶対値が所定の値より大きいなど適当な周波数帯域中のスペクトル成分Ｓj'= [Sj,ms' ，Sj,ms+1' ，… Sj,me']Tを用いて遅延時間の推定を行う。
MUSIC法では、受信信号Ｓj'が真の信号
と雑音Ｎjから構成されていると考えて、
【数７】
とモデル化する。ここで、Ａは混合行列、つまり遅延を与える行列である。またLは反射信号（反射点）の個数である。受信信号Ｓj'の共分散行列ＲSSを求めて固有値分解する。
【数８】
行列Ｄは共分散行列ＲSSの固有値の２乗を対角要素とする行列である。ここで、絶対値が小さな固有値が雑音Ｎjによるもの、絶対値が大きな固有値を真の信号Ｓjによるものと仮定し、(8)式をさらに
【数９】
と分解する。ここで、ＤSは絶対値が大きな固有値、つまり、真の信号
に対応する固有値、ＵSはこれらの固有値に対応する固有ベクトルである。ＤNは絶対値が小さな固有値、つまり、雑音Ｎjに対応する固有値、ＵNはこれらの固有値に対応する固有ベクトルである。
【００１７】
受信信号Ｓj'がなす空間を、信号空間（真の信号に対応する固有ベクトルなす空間）と雑音空間（雑音に対応する固有ベクトルがなす空間）に分割する。信号空間のベクトルは雑音空間に直交するから、混合行列Ａの列ベクトルａj=[ａj,0 ，ａj,1 ，… ａj,L'-1]Tで、雑音空間に直交するものを求めると、遅延時間を得ることができる。列ベクトルａj,l (l = 0，… L'-1)を信号空間を張る固有ベクトルで展開して、展開係数とＤSから列ベクトルａj,l (l = 0，… L'-1)に付随する電力を得ることができる。
【００１８】
［偽反射の抑制法］
次に、本発明で用いる偽の反射を抑制する方法について説明する。
MUSIC法では、雑音に対応する固有値の個数、あるいは真の信号に対応する固有値の個数を指定する必要がある。地中レーダの受信信号では、反射点の個数は未知である。適当に余裕を持った個数の反射点があるもとして処理を行う。すると、偽の反射が生ずる。予め反射点の個数が分かっている場合でも、付随する電力が大きなものから与えれた個数だけ反射を選択した場合に、検知すべき反射が全てその中に含まれるとは限らず、個数に余裕を持って遅延時間を推定する必要がある場合が多い。
図６（ａ）は測定された信号を表し、図６（ｂ）および図６（ｃ）は図６（ａ）に示される測定信号を処理して得られた信号を表している。この場合、図５（ｂ）の黒線で表される真の反射点以外に、破線で表される多くの偽の反射点が得られる。図から分かるようにこれらの偽の反射点は真の反射点の周辺に生ずることが多い。真の反射点に近い位置に多くの反射点が推定されることに変わりは無いが、実際の位置は雑音の状態等により変化する。これらの反射点に付随する電力は顕著な差が無く、電力の差だけで選択することは難しい。
【００１９】
そこで、本発明ではこれらの偽の反射点を電力を重みとして反射点の平均位置を求めることにより統合する。つまり、
(1) 互い位置が所定の範囲内の反射点をグループ化する。
(2) グループgに属する反射点を以下の方法で統合する
【数１０】
により、統合した反射点の位置（遅延時間）<ｔg>と、それに付随する電力Egを求める。ここで、ｔgiはグループgに属する反射点の遅延時間、Ｅgiは同反射点に付随する電力である。また、Ｌgは同グループに含まれる反射点の個数である。
(3) 統合後の反射点を付随する電力が大きなものから所定の個数選択する。
前述の統合前の反射点にこのような処理を施した結果が図６（ｃ）であり、黒線で表される真の反射点とともに、処理された結果が破線で表されている。真の反射点からずれた位置に統合される場合もあるが、概ね真の反射点に近い位置に反射点が得られている。
【００２０】
［対象物の位置の推定］
統合前の反射点の位置が雑音の状態によって変化することから、統合後の反射点の位置も雑音の状況に応じて変化する。そこで、位置推定の信頼性を高めるために、全受信データの組ｓ = [ｓ0，ｓ1，… ｓNx-1]を用いて対象物の位置を推定する。
【００２１】
［真の遅延時間の計算法］
位置推定法の説明の前提として、対象物の位置が与えられたときに測定点で得られる反射信号の遅延時間を計算する方法を説明する。
点ｐTから送信された信号が位置ｑ(ｘq,ｙq) にある点状対象物Ｑによって反射され、点ｐRで受信されるときの反射信号の遅延時間T (ｐT,ｐR,ｑ) は、送信機ＸＴを発射した信号が点状対象物Ｑに達し、さらに対象物Ｑで反射して受信機ＸＲに達するのに要する時間であるから
【数１１】
で与えられる。ここで|| ｖ || はベクトルｖのノルムを表す。
【００２２】
［遅延時間のグループ化］
全受信データの組ｓ = [ｓ0，ｓ1，… ｓNx-1]に対して、［反射の遅延時間と振幅の推定］の項で述べられている方法により遅延時間の推定を行う。すると、各受信信号ｓj，j=0，…Nx-1に対して遅延時間の組
【数１２ａ】
が得られる。適当なしきい値δt、δｙを用意して
【数１２ｂ】
である遅延時間を同一グループに帰属させ、遅延時間をグループ化する。この処理によりそれぞれの要素数がＬGm個(m=0，…NG-1)のNG個のグループが作成されたとする。グループｍ (m=0，…NG-1)は推定遅延時間の組
【数１２ｃ】
によって構成される。ここで、Ｔmj(ｐMkj)はｔm = [tM0，tM1，… tMM-1]T のいずれかであり、ｐMkjは測定点の座標ＰM = [ｐM0，ｐM1，… ｐMNy-1]の何れかである。
【００２３】
［推定法］
測定位置ｐMj ，j = 1，… Ｎy-1での測定における反射点Ｑからの反射信号の真の遅延時間を遅延時間ｔq(ＰMj)とする。同一グループjに属する推定遅延時間
に対して、遅延時間の推定誤差をｅ(ｐMj) とすると
【数１３】
である。推定誤差e(ｐMj)が平均値0、分散σの正規分布をなすとすると、真の遅延時間が
であるときに推定遅延時間が
である確率は
【数１４】
であるから、
【数１５ａ】
を最大にする
を求めれば、与えられた遅延時間の推定値の下で最も事後確率が高くなる対象物の位置が求められる。つまり、
【数１５ｂ】
真の遅延時間は未知であるから、分散σも未知である。そこで、本発明では分散 の適当な推定値を用いる。例えば、遅延時間の推定値の変動の分散を用いる。さらに、遅延時間の推定値に重みをつけることも可能である。そこで、本発明では遅延時間の推定値
に付随する電力
を重みとして、
【数１６ａ】
を最大とする方法も採用する。つまり、
【数１６ｂ】
本発明では遅延時間をグループ化し、グループごとに式(15)あるいは(16)を適用する。このとき、以下のような手順で処理を行う。
【００２４】
(1) グループに含まれる遅延時間推定値
の中で最小のものを見つける。この遅延時間推定値の最小値をｔgMINとし、この最小値が測定された測定位置をｐgMIN = (ｘgMIN,ｙgMIN)とする。
(2) 測定位置ｐgMINの直下（同じｙ座標）に対象物があると仮定して、つまり、 ｙObj = ｙgMINと仮定して式(11)により対象物のｘ座標ｘObjを求める。
(3) 以下の(a)または(b)の方法で対象物の位置を推定する。
(a) (2)で得た座標(ｘObj,ｙObj)を初期値として、ニュートン法、共役勾配法、最急降下法などで(15b)または(16b)を解いて対象物の位置ｑ(ｘq,ｙq) と評価値ｆ(q)を求める。評価値が所定の値より大きい場合は、これらを対象物の推定位置ｑE(ｘq,ｙq) と評価値ｆE (q)とする。
(b) 座標(ｘObj,ｙObj)の周辺に適当な間隔で格子点を与え、各格子点で式(15a)あるいは(16a)を評価し、評価値が所定の値より大きかった場合、最大の格子点を対象物の推定位置ｑE (ｘq,ｙq)とする。また、この格子点での評価値を対象物の推定位置での評価値ｆE (ｑ)とする。
(4) 各グループを処理して得られた対象物の推定位置に対して、推定位置間の距離が所定の距離より小さなものをグループ化する。各グループに属する指定位置ｑE (ｘq,ｙq)に対して評価値ｆE (ｑ)とする重み付け平均を求め、平均位置を対象物の位置、評価値の和を対象物の反射振幅とする。
【産業上の利用可能性】
【００２５】
本発明は、地中レーダ等の近距離レーダの画像を処理する方法であり、各測定位置における受信信号をデコンボリューション等により処理して反射信号の遅延時間を推定し、推定された遅延時間と測定位置の差が所定の範囲内にあるものを同じグループに帰属させることにより推定遅延時間をグループ化する。次に、グループごとに同一グループの画素の中で最も遅延時間が短い位置の近傍に反射点があると仮定して、同グループの画素による反射点の推定を行い、推定された反射点に同グループに属する遅延時間の評価値と推定誤差によって決まる評価値を与え、近距離にある推定位置を統合することにより、同一対象物による可能性が高い推定位置を統合して、推定位置の確度を高める。本発明の画像処理方法により、管路等の人工物を地中レーダなどにより探査する場合、測定現場で目的の人工物を容易に判定できる効率的な処理が可能となる。
【図面の簡単な説明】
【００２６】
【図１】受信アンテナのみを移動して測定する場合の測定方法
【図２】図１の測定した場合の受信信号の様子
【図３】送信／受信アンテナを一体で移動して測定する場合の測定方法
【図４】送信／受信アンテナを一体で移動して測定した場合の受信信号の様子
【図５】本発明に係る実施例の処理フローを示す図
【図６】本発明に係る実施例に基づいて処理された測定信号を表す図
【特許請求の範囲】
【請求項１】
レーダ信号の送信手段および受信手段の少なくとも一方を各測定位置ごとに送信信号を送出して移動し対象からの反射信号データから対象物の位置を推定する地中レーダ画像処理方法であって、
各測定位置における受信信号から反射信号の遅延時間を推定する第１のステップと、
受信信号を孤立反射体からの反射信号の重なりとみなし第１のステップで得られた遅延時間から孤立物体の位置を推定する第２のステップと、
第２のステップで推定された遅延時間を確からしさに基づいて統合し、確からしさの低い遅延時間を除外する第３のステップと、
第３のステップにより統合された推定遅延時間が所定の範囲内にあるものを同じグループにグループ化し、グループごとに同一グループの画素の中で最も遅延時間が短い位置の近傍に反射点があると仮定して、同グループの画素による反射点の推定を行い、推定された反射点に同グループに属する遅延時間の評価値と推定誤差によって決まる評価値を与える第４のステップと、
第４のステップにより与えられた評価値の低い推定位置を除外して偽の対象物を抑制しする第５のステップと、
近距離にある推定位置を統合することにより、同一対象物による可能性が高い推定位置を統合して、推定位置の確度を高める第６のステップからなる
地中レーダ画像処理方法。
【請求項２】
請求項１記載の地中レーダ画像処理法であって、
第４のステップにおいて対象物の位置を推定するための評価関数として、仮定した対象物の位置に対する反射信号の遅延時間と手順1で得られた遅延時間に対する評価関数
あるいは
を用いることを特徴とする地中レーダ画像処理法。
【請求項３】
請求項２記載の地中レーダ画像処理法であって、
評価関数
あるいは
において第３のステップで推定された遅延時間の分散を確率密度の分散として用いることを特徴とする地中レーダ画像処理法。
【請求項１】
レーダ信号の送信手段および受信手段の少なくとも一方を各測定位置ごとに送信信号を送出して移動し対象からの反射信号データから対象物の位置を推定する地中レーダ画像処理方法であって、
各測定位置における受信信号から反射信号の遅延時間を推定する第１のステップと、
受信信号を孤立反射体からの反射信号の重なりとみなし第１のステップで得られた遅延時間から孤立物体の位置を推定する第２のステップと、
第２のステップで推定された遅延時間を確からしさに基づいて統合し、確からしさの低い遅延時間を除外する第３のステップと、
第３のステップにより統合された推定遅延時間が所定の範囲内にあるものを同じグループにグループ化し、グループごとに同一グループの画素の中で最も遅延時間が短い位置の近傍に反射点があると仮定して、同グループの画素による反射点の推定を行い、推定された反射点に同グループに属する遅延時間の評価値と推定誤差によって決まる評価値を与える第４のステップと、
第４のステップにより与えられた評価値の低い推定位置を除外して偽の対象物を抑制しする第５のステップと、
近距離にある推定位置を統合することにより、同一対象物による可能性が高い推定位置を統合して、推定位置の確度を高める第６のステップからなる
地中レーダ画像処理方法。
【請求項２】
請求項１記載の地中レーダ画像処理法であって、
第４のステップにおいて対象物の位置を推定するための評価関数として、仮定した対象物の位置に対する反射信号の遅延時間と手順1で得られた遅延時間に対する評価関数
あるいは
を用いることを特徴とする地中レーダ画像処理法。
【請求項３】
請求項２記載の地中レーダ画像処理法であって、
評価関数
あるいは
において第３のステップで推定された遅延時間の分散を確率密度の分散として用いることを特徴とする地中レーダ画像処理法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図２】
【図３】
【図４】
【図５】
【図６】
【公開番号】特開２００７−１６３２７１（Ｐ２００７−１６３２７１Ａ）
【公開日】平成１９年６月２８日（２００７．６．２８）
【国際特許分類】
【出願番号】特願２００５−３５９２４１（Ｐ２００５−３５９２４１）
【出願日】平成１７年１２月１３日（２００５．１２．１３）
【出願人】（０００２０８８９１）ＫＤＤＩ株式会社 (2,700)
【Ｆターム（参考）】
【公開日】平成１９年６月２８日（２００７．６．２８）
【国際特許分類】
【出願日】平成１７年１２月１３日（２００５．１２．１３）
【出願人】（０００２０８８９１）ＫＤＤＩ株式会社 (2,700)
【Ｆターム（参考）】
[ Back to top ]