WO2024034000A1

WO2024034000A1 - データ処理方法、計測システム、及び、プログラム

Info

Publication number: WO2024034000A1
Application number: PCT/JP2022/030378
Authority: WO
Inventors: 康成森
Original assignee: 株式会社三井E&S
Priority date: 2022-08-09
Filing date: 2022-08-09
Publication date: 2024-02-15

Abstract

物体に放射した波動の散乱波を解析するデータ処理方法であって、ｘｙ平面上に２次元に配列された複数の送信点p₁(x'₁, y'₁, z'₁)から、物体に波動を放射し、物体上の反射点(x, y, z)において反射率f(x, y, z)で反射した散乱波を、ｘｙ平面上に２次元に配列された複数の受信点p₂(x'₂, y'₂, z'₂)で計測値s(x'₁, x'₂, y'₁, y'₂, z'₁, z'₂, k)として受信し、計測値s(x'₁, x'₂, y'₁, y'₂, z'₁, z'₂, k)を４重フーリエ変換してS(k'_x1, k'_x2, k'_y1, k'_y2, z'₁, z'₂, k)を求め、ｘ方向の計測間隔で定まるx'₁方向のナイキスト波数をk'_x1,nyqとし、x'₂方向のナイキスト波数をk'_x2,nyqとすると、-k_x1≦k'_x1≦k_x1、かつ、-k_x2≦k'_x2≦k_x2（但し、k'_x1,nyq≦k_x1、かつ、k'_x2,nyq≦k_x2）の範囲でｘ方向の変数置換処理を行い、ｙ方向の計測間隔で定まるy'₁方向のナイキスト波数をk'_y1,nyqとし、y'₂方向のナイキスト波数をk'_y2,nyqとすると、-k_y1≦k'_y1≦k_y1、かつ、-k_y2≦k'_y2≦k_y2（但し、k'_y1,nyq≦k_y1、かつ、k'_y2,nyq≦k_y2）の範囲でｙ方向の変数置換処理を行い、３重逆フーリエ変換して、反射率f(x, y, z)を求めるデータ処理方法。

Description

データ処理方法、計測システム、及び、プログラム

　本発明は、空間に生成する電磁波等の波動の周波数と空間の空間座標とによって値が定まる波動の計測データを、コンピュータを用いて処理するデータ処理方法、計測システム、及び、プログラムに関する。

　従来、コンクリートや木材等の非金属の構造物の内部を非破壊で検査するレーダ装置が知られている。従来のレーダ装置は、平面上に複数のアンテナが配置されたアレイアンテナを有する。アレイアンテナは、例えば、平面アンテナ等のアンテナが一方向に並んだ構成を有し、送信用アレイアンテナと受信用アレイアンテナが近接して配置される。また、レーダ装置は、構造物の内部を精度よく計測するために、電磁波の周波数を設定された周波数間隔で変更しながら、広帯域の周波数で測定対象物を計測する。

　アレイアンテナを有するレーダ装置に関して、例えば、複数の平面アンテナで構成された送信用アレイアンテナと受信用アレイアンテナが共通の誘電体基板に形成されたレーダ装置が知られている（特開２０１５－０９５８４０号公報、以下、「特許文献１」）。従来のレーダ装置では、送信用アレイアンテナの平面アンテナの配列方向は、受信用アレイアンテナの平面アンテナの配列方向と平行である。受信用アレイアンテナの平面アンテナの配列方向における位置は、送信用アレイアンテナの隣接する平面アンテナの２つの位置の中間にある。

　また、逆問題の解析を汎用的かつ高速に行い、物体内部の情報を簡便に映像化することができる散乱トモグラフィ方法が知られている（特許第６５５７７４７号公報、以下、「特許文献２」）。

　計測したデータから構造物の内部を映像化するために、合成開口処理が利用される。合成開口処理には、大きく、ディフラクションスタッキング法などの足し込み法と、Ｆ－Ｋマイグレーション法などのフーリエ変換を利用するものがある。
　実用的な演算時間を実現するためには、フーリエ変換を利用した合成開口処理が現実的である。

　レーダ装置では、構造物の内部を正確に検査するために、計測における空間分解能が高いことが望まれる。一般に、電磁波等の周波数を有する波動の放射によって得られるデータの空間分解能は、計測対象の構造物と送信用アレイアンテナ及び受信用アレイアンテナの計測面との間の距離が相対的に近接し、かつ、計測データの計測間隔が小さい場合、波動の中心周波数の波長によって定まる。ここで、計測対象の構造物と送信用アレイアンテナ及び受信用アレイアンテナの計測面との間の距離は、例えば、アレイアンテナの配列長さの４分の１以下である。また、空間分解能は、アレイアンテナの各アンテナが配置される平面内の分解能である。

　例えば、アレイアンテナの各アンテナが配置される平面に沿った計測データの計測間隔が十分に小さい場合の理論空間分解能は、電磁波の往復経路を考慮して、放射する波動の周波数が周波数帯域を持って掃引される場合、周波数帯域の中心周波数における波動の波長の４分の１になる。しかし、計測データの計測間隔が大きくなり、電磁波の最小波長の４分の１を超える場合、実際の計測における空間分解能は、理想空間分解能より大きくなる。場合によっては、実際の計測における空間分解能は、計測間隔になる。

　従来のレーダ装置では、低い周波数から高い周波数まで広い周波数帯域で電磁波を計測するため、電磁波の最大波長は長くなる。このため、アレイアンテナを構成する各アンテナは大きくなり、アレイアンテナにおけるアンテナの配列方向の長さは長くなる。その結果、送受信用アレイアンテナにおけるアンテナの配置間隔は長くなり、計測データの計測間隔は、放射される電磁波の最小波長の４分の１を超え易く、空間分解能は理論分解能より低下し、しかも、計測データにはエイリアシング成分が生じ易い。空間分解能を、理想とする理論空間分解能、すなわち、中心周波数における電磁波の波長の４分の１にするには、送受信用アレイアンテナ内でアンテナの配置数を増やして、配置間隔を短くしなければならない。しかし、前述の通り、広帯域のアンテナの大きさを小さくすることは困難であるため、配置間隔を短くすることも困難である。

　本発明は、アンテナの配置数を一定に維持したまま、計測における空間分解能を向上させることができるデータ処理方法、計測システム、及び、プログラムを提供することを目的とする。

　本発明の第１の観点は、
　物体に放射した波動の散乱波を解析するデータ処理方法であって、
　ｙ軸上に配列された複数の送信点p₁(x’₁, y’₁, z’₁)から、前記物体に前記波動を放射し、
　前記物体上の反射点(x, y, z)において反射率f(x, y, z)で反射した前記散乱波を、ｙ軸上に配列された複数の受信点p₂(x’₂, y’₂, z’₂)で計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)として受信し、
　前記計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)を式（１）より３重フーリエ変換してS(k’_x1, k’_x2, k’_y1, k’_y2, z’₁, z’₂, k)を求め、

　ｘ方向の計測間隔で定まるｘ方向のナイキスト波数をk_x,nyqとすると、-k_x1≦k_x≦k_x1（但し、k_x,nyq≦k_x1）の範囲でｘ方向の変数置換処理を行い、
　ｙ方向の計測間隔で定まるy’₁方向のナイキスト波数をk’_y1,nyqとし、y’₂方向のナイキスト波数をk’_y2,nyqとすると、-k_y1≦k’_y1≦k_y1、かつ、-k_y2≦k’_y2≦k_y2（但し、k’_y1,nyq≦k_y1、かつ、k’_y2,nyq≦k_y2）の範囲でｙ方向の変数置換処理を行い、
　式（２）より３重逆フーリエ変換して、前記反射率f(x, y, z)を求める、

　データ処理方法である。
但し、
x’ = x’₁ = x’₂
z’₁= z’₂ = 0
kは、伝播する前記波動の波数、
k’_x1, k’_y1, k’_z1は、前記送信点p₁(x’₁, y’₁, z’₁)から前記反射点(x, y, z)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k’_x2, k’_y2, k’_z2は、前記反射点(x, y, z)から前記受信点p₂(x’₂, y’₂, z’₂)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k_x = k’_x1 + k’_x2, u = k’_x1 - k’_x2, k_y = k’_y1 + k’_y2, v = k’_y1 - k’_y2
である。

　本発明の第２の観点は、
　物体に放射した波動の散乱波を解析する計測システムであって、
　送受信部であって、
　　ｙ軸上に配列された複数の送信点p₁(x’₁, y’₁, z’₁)から、前記物体に前記波動を放射する送信部と、
　　前記物体上の反射点(x, y, z)において反射率f(x, y, z)で反射した前記散乱波を、ｙ軸上に配列された複数の受信点p₂(x’₂, y’₂, z’₂)で計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)として受信する受信部と、
　を有する送受信部と、
　処理装置であって、
　　前記計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)を式（１）より３重フーリエ変換してS(k’_x1, k’_x2, k’_y1, k’_y2, z’₁, z’₂, k)を求める手順と、

　　ｘ方向の計測間隔で定まるｘ方向のナイキスト波数をk_x,nyqとすると、-k_x1≦k_x≦k_x1（但し、k_x,nyq≦k_x1）の範囲でｘ方向の変数置換処理を行う手順と、
　　ｙ方向の計測間隔で定まるy’₁方向のナイキスト波数をk’_y1,nyqとし、y’₂方向のナイキスト波数をk’_y2,nyqとすると、-k_y1≦k’_y1≦k_y1、かつ、-k_y2≦k’_y2≦k_y2（但し、k’_y1,nyq≦k_y1、かつ、k’_y2,nyq≦k_y2）の範囲でｙ方向の変数置換処理を行う手順と、
　式（２）より３重逆フーリエ変換して、前記反射率f(x, y, z)を求める手順と、

　を実行する処理装置と、
　を有する、計測システムである。
但し、
x’ = x’₁ = x’₂
z’₁= z’₂ = 0
kは、伝播する前記波動の波数、
k’_x1, k’_y1, k’_z1は、前記送信点p₁(x’₁, y’₁, z’₁)から前記反射点(x, y, z)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k’_x2, k’_y2, k’_z2は、前記反射点(x, y, z)から前記受信点p₂(x’₂, y’₂, z’₂)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k_x = k’_x1 + k’_x2, u = k’_x1 - k’_x2, k_y = k’_y1 + k’_y2, v = k’_y1 - k’_y2
である。

　本発明の第３の観点は、
　物体に放射した波動の散乱波を解析するプログラムであって、
　計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)を式（１）より３重フーリエ変換してS(k’_x1, k’_x2, k’_y1, k’_y2, z’₁, z’₂, k)を求める手順と、

　ｘ方向の計測間隔で定まるｘ方向のナイキスト波数をk_x,nyqとすると、-k_x1≦k_x≦k_x1（但し、k_x,nyq≦k_x1）の範囲でｘ方向の変数置換処理を行う手順と、
　ｙ方向の計測間隔で定まるy’₁方向のナイキスト波数をk’_y1,nyqとし、y’₂方向のナイキスト波数をk’_y2,nyqとすると、-k_y1≦k’_y1≦k_y1、かつ、-k_y2≦k’_y2≦k_y2（但し、k’_y1,nyq≦k_y1、かつ、k’_y2,nyq≦k_y2）の範囲でｙ方向の変数置換処理を行う手順と、
　式（２）より３重逆フーリエ変換して、反射率f(x, y, z)を求める手順と、

　をコンピュータに実行させるプログラムである。
但し、
x’ = x’₁ = x’₂
z’₁= z’₂ = 0
kは、伝播する前記波動の波数、
k’_x1, k’_y1, k’_z1は、前記送信点p₁(x’₁, y’₁, z’₁)から前記反射点(x, y, z)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k’_x2, k’_y2, k’_z2は、前記反射点(x, y, z)から前記受信点p₂(x’₂, y’₂, z’₂)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k_x = k’_x1 + k’_x2, u = k’_x1 - k’_x2, k_y = k’_y1 + k’_y2, v = k’_y1 - k’_y2
である。

　本発明の第４の観点は、
　物体に放射した波動の散乱波を解析するデータ処理方法であって、
　ｘｙ平面上に２次元に配列された複数の送信点p₁(x’₁, y’₁, z’₁)から、前記物体に前記波動を放射し、
　前記物体上の反射点(x, y, z)において反射率f(x, y, z)で反射した前記散乱波を、ｘｙ平面上に２次元に配列された複数の受信点p₂(x’₂, y’₂, z’₂)で計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)として受信し、
　前記計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)を式（１）より４重フーリエ変換してS(k’_x1, k’_x2, k’_y1, k’_y2, z’₁, z’₂, k)を求め、

　ｘ方向の計測間隔で定まるx’₁方向のナイキスト波数をk’_x1,nyqとし、x’₂方向のナイキスト波数をk’_x2,nyqとすると、-k_x1≦k’_x1≦k_x1、かつ、-k_x2≦k’_x2≦k_x2（但し、k’_x1,nyq≦k_x1、かつ、k’_x2,nyq≦k_x2）の範囲でｘ方向の変数置換処理を行い、
　ｙ方向の計測間隔で定まるy’₁方向のナイキスト波数をk’_y1,nyqとし、y’₂方向のナイキスト波数をk’_y2,nyqとすると、-k_y1≦k’_y1≦k_y1、かつ、-k_y2≦k’_y2≦k_y2（但し、k’_y1,nyq≦k_y1、かつ、k’_y2,nyq≦k_y2）の範囲でｙ方向の変数置換処理を行い、
　式（２）より３重逆フーリエ変換して、前記反射率f(x, y, z)を求める、

　データ処理方法である。
但し、
z’₁= z’₂ = 0
kは、伝播する前記波動の波数、
k’_x1, k’_y1, k’_z1は、前記送信点p₁(x’₁, y’₁, z’₁)から前記反射点(x, y, z)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k’_x2, k’_y2, k’_z2は、前記反射点(x, y, z)から前記受信点p₂(x’₂, y’₂, z’₂)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k_x = k’_x1 + k’_x2, u = k’_x1 - k’_x2, k_y = k’_y1 + k’_y2, v = k’_y1 - k’_y2
である。

　本発明の第５の観点は、
　物体に放射した波動の散乱波を解析する計測システムであって、
　送受信部であって、
　　ｘｙ平面上に２次元に配列された複数の送信点p₁(x’₁, y’₁, z’₁)から、前記物体に前記波動を放射する送信部と、
　　前記物体上の反射点(x, y, z)において反射率f(x, y, z)で反射した前記散乱波を、ｘｙ平面上に２次元に配列された複数の受信点p₂(x’₂, y’₂, z’₂)で計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)として受信する受信部と、
　を有する送受信部と、
　処理装置であって、
　　前記計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)を式（１）より４重フーリエ変換してS(k’_x1, k’_x2, k’_y1, k’_y2, z’₁, z’₂, k)を求める手順と、

　ｘ方向の計測間隔で定まるx’₁方向のナイキスト波数をk’_x1,nyqとし、x’₂方向のナイキスト波数をk’_x2,nyqとすると、-k_x1≦k’_x1≦k_x1、かつ、-k_x2≦k’_x2≦k_x2（但し、k’_x1,nyq≦k_x1、かつ、k’_x2,nyq≦k_x2）の範囲でｘ方向の変数置換処理を行う手順と、
　ｙ方向の計測間隔で定まるy’₁方向のナイキスト波数をk’_y1,nyqとし、y’₂方向のナイキスト波数をk’_y2,nyqとすると、-k_y1≦k’_y1≦k_y1、かつ、-k_y2≦k’_y2≦k_y2（但し、k’_y1,nyq≦k_y1、かつ、k’_y2,nyq≦k_y2）の範囲でｙ方向の変数置換処理を行う手順と、
　式（２）より３重逆フーリエ変換して、前記反射率f(x, y, z)を求める手順と、

　を実行する処理装置と、
　を有する、計測システムである。
但し、
z’₁= z’₂ = 0
kは、伝播する前記波動の波数、
k’_x1, k’_y1, k’_z1は、前記送信点p₁(x’₁, y’₁, z’₁)から前記反射点(x, y, z)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k’_x2, k’_y2, k’_z2は、前記反射点(x, y, z)から前記受信点p₂(x’₂, y’₂, z’₂)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k_x = k’_x1 + k’_x2, u = k’_x1 - k’_x2, k_y = k’_y1 + k’_y2, v = k’_y1 - k’_y2
である。

　本発明の第６の観点は、
　物体に放射した波動の散乱波を解析するプログラムであって、
　計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)を式（１）より４重フーリエ変換してS(k’_x1, k’_x2, k’_y1, k’_y2, z’₁, z’₂, k)を求める手順と、

　ｘ方向の計測間隔で定まるx’₁方向のナイキスト波数をk’_x1,nyqとし、x’₂方向のナイキスト波数をk’_x2,nyqとすると、-k_x1≦k’_x1≦k_x1、かつ、-k_x2≦k’_x2≦k_x2（但し、k’_x1,nyq≦k_x1、かつ、k’_x2,nyq≦k_x2）の範囲でｘ方向の変数置換処理を行う手順と、
　ｙ方向の計測間隔で定まるy’₁方向のナイキスト波数をk’_y1,nyqとし、y’₂方向のナイキスト波数をk’_y2,nyqとすると、-k_y1≦k’_y1≦k_y1、かつ、-k_y2≦k’_y2≦k_y2（但し、k’_y1,nyq≦k_y1、かつ、k’_y2,nyq≦k_y2）の範囲でｙ方向の変数置換処理を行う手順と、
　式（２）より３重逆フーリエ変換して、反射率f(x, y, z)を求める手順と、

　をコンピュータに実行させるプログラムである。
但し、
z’₁= z’₂ = 0
kは、伝播する前記波動の波数、
k’_x1, k’_y1, k’_z1は、送信点p₁(x’₁, y’₁, z’₁)から反射点(x, y, z)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k’_x2, k’_y2, k’_z2は、前記反射点(x, y, z)から受信点p₂(x’₂, y’₂, z’₂)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k_x = k’_x1 + k’_x2, u = k’_x1 - k’_x2, k_y = k’_y1 + k’_y2, v = k’_y1 - k’_y2
である。

　本発明のデータ処理方法、計測システム、及び、プログラムによれば、アンテナの配置数を一定に維持したまま、計測における空間分解能を向上させることができる。

第１実施形態のレーダ装置の構成を示す図図１に示すアレイアンテナの構成を示す図第１実施形態のアレイアンテナと測定対象物との位置関係を説明する図第１実施形態のデータ処理方法を示すフローチャートエイリアシングの概念図第１実施形態の計測波形の一例第１実施形態におけるｘ軸に関するフーリエ変換第１実施形態におけるｙ軸に関するフーリエ変換第２実施形態のアレイアンテナと測定対象物との位置関係を説明する図第２実施形態のデータ処理方法を示すフローチャート第２実施形態におけるｘ’_１軸、ｘ’_２軸に関するフーリエ変換第２実施形態におけるｙ’_１軸、ｙ’_２軸に関するフーリエ変換第３実施形態のアレイアンテナと測定対象物との位置関係を説明する図第３実施形態のデータ処理方法を示すフローチャート第３実施形態におけるｘ軸に関するフーリエ変換第３実施形態におけるｙ’_１軸、ｙ’_２軸に関するフーリエ変換比較例１のデータ処理方法でシミュレーションした点ターゲット１のｘｚ平面画像実施例１のデータ処理方法でシミュレーションした点ターゲット１のｘｚ平面画像比較例１のデータ処理方法でシミュレーションした点ターゲット１のｙｚ平面画像実施例１のデータ処理方法でシミュレーションした点ターゲット１のｙｚ平面画像比較例１と実施例１の点ターゲット１のｘ軸方向の波形比較例１と実施例１の点ターゲット１のｙ軸方向の波形比較例１と実施例１の点ターゲット１のｚ軸方向の波形比較例１のデータ処理方法でシミュレーションした点ターゲット２のｘｚ平面画像実施例１のデータ処理方法でシミュレーションした点ターゲット２のｘｚ平面画像比較例１のデータ処理方法でシミュレーションした点ターゲット２のｙｚ平面画像実施例１のデータ処理方法でシミュレーションした点ターゲット２のｙｚ平面画像比較例１と実施例１の点ターゲット２のｘ軸方向の波形比較例１と実施例１の点ターゲット２のｙ軸方向の波形比較例１と実施例１の点ターゲット２のｚ軸方向の波形十字型ｋ’_y1－ｋ’_y2フィルタの概念図４つの点ターゲットのシミュレーション結果点ターゲットの回復率送信アンテナと受信アンテナのレイアウト超分解能処理を行わない、試験体の３次元画像超分解能処理を行った、試験体の３次元画像十字型ｋ’_y1－ｋ’_y2フィルタを適用して超分解能処理を行った、試験体の３次元画像比較例２と実施例２の点ターゲット１のｘ軸方向の波形比較例２と実施例２の点ターゲット１のｙ軸方向の波形比較例２と実施例２の点ターゲット１のｚ軸方向の波形比較例２と実施例２の点ターゲット２のｘ軸方向の波形比較例２と実施例２の点ターゲット２のｙ軸方向の波形比較例２と実施例２の点ターゲット２のｚ軸方向の波形十字型ｋ’_ｘ1－ｋ’_ｘ2フィルタの概念図

＜第１実施形態＞
　以下、第１実施形態のデータ処理方法、計測システム、及び、プログラムについて、詳細に説明する。図１は、本実施形態のレーダ装置の構成を示す。図２は、図１に示すアレイアンテナの構成を示す。図３は、本実施形態のアレイアンテナと測定対象物との位置関係を説明する図である。本実施形態では、電磁波を空間に放射する波動として説明するが、電磁波の代わりにＸ線や超音波等の空間中に伝播する波動を用いてもよい。

　本実施形態の計測システム１は、送受信部と、処理装置と、を有する。処理装置は、送受信部と一体に設けられてもよいし、送受信部とネットワークで接続された別の場所に設けられてもよい。以下の実施形態では、処理装置が送受信部と一体に設けられる例を説明する。

　図１に示す本実施形態のレーダ装置６０は、送信用アレイアンテナ及び受信用アレイアンテナ（送受信部）を用いて、電磁波の周波数を掃引しながら、電磁波を送信アンテナから放射する。そして、レーダ装置６０は、測定対象物の反射波を受信アンテナで受信して、計測データｓ（ｘ’，ｙ’，ｚ’，ｋ）を得る。計測データｓ（ｘ’，ｙ’，ｚ’，ｋ）は、ｘ座標成分、ｙ座標成分、及びｚ座標成分と電磁波の周波数とを変数とするデータである。

　レーダ装置６０は、計測ユニット６１と、データ処理ユニット（処理装置）６６と、画像表示ユニット６８とを有する。計測ユニット６１は、送信用アレイアンテナ５０と、受信用アレイアンテナ５２と、高周波スイッチ５８，５９と、高周波回路６２と、システム制御回路６４とを有する。レーダ装置６０は、１０ＭＨｚ以上、例えば１０～２０ＧＨｚの電磁波を放射するが、電磁波の周波数は、特に制限されない。

　図２に示されるように、送信用アレイアンテナ５０は、一方向に配列された複数の送信アンテナ１０ａを有する。各送信アンテナ１０ａは、測定対象物に向けて電磁波を放射する。受信用アレイアンテナ５２は、送信アンテナ１０ａの配列方向に沿って配列された複数の受信アンテナ１０ｂを有する。各受信アンテナ１０ｂは、測定対象物から反射した電磁波を受信する。
　送信用アレイアンテナ５０の送信アンテナ１０ａと、受信用アレイアンテナ５２の受信アンテナ１０ｂは、一平面上に配置される。この平面に測定対象物が対向するように、送信用アレイアンテナ５０と受信用アレイアンテナ５２が配置される。

　データ処理ユニット６６は、複数の送信アンテナ１０ａによる測定対象物に向けた送信と、複数の受信アンテナ１０ｂによる受信とによって得られる複数の計測データを処理し、測定対象物に関する画像データを算出する。本実施形態の送信アンテナ１０ａ及び受信アンテナ１０ｂは、基板に平面的にアンテナパターンが形成された平面アンテナであるが、平面アンテナに制限されない。

　送信用アレイアンテナ５０と受信用アレイアンテナ５２は、測定対象物の面に平行に移動する。すなわち、送信用アレイアンテナ５０と受信用アレイアンテナ５２は、測定対象物の表面に沿って走査しながら計測する。送信用アレイアンテナ５０と受信用アレイアンテナ５２が移動するとき、システム制御回路６４は、高周波回路６２の動作を制御する。具体的には、送信用アレイアンテナ５０と受信用アレイアンテナ５２の移動距離の単位長さ毎に、送信アンテナ１０ａを高周波スイッチ５８により切り替えつつ、電磁波を放射するように、システム制御回路６４は、高周波回路６２の動作を制御する。

　レーダ装置６０は、エンコーダ６９を有する。エンコーダ６９は、一定の移動距離ごとにパルス信号を発生する。エンコーダ６９は、送信用アレイアンテナ５０及び受信用アレイアンテナ５２の移動を感知する。
　このとき、個々の送信アンテナ１０ａから電磁波の放射が行われる度に、高周波スイッチ５９は、複数の受信アンテナ１０ｂを順次切り替えて、各受信アンテナ１０ｂに受信させる。

　なお、送信用アレイアンテナ５０から放射される電磁波の周波数を、一定の時間に、例えば１０～２０ＧＨｚの範囲で、設定された周波数間隔で掃引して、電磁波が放射される。したがって、高周波回路６２から得られる計測データは送信アンテナ１０ａの送信した位置と、受信アンテナ１０ｂの受信した位置と、周波数と、ターゲットの位置とによって値が定まるデータである。
　このとき、送信アンテナ１０ａから放射された電磁波が測定対象物で反射したときの電磁波の反射波を、電磁波を放射した送信アンテナ１０aに最も近い受信アンテナ１０ｂで受信するように、高周波スイッチ５９の動作が制御される。受信用マイクロ波増幅器（ＲＦアンプ）は、送信する送信アンテナ１０ａと受信する受信アンテナ１０ｂの対毎にゲインを変化させるように設定される場合がある。このとき、高周波回路６２は、送信アンテナ１０ａと受信アンテナ１０ｂの対の選択に応じてゲインを切り替える可変ゲイン増幅機能を有する。これにより、測定対象物中の欠陥等の検査可能な深度を大きくできる。

　本実施形態では、送信アンテナ１０ａと受信アンテナ１０ｂの配列方向は平行であり、図２に示すように、配列方向をｙ方向とする。一方、送信用アレイアンテナ５０と受信用アレイアンテナ５２の移動方向（走査方向）を、ｘ方向とする。送信用アレイアンテナ５０と受信用アレイアンテナ５２からみて、測定対象物のある方向（電磁波の送信方向）をｚ方向とする。

　なお、送信用アレイアンテナ５０と受信用アレイアンテナ５２の移動方向（走査方向）をｙ方向としてもよい。すなわち、送信アンテナ１０ａと受信アンテナ１０ｂの配列方向と同じ方向に、移動（走査）してもよい。
　また、送信用アレイアンテナ５０が１つの送信アンテナ１０ａのみを有し、受信用アレイアンテナ５２が複数の受信アンテナ１０ｂを有してもよい。この場合も、送信用アレイアンテナ５０と受信用アレイアンテナ５２の移動方向（走査方向）をｙ方向としてもよい。すなわち、受信アンテナ１０ｂの配列方向と同じ方向に、移動（走査）してもよい。

　データ処理ユニット６６は、送信用アレイアンテナ５０及び受信用アレイアンテナ５２による電磁波の送受信によって得られる計測データｓ（ｘ’，ｙ’，ｚ’，ｋ）を処理して、測定対象物の内部を表す画像データを作成する。データ処理ユニット６６は、例えばコンピュータにより構成され、記憶部６６ａに記憶されているプログラムを呼び出して起動する。これにより、データ処理ユニット６６の機能を発揮できる。すなわち、データ処理ユニット６６は、ソフトウェアモジュールで構成される。画像表示ユニット６８は、作成された画像データを用いて、測定対象物の内部の画像を表示する。

　図２は、送信用アレイアンテナ５０と受信用アレイアンテナ５２を模式的に示す。送信アンテナ１０ａと受信アンテナ１０ｂは、ｘ方向の位置がΔＬだけずれているが、以降の説明では、送信アンテナ１０ａと受信アンテナ１０ｂのｘ方向の位置は、送信アンテナ１０ａと受信アンテナ１０ｂの間の中間の丸印の点にあるものする。この丸印の点を、送受信点と呼ぶ。
　なお、送信アンテナ１０ａと受信アンテナ１０ｂのｙ方向のずれが無い場合もある。すなわち、Δｙ＝０となる場合もある。また、送信アンテナ１０ａと受信アンテナ１０ｂが共有される場合もある。すなわち、Δｙ＝０、ΔＬ＝０となる場合もある。

　したがって、測定対象物と送信用アレイアンテナ５０と受信用アレイアンテナ５２との位置関係は、図３に示すように表すことができる。
　ここで、送受信点の座標をｐ（ｘ’，ｙ’，ｚ’）とする。測定対象物の反射点（ｘ，ｙ，ｚ）における反射率をｆ（ｘ，ｙ，ｚ）とする。送受信点ｐ（ｘ’，ｙ’，ｚ’）における計測データをｓ（ｘ’，ｙ’，ｚ’，ｋ）とする。真空中の電磁波の伝播波長をλ_０とする。媒質の比誘電率をε_ｒとする。伝播する電磁波の波数をｋとする。

　このとき、送受信点ｐ（ｘ’，ｙ’，ｚ’）における計測データｓ（ｘ’，ｙ’，ｚ’，ｋ）は、以下の式で表せる。

　但し、

である。

　式（１－１）では、電磁波を球面波で表しており、距離減衰は省略されている。この距離減衰は、以降の処理を行う上で影響が小さいため、省略されている。式（１－１）中の二段目の式の被積分関数の指数部をフーリエ変換の表記で表すと、以下の式となる。これは、式（１－１）の往復球面波を３次元の平面波に分解することに等しい。

　ここで、（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）は、送受信点ｐ（ｘ’，ｙ’，ｚ’）と反射点（ｘ，ｙ，ｚ）の間で伝播する波動の往復球面波の波数ベクトルの成分である。但し、

　を満たす。

　以下、式（１－３）に基づいて、計測データｓ（ｘ’，ｙ’，ｚ’，ｋ）から反射率ｆ（ｘ，ｙ，ｚ）を導出する。まず、式（１－３）を以下のように整理する。

　ここで、｛　｝の内側の積分は、（ｘ，ｙ，ｚ）に関する３重フーリエ変換である。また、［　］の内側の積分は、（ｋ_ｘ，ｋ_ｙ）に関する２重逆フーリエ変換である。そこで、式（１－５）の両辺を（ｘ’，ｙ’）に関して２重フーリエ変換を行う。関数ｆ（ｘ，ｙ，ｚ）の３重フーリエ変換後の関数をＦ（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）とする。計測データｓ（ｘ’，ｙ’，ｚ’，ｋ）の２重フーリエ変換後の関数をＳ（ｋ_ｘ，ｋ_ｙ，ｚ’，ｋ）とする。このとき、式（１－５）は、以下の式で表される。

　式（１－６）の２行目の式の両辺を（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）について３重逆フーリエ変換すると、反射率ｆ（ｘ，ｙ，ｚ）が以下のように得られる。

　ここで、本実施形態では、図３に示すように、送受信点ｐ（ｘ’，ｙ’，ｚ’）をｘｙ平面に配置し、ｚ’＝０となるため、式（１－７）は以下のように表せる。

　以上のように、データ処理ユニット６６は、計測データｓ（ｘ’，ｙ’，ｚ’，ｋ）に基づいて、反射率ｆ（ｘ，ｙ，ｚ）を求める。

　以下、図４を参照して、本実施形態のデータ処理方法、及び、プログラムについて説明する。図４は、本実施形態のデータ処理方法を示すフローチャートである。
　まず、計測ユニット６１が計測データｓ（ｘ’，ｙ’，０，ｋ）を取得する（ステップＳ１－１）。そして、データ処理ユニット６６は、計測データｓ（ｘ’，ｙ’，０，ｋ）に対して、ヒルベルト変換を行う（ステップＳ１－２）。これにより、各送受信点における周波数データの虚数成分が得られる。

　次に、データ処理ユニット６６は、計測データｓ（ｘ’，ｙ’，０，ｋ）に対して、（ｘ’，ｙ’）に関する２重フーリエ変換を行う（ステップＳ１－３）。これにより、式（１－６）に示されるように、Ｓ（ｋ_ｘ，ｋ_ｙ，０，ｋ）が得られる。
　次に、データ処理ユニット６６は、Ｓ（ｋ_ｘ，ｋ_ｙ，０，ｋ）に対して、変数置換を行う（ステップＳ１－４）。具体的には、式（１－４）を用いて、（ｋ_ｘ，ｋ_ｙ，ｋ）の関数を（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）の関数にする。これにより、Ｓ（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）が得られる。
　次に、データ処理ユニット６６は、Ｓ（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）に対して、（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）に対して３重逆フーリエ変換を行う（ステップＳ１－５）。これにより、式（１－８）に示されるように、反射率ｆ（ｘ，ｙ，ｚ）が得られる。

　記憶部６６ａは、本実施形態のデータ処理方法を実行するためのプログラムを記憶する。記憶部６６ａに記憶されたプログラムは、データ処理ユニット６６に、本実施形態のデータ処理方法を実行させる。

　ここで、フーリエ変換処理における折り返し雑音とされるエイリアスデータについて説明する。アレイアンテナで計測する場合、アンテナのサイズ、走査間隔をナイキスト基準以下にしなければ、計測データにエイリアシングが発生する。エイリアシングとは、異なる周波数成分の連続信号が標本化によって区別できなくなることをいう。

　図５は、エイリアシングの概念図を示す。空間的なサンプリング間隔から決定されるナイキスト波数をｋ_ｎｙｑと定義する。ナイキスト波数ｋ_ｎｙｑは、以下の式（１－１０）で表される。

　ここで、Δｘ、Δｙは、それぞれｘ軸方向、ｙ軸方向の計測間隔である。
　ナイキスト波数ｋ_ｎｙｑを超える信号は、エイリアス（aliases）と呼ばれ、信号処理の中で折り返し雑音として扱われる。

　図６は、第１実施形態の計測波形の一例を示す。ここで、電磁波の周波数を１５ＧＨｚ、比誘電率を１、計測間隔Δｙ = ７．５ｍｍ、計測幅５６２．５ｍｍ、点ターゲットの深さを原点直下１２０ｍｍとする。ここで、原点からの距離ｙの位置で測定される波形の波数ｋｙは、以下の式（１－１１）で表される。

　ここで、θ_tyは、最外縁送受信アンテナとターゲットを結ぶ線とｚ軸とのなす角度である。また、ナイキスト波数は、以下の式（１－１２）で表される。

　これより、ナイキスト基準は、以下の式（１－１３）で表される。

　この条件を満たさないｙ座標の波形、すなわち、図６におけるエイリアシング領域のデータは、全てエイリアスデータ（折り返し雑音）となる。

　上記の式から、電磁波の波長が短くなる、計測間隔が広がる、又は、ターゲットが計測面に近くなるほど、エイリアスデータが増加することが分かる。通常、エイリアスデータも信号処理の中で折り返し雑音として処理される。そのため、エイリアスデータの増加は、画像分解能の悪化、浅いターゲットの画像強度の低下、Ｓ／Ｎ比の悪化の原因となる。

　以下、計測間隔Δｘ＝１０ｍｍ、計測間隔Δｙ＝１９．２５ｍｍ、媒質中の最高周波数の波長λ_ｆｍａｘ＝２１．２４ｍｍとする。このとき、以下の式（１－１４）が得られる。

　このとき、ナイキスト基準を満たさず、ｘ軸方向、ｙ軸方向共に、エイリアシングが発生する。この場合のナイキスト波数は、以下の式（１－１５）で表される。

　最高周波数から求まる最大波数は、以下の式（１－１６）で表される。

　ここで、ステップＳ１－３におけるｘ軸に関するフーリエ変換は、図７で表される。また、ステップＳ１－３におけるｙ軸に関するフーリエ変換は、図８で表される。
　図７より、ｘ軸に関するフーリエ変換において、±ｋ_{ｘ，ｎｙｑ}上限とする波数が範囲Aに出力される。±ｋ_{ｘ，ｎｙｑ}を超える波数のデータは、エイリアスデータとして2ｋ_{ｘ，ｎｙｑ}を周期として折り返し現れる。よって、ステップＳ１－４の変数置換処理において、範囲Ａではなく、±ｋ_ｘｍａｘを上限とする範囲Ａ’で変数置換処理を行う。これにより、エイリアスデータを折り返し雑音として捨てるのではなく、実際に意味のあるデータとして変数置換処理に加えることにより、ｘ軸方向の分解能を向上させ、更には浅いターゲットの画像強度をより強くさせることが可能になる。
　ｙ軸に関しても、同様に、図８に示すように、範囲Ｂではなく範囲Ｂ’で変数置換処理を行う。

＜第２実施形態＞
　以下、第２実施形態のデータ処理方法、計測システム、及び、プログラムについて、詳細に説明する。第１実施形態の送信用アレイアンテナ５０及び受信用アレイアンテナ５２は、一方向（図３ではｙ方向）に配列されるが、本実施形態は、送信用アレイアンテナ５０及び受信用アレイアンテナ５２の配列が異なる。本実施形態の送信用アレイアンテナ５０及び受信用アレイアンテナ５２は、平面状に配置される。
　また、第１実施形態では、送信点と受信点の座標をいずれもｐ（ｘ’，ｙ’，ｚ’）としたが、本実施形態では、送信点と受信点の座標が異なる。本実施形態では、図９に示すように、送信点ｐ_１（ｘ’_１，ｙ’_１，ｚ’_１）、受信点ｐ_２（ｘ’_２，ｙ’_２，ｚ’_２）が、ｘｙ平面に配列される。

　ここで、測定対象物の反射点（ｘ，ｙ，ｚ）における反射率をｆ（ｘ，ｙ，ｚ）とする。送信点ｐ_１（ｘ’_１，ｙ’_１，ｚ’_１）及び受信点ｐ_２（ｘ’_２，ｙ’_２，ｚ’_２）における計測データをｓ（ｘ’_１，ｘ’_２，ｙ’_１，ｙ’_２，ｚ’_１，ｚ’_２，ｋ）とする。真空中の電磁波の伝播波長をλ_０とする。媒質の比誘電率をε_ｒとする。伝播する電磁波の波数をｋとする。

　このとき、計測データｓ（ｘ’_１，ｘ’_２，ｙ’_１，ｙ’_２，ｚ’_１，ｚ’_２，ｋ）は、以下の式で表せる。

　但し、

である。

　式（２－１）では、電磁波を球面波で表しており、距離減衰は省略されている。この距離減衰は、以降の処理を行う上で影響が小さいため、省略されている。式（２－１）中の被積分関数の指数部をフーリエ変換の表記で表すと、以下の式となる。これは、式（２－１）の球面波を３次元の平面波に分解することに等しい。

　ここで、（ｋ’_ｘ１，ｋ’_ｙ１，ｋ’_ｚ１）は、送信点から反射点までの間で伝搬する波動の球面波の波数ベクトルの成分である。また、（ｋ’_ｘ２，ｋ’_ｙ２，ｋ’_ｚ２）は、反射点から受信点までの間で伝搬する波動の球面波の波数ベクトルの成分である。但し、

　を満たす。

　以下、式（２－３）に基づいて、ｓ（ｘ’_１，ｘ’_２，ｙ’_１，ｙ’_２，ｚ’_１，ｚ’_２，ｋ）から反射率ｆ（ｘ，ｙ，ｚ）を導出する。まず、式（２－３）の両辺をｘ’_１，ｘ’_２，ｙ’_１，ｙ’_２に関して４重フーリエ変換を行う。

　式（２－５）の左辺を以下の式（２－６）のように書き換えて整理する。

　すると、式（２－５）は、式（２－７）で表される。

　式（２－７）の両辺に以下の積分を行う。

　ここで、以下の変数置換を行う。

　ここで、式（２－９）、式（２－１０）から、以下の式が得られる。

　この変数置換でのヤコビアンの絶対値｜Ｊ｜は式（２－１２）、式（２－１３）より以下の式でそれぞれ与えられる。

　式（２－９）、式（２－１０）、式（２－１４）、式（２－１５）を式（２－８）に代入して変数置換を行うと、以下の式が得られる。

　ここで、式（２－１６）の２行目右辺の（ｕ，ｖ）に関する積分は定数となるため、省略した。式（２－１６）の両辺に（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）について３重逆フーリエ変換を行うと、反射率ｆ（ｘ，ｙ，ｚ）が以下のように得られる。

　式（２－１７）を解くために、（ｋ’_ｚ１，ｋ’_ｚ２，ｋ）を（ｋ’_ｘ１，ｋ’_ｘ２，ｋ’_ｙ１，ｋ’_ｙ２，ｋ_ｚ）又は（ｋ_ｘ，ｕ，ｋ_ｙ，ｖ，ｋ_ｚ）で表す必要がある。
　式（２－４）、式（２－１１）を用いて整理し、また、送信点と受信点がいずれも原点を通るｘｙ平面上に位置する場合、ｚ’_１＝ｚ’_２＝０となるため、式（２－１７）は以下のように表せる。

　以上のように、データ処理ユニット６６は、計測データｓ（ｘ’_１，ｘ’_２，ｙ’_１，ｙ’_２，ｚ’_１，ｚ’_２，ｋ）に基づいて、反射率ｆ（ｘ，ｙ，ｚ）を求める。

　以下、図１０を参照して、本実施形態のデータ処理方法、及び、プログラムについて説明する。図１０は、本実施形態のデータ処理方法を示すフローチャートである。
　まず、計測ユニット６１が計測データｓ（ｘ’_１，ｘ’_２，ｙ’_１，ｙ’_２，ｚ’_１，ｚ’_２，ｋ）を取得する（ステップＳ２－１）。そして、データ処理ユニット６６は、計測データｓ（ｘ’_１，ｘ’_２，ｙ’_１，ｙ’_２，ｚ’_１，ｚ’_２，ｋ）に対して、ヒルベルト変換を行う（ステップＳ２－２）。これにより、各計測点における周波数データの虚数成分が得られる。

　次に、データ処理ユニット６６は、計測データｓ（ｘ’_１，ｘ’_２，ｙ’_１，ｙ’_２，ｚ’_１，ｚ’_２，ｋ）に対して、（ｘ’_１，ｘ’_２，ｙ’_１，ｙ’_２）に関する４重フーリエ変換を行う（ステップＳ２－３）。これにより、式（２－６）に示されるように、Ｓ（ｋ’_ｘ１，ｋ’_ｘ２，ｋ’_ｙ１，ｋ’_ｙ２，０，０，ｋ）が得られる。
　次に、データ処理ユニット６６は、Ｓ（ｋ’_ｘ１，ｋ’_ｘ２，ｋ’_ｙ１，ｋ’_ｙ２，０，０，ｋ）に対して、変数置換を行う（ステップＳ２－４）。具体的には、式（２－９）、式（２－１０）、式（２－１４）、式（２－１５）を用いて、（ｋ’_ｘ１，ｋ’_ｘ２，ｋ’_ｙ１，ｋ’_ｙ２，ｋ）の関数を（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）の関数にする。これにより、Ｓ（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）が得られる。
　次に、データ処理ユニット６６は、Ｓ（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）に対して、（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）に対して３重逆フーリエ変換を行う（ステップＳ２－５）。これにより、式（２－１８）に示されるように、反射率ｆ（ｘ，ｙ，ｚ）が得られる。

　ここで、電磁波の周波数を４．４６４８４ＧＨｚ、比誘電率を１０、計測間隔Δｘ’_１ = ３８．５ｍｍ、計測間隔Δｘ’_２ = ３８．５ｍｍ、計測間隔Δｙ’_１ = ３８．５ｍｍ、計測間隔Δｙ’_２ = ３８．５ｍｍ、計測幅６１６ｍｍ、媒質中の最高周波数の波長λ_ｆｍａｘ＝２１．２４ｍｍとする。このとき、以下の式（２－２０）が得られる。

　このとき、ナイキスト基準を満たさず、ｘ’_１軸方向、ｘ’_２軸方向、ｙ’_１軸方向、ｙ’_２軸方向のそれぞれに、エイリアシングが発生する。この場合のナイキスト波数は、以下の式（２－２１）で表される。

　最高周波数から求まる最大波数は、以下の式（２－２２）で表される。

　ここで、ステップＳ２－３におけるｘ’_１軸、ｘ’_２軸に関するフーリエ変換は、図１１で表される。また、ステップＳ２－３におけるｙ’_１軸、ｙ’_２軸に関するフーリエ変換は、図１２で表される。
　図１１より、ｘ’_１軸、ｘ’_２軸に関するフーリエ変換において、それぞれ、±ｋ’_{ｘ１，ｎｙｑ}、±ｋ’_{ｘ２，ｎｙｑ}を上限とする波数が範囲Ａに出力される。±ｋ’_{ｘ１，ｎｙｑ}、±ｋ’_{ｘ２，ｎｙｑ}を超える波数のデータは、エイリアスデータとして、それぞれ、２ｋ’_{ｘ１，ｎｙｑ}、２ｋ’_{ｘ２，ｎｙｑ}を周期として折り返し現れる。よって、ステップＳ２－４の変数置換処理において、範囲Ａではなく、それぞれ、±ｋ’_{ｘ１ｍａｘ}、±ｋ’_{ｘ２ｍａｘ}を上限とする範囲Ａ’で変数置換処理を行う。これにより、エイリアスデータを折り返し雑音として捨てるのではなく、実際に意味のあるデータとして変数置換処理に加えることにより、ｘ’_１軸、ｘ’_２軸方向の分解能を向上させ、更には浅いターゲットの画像強度をより強くさせることが可能になる。
　ｙ’_１軸、ｙ’_２軸に関しても、同様に、図１２に示すように、範囲Ｂではなく、それぞれ、±ｋ’_{ｙ１ｍａｘ}、±ｋ’_{ｙ２ｍａｘ}を上限とする範囲Ｂ’で変数置換処理を行う。

＜第３実施形態＞
　以下、第３実施形態のデータ処理方法、計測システム、及び、プログラムについて、詳細に説明する。第２実施形態では、送信用アレイアンテナ５０及び受信用アレイアンテナ５２は、平面状に配置されるが、本実施形態は、送信用アレイアンテナ５０及び受信用アレイアンテナ５２の配列が異なる。本実施形態の送信用アレイアンテナ５０及び受信用アレイアンテナ５２は、直線状に配置される。
　具体的には、本実施形態では、送信アンテナ１０ａ及び受信アンテナ１０ｂは、図１３に示すように、ｙ方向に配列される。送信用アレイアンテナ５０及び受信用アレイアンテナ５２の移動方向（走査方向）を、ｘ方向とする。送信用アレイアンテナ５０と受信用アレイアンテナ５２からみて、測定対象物のある方向（電磁波の送信方向）をｚ方向とする。
　なお、送信用アレイアンテナ５０及び受信用アレイアンテナ５２の移動方向（走査方向）をｙ方向としてもよい。

　したがって、測定対象物と送信用アレイアンテナ５０と受信用アレイアンテナ５２との位置関係は、図１３に示すように表すことができる。
　ここで、送信点の座標をｐ_１（ｘ’_１，ｙ’_１，ｚ’_１）、受信点の座標をｐ_２（ｘ’_２，ｙ’_２，ｚ’_２）とする。測定対象物の反射点（ｘ，ｙ，ｚ）における反射率をｆ（ｘ，ｙ，ｚ）とする。ｐ_２（ｘ’_２，ｙ’_２，ｚ’_２）における計測データをｓ（ｘ’_１，ｘ’_２，ｙ’_１，ｙ’_２，ｚ’_１，ｚ’_２，ｋ）とする。真空中の電磁波の伝播波長をλ_０とする。媒質の比誘電率をε_ｒとする。伝播する電磁波の波数をｋとする。

　このとき、計測データｓ（ｘ’_１，ｘ’_２，ｙ’_１，ｙ’_２，ｚ’_１，ｚ’_２，ｋ）は、以下の式で表せる。

　但し

である。

　式（３－１）では、電磁波を球面波で表しており、距離減衰は省略されている。この距離減衰は、以降の処理を行う上で影響が小さいため、省略されている。式（３－１）中の二段目の式の被積分関数の指数部をフーリエ変換の表記で表すと、以下の式となる。これは、式（３－１）の球面波を３次元の平面波に分解することに等しい。

　ここで、送信点ｐ_１（ｘ’_１，ｙ’_１，ｚ’_１）と受信点ｐ_２（ｘ’_２，ｙ’_２，ｚ’_２）のｘ座標が等しいことから、ｘ’＝ｘ’_１＝ｘ’_２とすると、式（３－３）は以下の式で表される。

　ここで、以下の変数置換を行う。

　式（３－６）から以下の式が得られる。

　式（３－７）からヤコビアンの絶対値を計算すると、以下の式が得られる。

　式（３－６）、式（３－８）を式（３－５）に代入して変数置換を行うと、以下の式が得られる。

　ここで、式（３－９）の２行目のｕに関する積分は、定数となるため省略した。式（３－９）の両辺に（ｘ，ｙ_１，ｙ_２）について３重逆フーリエ変換を行うと、以下の式が得られる。

　式（３－１０）の左辺を以下の式のように書き換えて整理する。

　すると、式（３－１０）は、以下の式で表される。

　式（３－１２）の両辺に以下の積分を行うと、以下の式が得られる。

　ここで、（ｋ’_ｙ１，ｋ’_ｙ２）、（ｋ’_ｚ１，ｋ’_ｚ２）に対して、以下の変数置換を定義する。

　ここで、式（３－１４）から、以下の式が得られる。

　この変数置換でのヤコビアンの絶対値は、以下の式で与えられる。

　式（３－１）、式（３－１５）、式（３－１７）を式（３－１３）に代入して、変数置換を行うと、以下の式が得られる。

　ここで、式（３－１８）の２行目右辺のｖに関する積分は、定数となるため省略した。式（３－１８）の両辺に（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）について３重逆フーリエ変換を行うと、反射率ｆ（ｘ，ｙ，ｚ）が以下のように得られる。

　ここで、原点を通るｘ’－ｙ’平面に計測面を合わせるため、ｚ’＝０とすると、式（３－１９）は、以下のように表される。

　式（３－２０）を解くために、ｋを（ｋ’_ｙ１，ｋ’_ｙ２，ｋ_ｚ）又は（ｋ_ｙ，ｖ，ｋ_ｚ）で表す必要がある。
　式（３－４）、式（３－６）、式（３－１５）、及び、仮定より得られる以下の式（３－２１）の４つの式の連立方程式を解く。

　これより、ｋは以下の式で表される。

　以下、図１４を参照して、本実施形態のデータ処理方法、及び、プログラムについて説明する。図１４は、本実施形態のデータ処理方法を示すフローチャートである。
　まず、計測ユニット６１が計測データｓ（ｘ’，ｙ’_１，ｙ’_２，０，０，ｋ）を取得する（ステップＳ３－１）。そして、データ処理ユニット６６は、計測データｓ（ｘ’，ｙ’_１，ｙ’_２，０，０，ｋ）に対して、ヒルベルト変換を行う（ステップＳ３－２）。これにより、各計測点における周波数データの虚数成分が得られる。

　次に、データ処理ユニット６６は、計測データｓ（ｘ’，ｙ’_１，ｙ’_２，０，０，ｋ）に対して、（ｘ’，ｙ’_１，ｙ’_２）に関する３重フーリエ変換を行う（ステップＳ３－３）。これにより、式（３－１１）に示されるように、Ｓ（ｋ_ｘ，ｋ’_ｙ１，ｋ’_ｙ２，０，０，ｋ）が得られる。
　次に、データ処理ユニット６６は、Ｓ（ｋ_ｘ，ｋ’_ｙ１，ｋ’_ｙ２，０，０，ｋ）に対して、変数置換を行う（ステップＳ３－４）。具体的には、式（３－１４）、式（３－１５）を用いて、（ｋ_ｘ，ｋ’_ｙ１，ｋ’_ｙ２，ｋ）の関数を（ｋ_ｘ，ｋ_ｙ，ｖ，ｋ）の関数にする。これにより、Ｓ（ｋ_ｘ，ｋ_ｙ，ｖ，０，０，ｋ）が得られる。
　次に、データ処理ユニット６６は、Ｓ（ｋ_ｘ，ｋ_ｙ，ｖ，０，０，ｋ）に対して、（ｋ_ｘ，ｋ_ｙ，ｋ_ｚ）に対して３重逆フーリエ変換を行う（ステップＳ３－５）。これにより、式（３－２０）に示されるように、反射率ｆ（ｘ，ｙ，ｚ）が得られる。

　ここで、電磁波の周波数を４．４６４８４ＧＨｚ、比誘電率を１０、計測間隔Δｘ = １０ｍｍ、計測間隔Δｙ’_１ = ３８．５ｍｍ、計測間隔Δｙ’_２ = ３８．５ｍｍ、計測幅６１６ｍｍ、媒質中の最高周波数の波長λ_ｆｍａｘ＝２１．２４ｍｍとする。このとき、以下の式（３－２３）が得られる。

　このとき、ナイキスト基準を満たさず、ｘ軸方向、ｙ’_１軸方向、ｙ’_２軸方向のそれぞれに、エイリアシングが発生する。この場合のナイキスト波数は、以下の式（３－２４）で表される。

　最高周波数から求まる最大波数は、以下の式（３－２５）で表される。

　ここで、ステップＳ３－３におけるｘ軸に関するフーリエ変換は、図１５で表される。また、ステップＳ３－３におけるｙ’_１軸、ｙ’_２軸に関するフーリエ変換は、図１６で表される。
　図１５より、ｘ軸に関するフーリエ変換において、それぞれ、±ｋ_{ｘ，ｎｙｑ}、±ｋ’_{ｘ２，ｎｙｑ}を上限とする波数が範囲Ａに出力される。±ｋ_{ｘ，ｎｙｑ}を超える波数のデータは、エイリアスデータとして、２ｋ_{ｘ，ｎｙｑ}を周期として折り返し現れる。よって、ステップＳ３－４の変数置換処理において、範囲Ａではなく、それぞれ、±ｋ_ｘｍａｘを上限とする範囲Ａ’で変数置換処理を行う。これにより、エイリアスデータを折り返し雑音として捨てるのではなく、実際に意味のあるデータとして変数置換処理に加えることにより、ｘ軸方向の分解能を向上させ、更には浅いターゲットの画像強度をより強くさせることが可能になる。
　ｙ’_１軸、ｙ’_２軸に関しても、第２実施形態と同様に、図１６に示すように、範囲Ｂではなく、それぞれ、±ｋ’_{ｙ１ｍａｘ}、±ｋ’_{ｙ２ｍａｘ}を上限とする範囲Ｂ’で変数置換処理を行う。

（シミュレーション結果）
　以下、第３実施形態のデータ処理方法をコンピュータシミュレーションした結果について説明する。シミュレーション条件は、以下の通りである。

・使用周波数帯域ｆ_ｍｉｎ～ｆ_ｍａｘ：ＤＣ～４．４６４８４ＧＨｚ
・中心周波数ｆｃ：２．２３２４２ＧＨｚ
・中心周波数の真空中の波長λ_０ｃ：１３４．３８ｍｍ
・媒質の比誘電率ε_ｒ：１０
・媒質中の中心周波数の波長λｃ：４２．５０ｍｍ
・媒質中の最高周波数の波長λ_ｆｍａｘ：２１．２４ｍｍ
・計測数Ｎｘ：１２８ｐｔ
・計測間隔Δｘ：１０ｍｍ
・計測幅：１２８０ｍｍ
・計測数Ｎ_ｙ’１：１６ｐｔ
・計測間隔Δｙ’１：３８．５ｍｍ
・計測幅：６１６ｍｍ
・計測数Ｎ_ｙ’２：１６ｐｔ
・計測間隔Δｙ’２：３８．５ｍｍ
・計測幅：６１６ｍｍ
・点ターゲット１の座標（単位：ｍｍ）：（６４０，３０７，２０）
・点ターゲット２の座標（単位：ｍｍ）：（６４０，３０７，２００）

　図１７は、比較例１のデータ処理方法でシミュレーションした点ターゲット１のｘｚ平面画像を示す。図１８は、実施例１のデータ処理方法でシミュレーションした点ターゲット１のｘｚ平面画像を示す。図１９は、比較例１のデータ処理方法でシミュレーションした点ターゲット１のｙｚ平面画像を示す。図２０は、実施例１のデータ処理方法でシミュレーションした点ターゲット１のｙｚ平面画像を示す。
　図２１は、比較例１と実施例１の点ターゲット１のｘ軸方向の波形を示す。図２２は、比較例１と実施例１の点ターゲット１のｙ軸方向の波形を示す。図２３は、比較例１と実施例１の点ターゲット１のｚ軸方向の波形を示す。図２１～図２３において、いずれも実施例１は実線で示され、比較例１は破線で示される。

　図２４は、比較例１のデータ処理方法でシミュレーションした点ターゲット２のｘｚ平面画像を示す。図２５は、実施例１のデータ処理方法でシミュレーションした点ターゲット２のｘｚ平面画像を示す。図２６は、比較例１のデータ処理方法でシミュレーションした点ターゲット２のｙｚ平面画像を示す。図２７は、実施例１のデータ処理方法でシミュレーションした点ターゲット２のｙｚ平面画像を示す。
　図２８は、比較例１と実施例１の点ターゲット２のｘ軸方向の波形を示す。図２９は、比較例１と実施例１の点ターゲット２のｙ軸方向の波形を示す。図３０は、比較例１と実施例１の点ターゲット２のｚ軸方向の波形を示す。図２８～図３０において、いずれも実施例１は実線で示され、比較例１は破線で示される。

　比較例１では、本実施形態におけるステップＳ３－４の変数置換処理において、範囲Ａで変数置換処理を行った。実施例１では、本実施形態におけるステップＳ３－４の変数置換処理において、範囲Ａではなく、それぞれ、±ｋ_ｘｍａｘを上限とする範囲Ａ’で変数置換処理を行った。以下、実施例１における±ｋ_ｘｍａｘを上限とする範囲Ａ’での変数置換処理を超分解能処理と呼ぶ。

　図１７～図３０より、比較例1に比べて、超分解能処理を行う実施例１によれば、点ターゲット１、２のいずれにおいても、ｘ軸方向、ｙ軸方向の分解能が向上することが確認された。

＜第４実施形態＞
　本実施形態では、第３実施形態において生じ得るノイズを低減する。ｙ軸方向に関する計測データｓは、送信アンテナからターゲットまで伝播する電磁波の振幅をｔ、ターゲットの反射率を１、ターゲットから受信アンテナまで伝播する電磁波の振幅をｒとすると、以下の式（４－１）で表される。

　このとき、ナイキスト周波数以下の振幅成分をｔ_０、ｒ_０とし、ｎ回の折り返し振幅成分（エイリアスデータ）をｔ_ｎ、ｒ_ｎとすると、ｔとｒは、それぞれ以下の式（４－２）で表される。

　式（４－２）を式（４－１）に代入すると、以下の式（４－３）が得られる。

　第３実施形態で示した超分解能処理は、例えば、ステップＳ３－４の変数置換時に、図１６で示される（ｋ’_y1，ｋ’_y2）空間において、以下の式（４－４）で示されるように、送信アンテナの切替えに起因するｋ’_y1と、受信アンテナに起因するｋ’_y2がどちらもナイキスト波数以下の範囲内にある場合には、エイリアスデータ（折り返し雑音成分）を含む成分は、信号処理の中で全て雑音成分となる。

　このとき、式（４－３）は、式（４－５）のように表される。

　ここで、実線で囲んだｒ_０ｔ_０は信号成分であり、破線で囲んだ部分はノイズ成分である。

　次に、変数置換時にｋ’_y1のみが１回の折り返し雑音領域に入る場合には、式（４－６）の条件となり、式（４－３）は式（４－７）で表される。

　ここで、実線で囲んだｒ_０ｔ_１は信号成分であり、破線で囲んだ部分はノイズ成分である。

　次に、変数置換時にｋ’_y２が２回の折り返し雑音領域に入る場合には、式（４－８）の条件となり、式（４－３）は式（４－９）で表される。

　ここで、実線で囲んだｒ_２ｔ_１は信号成分であり、破線で囲んだ部分はノイズ成分である。

　ここで、図５に示したのと同様、送信アンテナ計測点のy’₁座標、受信アンテナ計測点のy’₂座標が、ターゲットに近い場合、すなわち、ターゲットの直上付近にある場合、y’₁、y’₂軸方向の計測データの空間周波数は、それぞれ低くなる。反対に、y’₁座標、y’₂座標が、ターゲットから遠くなるに従い、y’₁、y’₂軸方向の計測データの空間周波数は、それぞれ高くなる。そして、その地点の波数がナイキスト波数を超えたときに、折り返し雑音（エイリアス）成分として計測データに含まれる。このように、エイリアスデータ成分（折り返し雑音成分）は、相対的にターゲットから遠い地点で発生する傾向がある。そのため、距離減衰によりナイキスト波数以下のデータ成分と比較して、振幅の大きさが小さくなる。そのため、以下の式（４－１０）が一般的に成立する。

　これより、ｒ_ｎとｔ_ｎに関する２次の項の大きさは、ｒ_０或いはｔ_０含む項に対して小さくなる。これは、以下の式（４－１１）で表される。

　更に、減衰率が小さな空中とは異なり、減衰率が大きなコンクリート内部や地中を電磁波が伝播する場合には、距離減衰に加えて媒質に起因する減衰率によりｒ_ｎとｔ_ｎの大きさは、より小さくなる。このとき、式（４－１１）は、以下の式（４－１２）のように表される。

　よって、式（４－５）、式（４－７）、式（４－９）、式（４－１２）からｔ_ｎとｒ_ｎに関する２次の項は、それ以外の項に対して相対的にＳ／Ｎ比が悪いデータとなる。そのため、Ｓ／Ｎ比を改善するためには、ｒ_ｎとｔ_ｎに関する２次の項を、ステップＳ３－４の変数置換処理において、選択的に除外すれば良い。このとき、式（４－３）は、以下の式（４－１３）で表される。

　ここで、実線で囲んだ成分のみを変数置換処理で利用し、破線で囲んだ部分は変数置換処理で利用しない。

　このように、減衰率が大きな媒質を電磁波が伝播する場合には、変数置換処理時に式（４－１３）で表されるフィルタを利用することにより、3次元画像のＳ／Ｎ比の劣化を抑制できる。これは（ｋ’_y1，ｋ’_y2）空間における「十字型ｋ’_y1－ｋ’_y2フィルタ」として表される。図３１は、十字型ｋ’_y1－ｋ’_y2フィルタの概念図である。

　図１５、図１６に示されるｋ_ｘｍａｘ、ｋ’_{ｙ１ｍａｘ}、ｋ’_{ｙ２ｍａｘ}は、ターゲットの位置あるいはアンテナの指向性に依存する。ここで、計測面近傍（ｚ＝０）のターゲットまで計測対象に含めば、式（４－１４）で示すように、アンテナの指向性θ_ｂｘ，θ_ｂｙによって最大値が決定されてもよい。

　但し、

である。

（シミュレーション結果１）
　以下、第４実施形態のデータ処理方法をコンピュータシミュレーションした結果について説明する。まず、第３実施形態と同じシミュレーション条件で、浅深度ターゲットの画像強度の改善効果を検証する。
　超分解能処理では、失われていたエイリアスデータを用いる。そのため、超分解能処理が無い場合と比較して、ターゲット画像強度（画像振幅）が増加する。この効果は、深深度ターゲットよりも、エイリアスデータの多い浅深度ターゲットの方が大きい。すなわち、エイリアシングによって低減していた浅深度ターゲットの画像強度は、超分解能処理によって大きく回復する。

　図３２は、４つの点ターゲットのシミュレーション結果を示す。各点ターゲットの水平位置は計測面中央で、深さは、２０ｍｍ、２００ｍｍ、４００ｍｍ、６００ｍｍである。図３２において、超分解能処理を行った実施例１は実線で示され、超分解能処理を行わない比較例１は破線で示される。3次元画像振幅の大きさは、深さ２０ｍｍのターゲットの超分解能処理を行わない振幅を１として正規化されている。

　ターゲットの深さが浅いほど、超分解能処理の効果でターゲットの画像振幅が強くなる。シミュレーションデータでは距離減衰を考慮しているため、深いターゲットの振幅が非常に小さくなり、回復効果を直接比較することが難しい。そのため、同じ深さのターゲットの振幅について、超分解能処理を行わない３次元画像振幅に対する超分解能処理を行った３次元画像振幅の比を回復率として定義する。図３３は、各点ターゲットの回復率を示す。図３３より、ターゲットが浅いほど、エイリアスデータが多いため、超分解能処理による回復率が大きいことが分かる。深さ２０ｍｍでは、回復率は６倍に達する。一方、ターゲットが深いほど、エイリアスデータは少なくなり、回復率は小さくなる。深さ６００ｍｍでは、回復率は１．２倍程度である。

　以上より、超分解能処理によってエイリアスデータを用いることにより、画像振幅強度が改善することが確認された。特に、エイリアスデータの多い浅いターゲットの画像において、顕著な改善が見られる。これは、主に、浅深度ターゲット画像のＳ／Ｎ比の改善として３次元画像に反映される。

（シミュレーション結果２）
　以下、第３実施形態、第４実施形態のデータ処理方法をコンピュータシミュレーションした結果について説明する。シミュレーション条件は、以下の通りである。
・使用周波数帯域ｆ_ｍｉｎ～ｆ_ｍａｘ：ＤＣ～４．４６４８４ＧＨｚ
・中心周波数ｆｃ：２．２３２４２ＧＨｚ
・中心周波数の真空中の波長λ_０ｃ：１３４．３８ｍｍ
・媒質の比誘電率ε_ｒ：５
・媒質中の中心周波数の波長λｃ：６０．１ｍｍ
・媒質中の最高周波数の波長λ_ｆｍａｘ：３０．０５ｍｍ
・走査方向の計測間隔Δｘ：１０ｍｍ
・計測数Ｎｘ：１１５ｐｔ
・計測幅：１１５０ｍｍ
・送信アンテナの計測間隔Δｙ１：３８．５ｍｍ
・送信アンテナ数Ｎ_ｙ１：１２
・受信アンテナの計測間隔Δｙ２：３８．５ｍｍ
・受信アンテナ数Ｎ_ｙ２：１２ｐｔ
・計測幅：４８１．２５ｍｍ

　以下、鉄筋コンクリート試験体を計測した結果を示す。計測に用いた送信アンテナと受信アンテナのレイアウトを、図３４に示す。試験体は、ｚ座標の異なる複数の鉄筋を内部に有する。図３５は、超分解能処理を行わない、試験体の３次元画像である。図３６は、超分解能処理を行った、試験体の３次元画像である。図３５、図３６において、複数の鉄筋のｚ座標が示されている。図３５と図３６を比較すると、超分解能処理を行った図３６の方が、全体的に鉄筋の画像が細くなっており、画像分解能の向上が確認できた。また、深さ３５ｍｍ、４５ｍｍの鉄筋の画像強度が強くなっており、浅深度ターゲットの画像振幅の回復効果も確認できた。一方、鉄筋と鉄筋の間に干渉縞のような不要なノイズ（エイリアシングノイズ）が新たに発生していることも確認できた。

　次に、第４実施形態、図３１で説明した十字型ｋ’_y1－ｋ’_y2フィルタを適用した。図３７は、十字型ｋ’_y1－ｋ’_y2フィルタを適用して超分解能処理を行った、試験体の３次元画像である。十字型ｋ’_y1－ｋ’_y2フィルタを適用しないで超分解能処理を行った、試験体の３次元画像である図３６と、十字型ｋ’_y1－ｋ’_y2フィルタを適用して超分解能処理を行った図３７を比較すると、十字型ｋ’_y1－ｋ’_y2フィルタにより、鉄筋画像の分解能や浅深度鉄筋画像の強度の劣化を抑制しつつ、干渉縞ノイズ（エイリアシングノイズ）を低減できることが確認できた。

　続いて、十字型ｋ’_y1－ｋ’_y2フィルタが分解能に与える影響を、シミュレーションデータを用いて評価した。シミュレーションの条件は、第３実施形態と同様である。
　図３８は、比較例２と実施例２の点ターゲット１のｘ軸方向の波形を示す。図３９は、比較例２と実施例２の点ターゲット１のｙ軸方向の波形を示す。図４０は、比較例２と実施例２の点ターゲット１のｚ軸方向の波形を示す。図４１は、比較例２と実施例２の点ターゲット２のｘ軸方向の波形を示す。図４２は、比較例２と実施例２の点ターゲット２のｙ軸方向の波形を示す。図４３は、比較例２と実施例２の点ターゲット２のｚ軸方向の波形を示す。
　図３８～図４３において、いずれも実施例２は実線で示され、比較例２は破線で示される。

　比較例２では、第３実施形態におけるステップＳ３－４の変数置換処理において、十字型ｋ’_y1－ｋ’_y2フィルタを適用しないで、それぞれ、±ｋ_ｘｍａｘを上限とする範囲Ａ’で変数置換処理を行った。実施例２では、第３実施形態におけるステップＳ３－４の変数置換処理において、十字型ｋ’_y1－ｋ’_y2フィルタを適用して、それぞれ、±ｋ_ｘｍａｘを上限とする範囲Ａ’で変数置換処理を行った。
　図３８～図４３より、十字型ｋ’_y1－ｋ’_y2フィルタを適用して超分解能処理を行っても、ｘ軸、ｙ軸、ｚ軸のいずれも画像分解能にほとんど影響を与えていないことが確認された。

＜第５実施形態＞
　本実施形態では、第２実施形態において生じ得るノイズを低減する。第４実施形態では、第３実施形態において、十字型ｋ’_y1－ｋ’_y2フィルタを適用して超分解能処理を行った。本実施形態では、第２実施形態において、十字型ｋ’_ｘ1－ｋ’_ｘ2フィルタと十字型ｋ’_y1－ｋ’_y2フィルタを適用して、超分解能処理を行う。図４４は、十字型ｋ’_ｘ1－ｋ’_ｘ2フィルタの概念図である。すなわち、本実施形態では、第２実施形態において、図４４に示す十字型ｋ’_ｘ1－ｋ’_ｘ2フィルタと図３１に示す十字型ｋ’_y1－ｋ’_y2フィルタの両方を適用して超分解能処理を行う。

　ｋ’_{ｘ１ｍａｘ}、ｋ’_{ｘ２ｍａｘ}、ｋ’_{ｙ１ｍａｘ}、ｋ’_{ｙ２ｍａｘ}は、ターゲットの位置あるいはアンテナの指向性に依存する。ここで、計測面近傍（ｚ＝０）のターゲットまで計測対象に含めば、式（５－１）で示すように、アンテナの指向性θ_ｂｘ，θ_ｂｙによって最大値が決定されてもよい。

　但し、

である。

１０ａ　送信アンテナ
１０ｂ　受信アンテナ
５０　送信用アレイアンテナ
５２　受信用アレイアンテナ
６０　レーダ装置
６１　計測ユニット
６４　システム制御回路
６６　データ処理ユニット
６８　画像表示ユニット
５８、５９　高周波スイッチ
６２　高周波回路
６９　エンコーダ

Claims

　物体に放射した波動の散乱波を解析するデータ処理方法であって、
　ｘｙ平面上に２次元に配列された複数の送信点p₁(x’₁, y’₁, z’₁)から、前記物体に前記波動を放射し、
　前記物体上の反射点(x, y, z)において反射率f(x, y, z)で反射した前記散乱波を、ｘｙ平面上に２次元に配列された複数の受信点p₂(x’₂, y’₂, z’₂)で計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)として受信し、
　前記計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)を式（１）より４重フーリエ変換してS(k’_x1, k’_x2, k’_y1, k’_y2, z’₁, z’₂, k)を求め、

　ｘ方向の計測間隔で定まるx’₁方向のナイキスト波数をk’_x1,nyqとし、x’₂方向のナイキスト波数をk’_x2,nyqとすると、-k_x1≦k’_x1≦k_x1、かつ、-k_x2≦k’_x2≦k_x2（但し、k’_x1,nyq≦k_x1、かつ、k’_x2,nyq≦k_x2）の範囲でｘ方向の変数置換処理を行い、
　ｙ方向の計測間隔で定まるy’₁方向のナイキスト波数をk’_y1,nyqとし、y’₂方向のナイキスト波数をk’_y2,nyqとすると、-k_y1≦k’_y1≦k_y1、かつ、-k_y2≦k’_y2≦k_y2（但し、k’_y1,nyq≦k_y1、かつ、k’_y2,nyq≦k_y2）の範囲でｙ方向の変数置換処理を行い、
　式（２）より３重逆フーリエ変換して、前記反射率f(x, y, z)を求める、

　データ処理方法。
但し、
z’₁= z’₂ = 0
kは、伝播する前記波動の波数、
k’_x1, k’_y1, k’_z1は、前記送信点p₁(x’₁, y’₁, z’₁)から前記反射点(x, y, z)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k’_x2, k’_y2, k’_z2は、前記反射点(x, y, z)から前記受信点p₂(x’₂, y’₂, z’₂)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k_x = k’_x1 + k’_x2, u = k’_x1 - k’_x2, k_y = k’_y1 + k’_y2, v = k’_y1 - k’_y2
である。
　式（３）で定まるx’₁方向の最大波数をk’_x1maxとし、x’₂方向の最大波数をk’_x2maxとすると、-k’_x1max≦k’_x1≦k’_x1max、かつ、-k’_x2max≦k’_x2≦k’_x2maxの範囲でｘ方向の変数置換処理を行い、

　式（４）で定まるy’₁方向の最大波数をk’_y1maxとし、y’₂方向の最大波数をk’_y2maxとすると、-k’_y1max≦k’_y1≦k’_y1max、かつ、-k’_y2max≦k’_y2≦k’_y2maxの範囲でｙ方向の変数置換処理を行う、

　請求項１に記載のデータ処理方法。
　前記ｘ方向の変数置換処理において、k’_x1,nyq≦|k’_x1|≦k’_x1max、かつ、k’_x2,nyq≦|k’_x2|≦k’_x2maxの範囲を選択的に除外して、ｘ方向の前記変数置換処理を行い、
　前記ｙ方向の変数置換処理において、k’_y1,nyq≦|k’_y1|≦k’_y1max、かつ、k’_y2,nyq≦|k’_y2|≦k’_y2maxの範囲を選択的に除外して、ｙ方向の前記変数置換処理を行う、
　請求項２に記載のデータ処理方法。
　前記波動の波数と波数ベクトルの成分は、式（５）を満たす、
　請求項１～３のいずれかに記載のデータ処理方法。
　物体に放射した波動の散乱波を解析する計測システムであって、
　送受信部であって、
　　ｘｙ平面上に２次元に配列された複数の送信点p₁(x’₁, y’₁, z’₁)から、前記物体に前記波動を放射する送信部と、
　　前記物体上の反射点(x, y, z)において反射率f(x, y, z)で反射した前記散乱波を、ｘｙ平面上に２次元に配列された複数の受信点p₂(x’₂, y’₂, z’₂)で計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)として受信する受信部と、
　を有する送受信部と、
　処理装置であって、
　　前記計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)を式（１）より４重フーリエ変換してS(k’_x1, k’_x2, k’_y1, k’_y2, z’₁, z’₂, k)を求める手順と、

　ｘ方向の計測間隔で定まるx’₁方向のナイキスト波数をk’_x1,nyqとし、x’₂方向のナイキスト波数をk’_x2,nyqとすると、-k_x1≦k’_x1≦k_x1、かつ、-k_x2≦k’_x2≦k_x2（但し、k’_x1,nyq≦k_x1、かつ、k’_x2,nyq≦k_x2）の範囲でｘ方向の変数置換処理を行う手順と、
　ｙ方向の計測間隔で定まるy’₁方向のナイキスト波数をk’_y1,nyqとし、y’₂方向のナイキスト波数をk’_y2,nyqとすると、-k_y1≦k’_y1≦k_y1、かつ、-k_y2≦k’_y2≦k_y2（但し、k’_y1,nyq≦k_y1、かつ、k’_y2,nyq≦k_y2）の範囲でｙ方向の変数置換処理を行う手順と、
　式（２）より３重逆フーリエ変換して、前記反射率f(x, y, z)を求める手順と、

　を実行する処理装置と、
　を有する、計測システム。
但し、
z’₁= z’₂ = 0
kは、伝播する前記波動の波数、
k’_x1, k’_y1, k’_z1は、前記送信点p₁(x’₁, y’₁, z’₁)から前記反射点(x, y, z)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k’_x2, k’_y2, k’_z2は、前記反射点(x, y, z)から前記受信点p₂(x’₂, y’₂, z’₂)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k_x = k’_x1 + k’_x2, u = k’_x1 - k’_x2, k_y = k’_y1 + k’_y2, v = k’_y1 - k’_y2
である。
　前記処理装置は、
　　式（３）で定まるx’₁方向の最大波数をk’_x1maxとし、x’₂方向の最大波数をk’_x2maxとすると、-k’_x1max≦k’_x1≦k’_x1max、かつ、-k’_x2max≦k’_x2≦k’_x2maxの範囲でｘ方向の変数置換処理を行い、

　　式（４）で定まるy’₁方向の最大波数をk’_y1maxとし、y’₂方向の最大波数をk’_y2maxとすると、-k’_y1max≦k’_y1≦k’_y1max、かつ、-k’_y2max≦k’_y2≦k’_y2maxの範囲でｙ方向の変数置換処理を行う、

　請求項５に記載の計測システム。
　前記処理装置は、
　　前記ｘ方向の変数置換処理において、k’_x1,nyq≦|k’_x1|≦k’_x1max、かつ、k’_x2,nyq≦|k’_x2|≦k’_x2maxの範囲を選択的に除外して、ｘ方向の前記変数置換処理を行う手順と、
　　前記ｙ方向の変数置換処理において、k’_y1,nyq≦|k’_y1|≦k’_y1max、かつ、k’_y2,nyq≦|k’_y2|≦k’_y2maxの範囲を選択的に除外して、ｙ方向の前記変数置換処理を行う手順と、
　を実行する、
　請求項６に記載の計測システム。
　前記波動の波数と波数ベクトルの成分は、式（５）を満たす、
　請求項５～７のいずれかに記載の計測システム。
　物体に放射した波動の散乱波を解析するプログラムであって、
　計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)を式（１）より４重フーリエ変換してS(k’_x1, k’_x2, k’_y1, k’_y2, z’₁, z’₂, k)を求める手順と、

　ｘ方向の計測間隔で定まるx’₁方向のナイキスト波数をk’_x1,nyqとし、x’₂方向のナイキスト波数をk’_x2,nyqとすると、-k_x1≦k’_x1≦k_x1、かつ、-k_x2≦k’_x2≦k_x2（但し、k’_x1,nyq≦k_x1、かつ、k’_x2,nyq≦k_x2）の範囲でｘ方向の変数置換処理を行う手順と、
　ｙ方向の計測間隔で定まるy’₁方向のナイキスト波数をk’_y1,nyqとし、y’₂方向のナイキスト波数をk’_y2,nyqとすると、-k_y1≦k’_y1≦k_y1、かつ、-k_y2≦k’_y2≦k_y2（但し、k’_y1,nyq≦k_y1、かつ、k’_y2,nyq≦k_y2）の範囲でｙ方向の変数置換処理を行う手順と、
　式（２）より３重逆フーリエ変換して、反射率f(x, y, z)を求める手順と、

　をコンピュータに実行させるプログラム。
但し、
z’₁= z’₂ = 0
kは、伝播する前記波動の波数、
k’_x1, k’_y1, k’_z1は、送信点p₁(x’₁, y’₁, z’₁)から反射点(x, y, z)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k’_x2, k’_y2, k’_z2は、前記反射点(x, y, z)から受信点p₂(x’₂, y’₂, z’₂)の間で伝播する前記波動の球面波の波数ベクトルの成分、
k_x = k’_x1 + k’_x2, u = k’_x1 - k’_x2, k_y = k’_y1 + k’_y2, v = k’_y1 - k’_y2
である。
　　式（３）で定まるx’₁方向の最大波数をk’_x1maxとし、x’₂方向の最大波数をk’_x2maxとすると、-k’_x1max≦k’_x1≦k’_x1max、かつ、-k’_x2max≦k’_x2≦k’_x2maxの範囲でｘ方向の変数置換処理を行う手順と、

　　式（４）で定まるy’₁方向の最大波数をk’_y1maxとし、y’₂方向の最大波数をk’_y2maxとすると、-k’_y1max≦k’_y1≦k’_y1max、かつ、-k’_y2max≦k’_y2≦k’_y2maxの範囲でｙ方向の変数置換処理を行う手順と、

　をコンピュータに更に実行させる、請求項９に記載のプログラム。
　前記ｘ方向の変数置換処理において、k’_x1,nyq≦|k’_x1|≦k’_x1max、かつ、k’_x2,nyq≦|k’_x2|≦k’_x2maxの範囲を選択的に除外して、ｘ方向の前記変数置換処理を行う手順と、
　前記ｙ方向の変数置換処理において、k’_y1,nyq≦|k’_y1|≦k’_y1max、かつ、k’_y2,nyq≦|k’_y2|≦k’_y2maxの範囲を選択的に除外して、ｙ方向の前記変数置換処理を行う手順と、
　をコンピュータに更に実行させる、請求項１０に記載のプログラム。
　前記計測値s(x’₁, x’₂, y’₁, y’₂, z’₁, z’₂, k)は、ｘｙ平面上に２次元に配列された複数の送信点p₁(x’₁, y’₁, z’₁)から、前記物体に放射された波動が、前記物体上の前記反射点(x, y, z)において前記反射率f(x, y, z)で反射した前記散乱波を、ｘｙ平面上に２次元に配列された複数の受信点p₂(x’₂, y’₂, z’₂)で受信した値である、
　請求項９～１１のいずれかに記載のプログラム。
　前記波動の波数と波数ベクトルの成分は、式（５）を満たす、
　請求項９～１２のいずれかに記載のプログラム。