JP7802182B2 - Data processing method, measurement system, and program - Google Patents

Description

The present invention relates to a data processing method, measurement system, and program for using a computer to process measurement data of waves whose values are determined by the frequency of waves such as electromagnetic waves generated in space and the spatial coordinates of the space.

Radar devices are known for non-destructively inspecting the interior of non-metallic structures such as concrete and wood. Conventional radar devices have an array antenna, in which multiple antennas are arranged on a flat surface. The array antenna, for example, has a configuration in which antennas such as planar antennas are aligned in one direction, with the transmitting array antenna and the receiving array antenna arranged in close proximity. Furthermore, to accurately measure the interior of a structure, radar devices measure the object at a wide frequency band while changing the frequency of the electromagnetic waves at set frequency intervals.

Regarding radar devices with array antennas, for example, radar devices in which a transmitting array antenna and a receiving array antenna composed of multiple planar antennas are formed on a common dielectric substrate are known (JP 2015-095840 A, hereinafter referred to as "Patent Document 1"). In conventional radar devices, the arrangement direction of the planar antennas of the transmitting array antenna is parallel to the arrangement direction of the planar antennas of the receiving array antenna. The position of the planar antennas of the receiving array antenna in the arrangement direction is midway between the positions of two adjacent planar antennas of the transmitting array antenna.

In addition, a scattering tomography method is known that can perform general-purpose and high-speed analysis of inverse problems and easily visualize information about the interior of an object (Patent Publication No. 6557747, hereinafter referred to as "Patent Document 2").

Synthetic aperture processing is used to visualize the interior of structures from the measured data. There are two main types of synthetic aperture processing: addition methods such as the diffraction stacking method, and Fourier transform methods such as the FK migration method.
To achieve a practical calculation time, synthetic aperture processing using Fourier transform is practical.

In radar equipment, high spatial resolution in measurements is desirable to accurately inspect the interior of structures. Generally, the spatial resolution of data obtained by emitting waves with frequencies such as electromagnetic waves is determined by the wavelength of the center frequency of the waves when the distance between the structure to be measured and the measurement surfaces of the transmitting array antenna and receiving array antenna is relatively short and the measurement interval for the measurement data is small. Here, the distance between the structure to be measured and the measurement surfaces of the transmitting array antenna and receiving array antenna is, for example, less than one-quarter of the array antenna arrangement length. Furthermore, spatial resolution is the resolution within the plane in which each antenna of the array antenna is arranged.

For example, when the measurement interval of measurement data along the plane on which each antenna of an array antenna is arranged is sufficiently small, the theoretical spatial resolution is one-fourth of the wavelength of the wave at the center frequency of the frequency band when the frequency of the radiated wave is swept over a frequency band, taking into account the round-trip path of the electromagnetic wave. However, when the measurement interval of the measurement data becomes large and exceeds one-fourth of the minimum wavelength of the electromagnetic wave, the spatial resolution in actual measurements becomes larger than the ideal spatial resolution. In some cases, the spatial resolution in actual measurements becomes the measurement interval.

Conventional radar devices measure electromagnetic waves over a wide frequency band, from low to high frequencies, resulting in a long maximum wavelength. This requires each antenna in an array antenna to be large, and the length of the antenna arrangement in the array antenna is long. As a result, the spacing between antennas in a transmitting/receiving array antenna becomes long, and the measurement interval for measurement data tends to exceed one-quarter of the minimum wavelength of the emitted electromagnetic waves. This reduces spatial resolution below the theoretical resolution, and the measurement data is prone to aliasing. To achieve the ideal theoretical spatial resolution, i.e., one-quarter of the wavelength of the electromagnetic waves at the center frequency, the number of antennas in the transmitting/receiving array antenna must be increased and the spacing between them shortened. However, as mentioned above, it is difficult to reduce the size of a wideband antenna, making it also difficult to shorten the spacing between antennas.

The present invention aims to provide a data processing method, measurement system, and program that can improve spatial resolution in measurements while maintaining a constant number of antennas.

A first aspect of the present invention is
A data processing method for analyzing scattered waves of waves radiated to an object, comprising:
The wave is emitted to the object from a plurality of transmitting points p ₁ (x' ₁ , y' ₁ , z' ₁ ) arranged on the y-axis;
The scattered waves reflected with reflectivity f(x, y, z) at a reflection point (x, y, z) on the object are received as measurement values s(x' ₁ , x' 2 , y' ₁ , y' 2 , z' ₁ , z' ₂ , k) at a plurality of reception points p ₂ ( _x ' ₂ , y' ₂ , z' ₂ ) arranged on the _y axis,
The measurement value s( _x'1 , _x'2 , _y'1 , _y'2 , _z'1 , _z'2 , k) is triple-Fourier transformed using equation (1) to obtain S( _k'x1 , _k'x2 , _k'y1 , _k'y2 , _z'1 , _z'2 , k),

If the Nyquist wave number in the x direction determined by the measurement interval in the x direction is kx _,nyq , then variable substitution processing in the x direction is performed within the range of _-kx1 ≦ _kx ≦ _kx1 (where _kx,nyq ≦ _kx1 ),
If the Nyquist wave number in the _y'1 direction determined by the measurement interval in the y direction is k'y1 _,nyq and the Nyquist wave number in the _y'2 direction is _{k'y2 ,nyq} , then variable substitution processing in the y direction is performed within the range of _-ky1 ≦ _k'y1 ≦ _ky1 and -ky2 ≦ _k'y2 ≦ _ky2 (where k'y1 _,nyq ≦ _ky1 and _k'y2, nyq ≦ _ky2 ),
The reflectance f(x, y, z) is calculated by performing a triple inverse Fourier transform using equation (2).

It is a data processing method.
however,
x' = x' ₁ = x' ₂
z' ₁ = z' ₂ = 0
k is the wave number of the propagating wave,
k' _x1 , k' _y1 , k' _z1 are components of the wave vector of the spherical wave of the wave propagating between the transmitting point p ₁ (x' ₁ , y' ₁ , z' ₁ ) and the reflecting point (x, y, z),
k' _x2 , k' _y2 , k' _z2 are components of the wave vector of the spherical wave of the wave propagating between the reflection point (x, y, z) and the reception point 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
is.

A second aspect of the present invention is
A measurement system for analyzing scattered waves of waves radiated to an object, comprising:
A transceiver unit,
a transmitting unit that radiates the wave toward the object from a plurality of transmitting points p ₁ (x' ₁ , y' ₁ , z' ₁ ) arranged on the y-axis;
a receiving unit that receives the scattered waves reflected with reflectivity f(x, y, z) at a reflection point (x, y, z) on the object as measurement values s(x' ₁ , x' ₂ , y' 1 , y' ₂ , z' _{1 , z' 2 , k) at a plurality of receiving points p 2 (x' 2 , y ' 2} , z' ₂ ) arranged on the _y axis;
a transceiver unit having
A processing device,
a step of performing a triple Fourier transform on the measured values s( _x'1 , _x'2 , _y'1 , _y'2 , _z'1 , _z'2 , k) using equation (1) to obtain S( _k'x1 , _k'x2 , _k'y1 , _k'y2 , _z'1 , _z'2 , k);

When the Nyquist wave number in the x direction determined by the measurement interval in the x direction is kx _,nyq , a procedure for performing variable substitution processing in the x direction within the range of _-kx1 ≦ _kx ≦ _kx1 (where _kx,nyq ≦ _kx1 ) is performed.
where the Nyquist wave number in the _y'1 direction determined by the measurement interval in the y direction is k'y1 _,nyq and the Nyquist wave number in the _y'2 direction is _{k'y2 ,nyq} , a procedure for performing variable substitution processing in the _y direction within the range of _-ky1 ≦ _k'y1 ≦ _ky1 and -ky2 ≦k'y2 ≦ _ky2 (where k'y1 _,nyq ≦ _{ky1 and k'y2,} nyq ≦ _ky2 );
A procedure of calculating the reflectance f(x, y, z) by triple inverse Fourier transform using equation (2);

a processing unit for executing the
The measurement system has the following features:
however,
x' = x' ₁ = x' ₂
z' ₁ = z' ₂ = 0
k is the wave number of the propagating wave,
k' _x1 , k' _y1 , k' _z1 are components of the wave vector of the spherical wave of the wave propagating between the transmitting point p ₁ (x' ₁ , y' ₁ , z' ₁ ) and the reflecting point (x, y, z),
k' _x2 , k' _y2 , k' _z2 are components of the wave vector of the spherical wave of the wave propagating between the reflection point (x, y, z) and the reception point 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
is.

A third aspect of the present invention is
A program for analyzing scattered waves of waves radiated to an object,
The procedure involves triple Fourier transforming the measured values s( _x'1 , _x'2 , _y'1 , _y'2 , _z'1 , _z'2 , k) using equation (1) to obtain S( _k'x1 , _k'x2 , _k'y1 , _k'y2 , _z'1 , _z'2 , k), and

When the Nyquist wave number in the x direction determined by the measurement interval in the x direction is kx _,nyq , a procedure for performing variable substitution processing in the x direction within the range of _-kx1 ≦ _kx ≦ _kx1 (where _kx,nyq ≦ _kx1 ) is performed.
where the Nyquist wave number in the _y'1 direction determined by the measurement interval in the y direction is k'y1 _,nyq and the Nyquist wave number in the _y'2 direction is _{k'y2 ,nyq} , a procedure for performing variable substitution processing in the _y direction within the range of _-ky1 ≦ _k'y1 ≦ _ky1 and -ky2 ≦k'y2 ≦ _ky2 (where k'y1 _,nyq ≦ _{ky1 and k'y2,} nyq ≦ _ky2 );
A procedure for calculating the reflectance f(x, y, z) by triple inverse Fourier transform using equation (2),

It is a program that causes a computer to execute the above.
however,
x' = x' ₁ = x' ₂
z' ₁ = z' ₂ = 0
k is the wave number of the propagating wave,
k' _x1 , k' _y1 , k' _z1 are components of the wave vector of the spherical wave of the wave propagating between the transmitting point p ₁ (x' ₁ , y' ₁ , z' ₁ ) and the reflecting point (x, y, z),
k' _x2 , k' _y2 , k' _z2 are components of the wave vector of the spherical wave of the wave propagating between the reflection point (x, y, z) and the reception point 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
is.

A fourth aspect of the present invention is
A data processing method for analyzing scattered waves of waves radiated to an object, comprising:
The wave is emitted to the object from a plurality of transmitting points p ₁ (x' ₁ , y' ₁ , z' ₁ ) arranged two-dimensionally on an xy plane;
The scattered waves reflected with reflectivity f(x, y, z) at a reflection point (x, y, z) on the object are received as measurement values s(x' ₁ , x' 2 , y' ₁ , y' 2 , z' 1 , z' ₂ , k) at a plurality of reception points p ₂ (x' ₂ , y' ₂ , z _{' 2 )} arranged two-dimensionally on _{the xy} plane,
The measurement value s( _x'1 , _x'2 , _y'1 , _y'2 , _z'1 , _z'2 , k) is subjected to a quadruple Fourier transform using equation (1) to obtain S( _k'x1 , _k'x2 , _k'y1 , _k'y2 , _z'1 , _z'2 , k),

Let _k'x1,nyq be the Nyquist wave number in the _x'1 direction determined by the measurement interval in the x direction, and _k'x2,nyq be the Nyquist wave number in the _x'2 direction. Then, perform variable substitution processing in the x direction within the range of _-kx1 ≦ _k'x1 ≦ _kx1 and _-kx2 ≦ _k'x2 ≦ _kx2 (where k'x1 _,nyq ≦ _kx1 and _k'x2,nyq ≦ _kx2 ),
If the Nyquist wave number in the _y'1 direction determined by the measurement interval in the y direction is k'y1 _,nyq and the Nyquist wave number in the _y'2 direction is _{k'y2 ,nyq} , then variable substitution processing in the y direction is performed within the range of _-ky1 ≦ _k'y1 ≦ _ky1 and -ky2 ≦ _k'y2 ≦ _ky2 (where k'y1 _,nyq ≦ _ky1 and _k'y2, nyq ≦ _ky2 ),
The reflectance f(x, y, z) is calculated by performing a triple inverse Fourier transform using equation (2).

It is a data processing method.
however,
z' ₁ = z' ₂ = 0
k is the wave number of the propagating wave,
k' _x1 , k' _y1 , k' _z1 are components of the wave vector of the spherical wave of the wave propagating between the transmitting point p ₁ (x' ₁ , y' ₁ , z' ₁ ) and the reflecting point (x, y, z),
k' _x2 , k' _y2 , k' _z2 are components of the wave vector of the spherical wave of the wave propagating between the reflection point (x, y, z) and the reception point 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
is.

A fifth aspect of the present invention is
A measurement system for analyzing scattered waves of waves radiated to an object, comprising:
A transceiver unit,
a transmitting unit that radiates the wave toward the object from a plurality of transmitting points p ₁ (x' ₁ , y' ₁ , z' ₁ ) that are two-dimensionally arranged on an xy plane;
a receiving unit that receives the scattered waves reflected with reflectivity f(x, y, z) at a reflection point (x, y, z) on the object as measurement values s(x' ₁ , x' ₂ , y' 1 , y' ₂ , z' 1 , z' ₂ , k) at a plurality of receiving points p _{2 (} x' ₂ , y' ₂ , z' ₂ ) that are two-dimensionally arranged on _{the xy} plane;
a transceiver unit having
A processing device,
a step of performing a quadruple Fourier transform on the measured values s( _x'1 , _x'2 , _y'1 , _y'2 , _z'1 , _z'2 , k) using equation (1) to obtain S( _k'x1 , _k'x2 , _k'y1 , _k'y2 , _z'1 , _z'2 , k);

Let _k'x1,nyq be the Nyquist wave number in the _x'1 direction determined by the measurement interval in the x direction, and _k'x2,nyq be the Nyquist wave number in the _x'2 direction. Then, perform variable substitution processing in the x direction within the range of _-kx1 ≦ _k'x1 ≦ _kx1 and _-kx2 ≦ _k'x2 ≦ _kx2 (where k'x1 _,nyq ≦ _kx1 and _k'x2,nyq ≦ _kx2 );
where the Nyquist wave number in the _y'1 direction determined by the measurement interval in the y direction is k'y1 _,nyq and the Nyquist wave number in the _y'2 direction is _{k'y2 ,nyq} , a procedure for performing variable substitution processing in the _y direction within the range of _-ky1 ≦ _k'y1 ≦ _ky1 and -ky2 ≦k'y2 ≦ _ky2 (where k'y1 _,nyq ≦ _{ky1 and k'y2,} nyq ≦ _ky2 );
A procedure of calculating the reflectance f(x, y, z) by triple inverse Fourier transform using equation (2);

a processing unit for executing the
The measurement system has the following features:
however,
z' ₁ = z' ₂ = 0
k is the wave number of the propagating wave,
k' _x1 , k' _y1 , k' _z1 are components of the wave vector of the spherical wave of the wave propagating between the transmitting point p ₁ (x' ₁ , y' ₁ , z' ₁ ) and the reflecting point (x, y, z),
k' _x2 , k' _y2 , k' _z2 are components of the wave vector of the spherical wave of the wave propagating between the reflection point (x, y, z) and the reception point 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
is.

A sixth aspect of the present invention is
A program for analyzing scattered waves of waves radiated to an object,
The procedure involves performing a quadruple Fourier transform on the measured values s( _x'1 , _x'2 , _y'1 , _y'2 , _z'1 , _z'2 , k) using equation (1) to obtain _S ( _k'x1 , k'x2, _k'y1 , _k'y2 , _z'1 , _z'2 , k), and

Let _k'x1,nyq be the Nyquist wave number in the _x'1 direction determined by the measurement interval in the x direction, and _k'x2,nyq be the Nyquist wave number in the _x'2 direction. Then, perform variable substitution processing in the x direction within the range of _-kx1 ≦ _k'x1 ≦ _kx1 and _-kx2 ≦ _k'x2 ≦ _kx2 (where k'x1 _,nyq ≦ _kx1 and _k'x2,nyq ≦ _kx2 );
where the Nyquist wave number in the _y'1 direction determined by the measurement interval in the y direction is k'y1 _,nyq and the Nyquist wave number in the _y'2 direction is _{k'y2 ,nyq} , a procedure for performing variable substitution processing in the _y direction within the range of _-ky1 ≦ _k'y1 ≦ _ky1 and -ky2 ≦k'y2 ≦ _ky2 (where k'y1 _,nyq ≦ _{ky1 and k'y2,} nyq ≦ _ky2 );
A procedure for calculating the reflectance f(x, y, z) by triple inverse Fourier transform using equation (2),

It is a program that causes a computer to execute the above.
however,
z' ₁ = z' ₂ = 0
k is the wave number of the propagating wave,
k' _x1 , k' _y1 , k' _z1 are the components of the wave vector of the spherical wave of the wave propagating between the transmitting point p ₁ (x' ₁ , y' ₁ , z' ₁ ) and the reflecting point (x, y, z),
k' _x2 , k' _y2 , k' _z2 are components of the wave vector of the spherical wave of the wave propagating between the reflection point (x, y, z) and the reception point 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
is.

The data processing method, measurement system, and program of the present invention make it possible to improve the spatial resolution in measurements while maintaining the number of antennas constant.

FIG. 1 is a diagram showing a configuration of a radar device according to a first embodiment; FIG. 2 is a diagram showing the configuration of the array antenna shown in FIG. 1; FIG. 1 is a diagram illustrating the positional relationship between the array antenna and the object to be measured in the first embodiment. 1 is a flowchart showing a data processing method according to a first embodiment; Aliasing concept 1. An example of a measured waveform according to the first embodiment Fourier transform for the x-axis in the first embodiment Fourier transform for the y-axis in the first embodiment FIG. 10 is a diagram illustrating the positional relationship between the array antenna and the object to be measured in the second embodiment. 10 is a flowchart showing a data processing method according to a second embodiment. Fourier transform for _x'1 axis and _x'2 axis in the second embodiment Fourier transform for _y'1 axis and _y'2 axis in the second embodiment FIG. 10 is a diagram illustrating the positional relationship between the array antenna and the object to be measured according to the third embodiment. 10 is a flowchart showing a data processing method according to a third embodiment. Fourier transform for the x-axis in the third embodiment Fourier transform for _y'1 axis and _y'2 axis in the third embodiment Xz plane image of point target 1 simulated by the data processing method of Comparative Example 1 Xz plane image of point target 1 simulated by the data processing method of Example 1 YZ plane image of point target 1 simulated by the data processing method of Comparative Example 1 yz-plane image of point target 1 simulated by the data processing method of Example 1 Waveforms in the x-axis direction of point target 1 in Comparative Example 1 and Example 1 Waveforms in the y-axis direction of point target 1 in Comparative Example 1 and Example 1 Waveforms in the z-axis direction of point target 1 in Comparative Example 1 and Example 1 Xz plane image of point target 2 simulated by the data processing method of Comparative Example 1 Xz plane image of point target 2 simulated by the data processing method of Example 1 YZ plane image of point target 2 simulated by the data processing method of Comparative Example 1 yz-plane image of point target 2 simulated by the data processing method of Example 1 Waveforms in the x-axis direction of point target 2 in Comparative Example 1 and Example 1 Waveforms in the y-axis direction of point target 2 in Comparative Example 1 and Example 1 Waveforms in the z-axis direction of the point target 2 in Comparative Example 1 and Example 1 Conceptual diagram of a cross-shaped k' _y1 -k' _y2 filter Simulation results for four point targets Point target recovery rate Transmit and receive antenna layout 3D image of the test specimen without super-resolution processing 3D image of the test specimen after super-resolution processing 3D image of the test specimen after super-resolution processing using a cross-shaped k' _y1 -k' _y2 filter Waveforms in the x-axis direction of point target 1 in Comparative Example 2 and Example 2 Waveforms in the y-axis direction of point target 1 in Comparative Example 2 and Example 2 Waveforms in the z-axis direction of point target 1 in Comparative Example 2 and Example 2 Waveforms in the x-axis direction of point target 2 in Comparative Example 2 and Example 2 Waveforms in the y-axis direction of point target 2 in Comparative Example 2 and Example 2 Waveforms in the z-axis direction of point target 2 in Comparative Example 2 and Example 2 Conceptual diagram of cross-shaped k' _x1 -k' _x2 filter

First Embodiment
The data processing method, measurement system, and program of the first embodiment will be described in detail below. Fig. 1 shows the configuration of a radar device of this embodiment. Fig. 2 shows the configuration of the array antenna shown in Fig. 1. Fig. 3 is a diagram explaining the positional relationship between the array antenna of this embodiment and a measurement target. In this embodiment, electromagnetic waves are described as waves radiated into space, but waves that propagate through space, such as X-rays or ultrasound, may be used instead of electromagnetic waves.

The measurement system 1 of this embodiment has a transceiver unit and a processing device. The processing device may be provided integrally with the transceiver unit, or may be provided in a separate location connected to the transceiver unit via a network. In the following embodiment, an example will be described in which the processing device is provided integrally with the transceiver unit.

The radar device 60 of this embodiment shown in Figure 1 uses a transmitting array antenna and a receiving array antenna (transmitter/receiver) to emit electromagnetic waves from the transmitting antenna while sweeping the frequency of the electromagnetic waves. The radar device 60 then receives the reflected waves from the object to be measured using the receiving antenna to obtain measurement data s(x', y', z', k). The measurement data s(x', y', z', k) is data whose variables are the x-coordinate component, y-coordinate component, z-coordinate component, and the frequency of the electromagnetic waves.

The radar device 60 has a measurement unit 61, a data processing unit (processing device) 66, and an image display unit 68. The measurement unit 61 has a transmitting array antenna 50, a receiving array antenna 52, high-frequency switches 58 and 59, a high-frequency circuit 62, and a system control circuit 64. The radar device 60 emits electromagnetic waves of 10 MHz or higher, for example, 10 to 20 GHz, although the frequency of the electromagnetic waves is not particularly limited.

2, the transmitting array antenna 50 has a plurality of transmitting antennas 10a arranged in one direction. Each transmitting antenna 10a radiates electromagnetic waves toward the object to be measured. The receiving array antenna 52 has a plurality of receiving antennas 10b arranged in the same direction as the transmitting antennas 10a. Each receiving antenna 10b receives electromagnetic waves reflected from the object to be measured.
The transmitting antennas 10a of the transmitting array antenna 50 and the receiving antennas 10b of the receiving array antenna 52 are arranged on the same plane. The transmitting array antenna 50 and the receiving array antenna 52 are arranged so that the object to be measured faces this plane.

The data processing unit 66 processes multiple measurement data obtained by transmission to the measurement object by the multiple transmitting antennas 10a and reception by the multiple receiving antennas 10b, and calculates image data related to the measurement object. In this embodiment, the transmitting antennas 10a and receiving antennas 10b are planar antennas in which antenna patterns are formed in a plane on a substrate, but are not limited to planar antennas.

The transmitting array antenna 50 and the receiving array antenna 52 move parallel to the surface of the object to be measured. That is, the transmitting array antenna 50 and the receiving array antenna 52 perform measurements while scanning along the surface of the object to be measured. As the transmitting array antenna 50 and the receiving array antenna 52 move, the system control circuit 64 controls the operation of the high-frequency circuit 62. Specifically, the system control circuit 64 controls the operation of the high-frequency circuit 62 so that the transmitting antenna 10a is switched by the high-frequency switch 58 and electromagnetic waves are radiated for each unit length of movement of the transmitting array antenna 50 and the receiving array antenna 52.

The radar device 60 has an encoder 69. The encoder 69 generates a pulse signal every time the radar device 60 moves a certain distance. The encoder 69 detects the movement of the transmitting array antenna 50 and the receiving array antenna 52.
At this time, every time an electromagnetic wave is radiated from each of the transmitting antennas 10a, the high frequency switch 59 sequentially switches among the plurality of receiving antennas 10b, causing each receiving antenna 10b to receive the wave.

The frequency of the electromagnetic waves radiated from the transmitting array antenna 50 is swept at set frequency intervals within a range of, for example, 10 to 20 GHz for a fixed period of time, and the electromagnetic waves are radiated. Therefore, the measurement data obtained from the high-frequency circuit 62 is data whose value is determined by the position from which the transmitting antenna 10a transmits, the position from which the receiving antenna 10b receives, the frequency, and the position of the target.
At this time, the operation of the high-frequency switch 59 is controlled so that when an electromagnetic wave radiated from the transmitting antenna 10a is reflected by the object to be measured, the reflected wave of the electromagnetic wave is received by the receiving antenna 10b closest to the transmitting antenna 10a that radiated the electromagnetic wave. The receiving microwave amplifier (RF amplifier) may be set to change the gain for each pair of the transmitting antenna 10a that transmits and the receiving antenna 10b that receives the electromagnetic wave. In this case, the high-frequency circuit 62 has a variable gain amplification function that switches the gain depending on the selected pair of the transmitting antenna 10a and the receiving antenna 10b. This allows the depth that can be inspected for defects, etc. in the object to be measured to be increased.

In this embodiment, the arrangement directions of the transmitting antenna 10a and the receiving antenna 10b are parallel, and as shown in Figure 2, the arrangement direction is defined as the y direction. Meanwhile, the movement direction (scanning direction) of the transmitting array antenna 50 and the receiving array antenna 52 is defined as the x direction. When viewed from the transmitting array antenna 50 and the receiving array antenna 52, the direction in which the object to be measured is located (the direction in which electromagnetic waves are transmitted) is defined as the z direction.

The movement direction (scanning direction) of the transmitting array antenna 50 and the receiving array antenna 52 may be the y direction. That is, they may move (scan) in the same direction as the arrangement direction of the transmitting antennas 10a and the receiving antennas 10b.
Alternatively, the transmitting array antenna 50 may have only one transmitting antenna 10a, and the receiving array antenna 52 may have multiple receiving antennas 10b. In this case, the movement direction (scanning direction) of the transmitting array antenna 50 and the receiving array antenna 52 may also be the y direction. In other words, they may move (scan) in the same direction as the arrangement direction of the receiving antennas 10b.

The data processing unit 66 processes the measurement data s(x', y', z', k) obtained by transmitting and receiving electromagnetic waves using the transmitting array antenna 50 and the receiving array antenna 52 to create image data representing the interior of the object to be measured. The data processing unit 66 is configured, for example, by a computer, and invokes and starts a program stored in the memory unit 66a. This enables the functions of the data processing unit 66 to be fulfilled. In other words, the data processing unit 66 is configured as a software module. The image display unit 68 uses the created image data to display an image of the interior of the object to be measured.

2 shows a schematic diagram of a transmitting array antenna 50 and a receiving array antenna 52. The transmitting antenna 10a and the receiving antenna 10b are offset in the x direction by ΔL, but in the following description, the x direction positions of the transmitting antenna 10a and the receiving antenna 10b are assumed to be at a circular point midway between the transmitting antenna 10a and the receiving antenna 10b. This circular point is called the transmitting and receiving point.
Note that there may be no misalignment in the y direction between the transmitting antenna 10a and the receiving antenna 10b. That is, there may be cases where Δy = 0. There may also be cases where the transmitting antenna 10a and the receiving antenna 10b are shared. That is, there may be cases where Δy = 0 and ΔL = 0.

Therefore, the positional relationship between the object to be measured, the transmitting array antenna 50, and the receiving array antenna 52 can be expressed as shown in FIG.
Here, the coordinates of the transmitting and receiving points are assumed to be p(x', y', z'). The reflectance at the reflection point (x, y, z) of the object to be measured is assumed to be f(x, y, z). The measurement data at the transmitting and receiving point p(x', y', z') is assumed to be s(x', y', z', k). The propagation wavelength of the electromagnetic wave in a vacuum is assumed to be _λ0 . The relative dielectric constant of the medium is assumed to be _εr . The wave number of the propagating electromagnetic wave is assumed to be k.

In this case, the measurement data s(x', y', z', k) at the transmitting/receiving point p(x', y', z') can be expressed by the following equation.

however,

is.

In equation (1-1), electromagnetic waves are expressed as spherical waves, and distance attenuation is omitted. This distance attenuation is omitted because it has little effect on subsequent processing. When the exponent part of the integrand function in the second equation in equation (1-1) is expressed in Fourier transform notation, it becomes the following equation. This is equivalent to decomposing the round-trip spherical wave of equation (1-1) into three-dimensional plane waves.

Here, (k _x , k _y , k _z ) are components of the wave vector of the round-trip spherical wave of the wave propagating between the transmitting/receiving point p(x', y', z') and the reflecting point (x, y, z).

Meet the following.

Hereinafter, the reflectance f(x, y, z) is derived from the measurement data s(x', y', z', k) based on the formula (1-3). First, the formula (1-3) is rearranged as follows:

Here, the integral inside { } is a triple Fourier transform with respect to (x, y, z). Also, the integral inside [ ] is a double inverse Fourier transform with respect to (k _x , k _y ). Therefore, a double Fourier transform is performed on both sides of equation (1-5) with respect to (x', y'). Let F(k x , k y , k z ) be the function after the triple Fourier transform of function f( _x , _y , _z ). Let S(k x , k y , z', k) be the function after the double Fourier transform of measurement data s( _x' , _y' , z', k). In this case, equation (1-5) is expressed as follows:

When both sides of the second line of equation (1-6) are subjected to a triple inverse Fourier transform with respect to (k _x , k _y , k _z ), the reflectance f(x, y, z) is obtained as follows:

In this embodiment, as shown in FIG. 3, the transmitting and receiving point p(x', y', z') is arranged on the xy plane, and z'=0, so equation (1-7) can be expressed as follows:

As described above, the data processing unit 66 calculates the reflectance f(x, y, z) based on the measurement data s(x', y', z', k).

The data processing method and program of this embodiment will be described below with reference to Fig. 4. Fig. 4 is a flowchart showing the data processing method of this embodiment.
First, the measurement unit 61 acquires the measurement data s(x', y', 0, k) (step S1-1). Then, the data processing unit 66 performs a Hilbert transform on the measurement data s(x', y', 0, k) (step S1-2). This obtains the imaginary component of the frequency data at each transmitting and receiving point.

Next, the data processing unit 66 performs a double Fourier transform with respect to (x', y') on the measurement data s(x', y', 0, k) (step S1-3), thereby obtaining S(k _x , k _y , 0, k) as shown in equation (1-6).
Next, the data processing unit 66 performs variable substitution on S( _kx , _ky , 0, k) (step S1-4). Specifically, using equation (1-4), the function of ( _kx , _ky , k) is converted into a function of ( _kx , _ky , _kz ). This results in S( _kx , _ky , _kz ).
Next, the data processing unit 66 performs a triple inverse Fourier transform on S( _kx , _ky , _kz ) (step S1-5), thereby obtaining the reflectance f( _x , _y , _z ) as shown in equation (1-8).

The memory unit 66a stores a program for executing the data processing method of this embodiment. The program stored in the memory unit 66a causes the data processing unit 66 to execute the data processing method of this embodiment.

Here, we will explain alias data, which is considered to be aliasing noise in Fourier transform processing. When measuring with an array antenna, aliasing will occur in the measurement data unless the antenna size and scanning interval are below the Nyquist criterion. Aliasing occurs when continuous signals with different frequency components become indistinguishable due to sampling.

5 shows a conceptual diagram of aliasing. The Nyquist wave number determined from the spatial sampling interval is defined as k _nyq . The Nyquist wave number k _nyq is expressed by the following equation (1-10).

Here, Δx and Δy are the measurement intervals in the x-axis and y-axis directions, respectively.
Signals exceeding the Nyquist wave number k _nyq are called aliases and are treated as aliasing noise in signal processing.

6 shows an example of a measured waveform in the first embodiment. Here, the frequency of the electromagnetic wave is 15 GHz, the relative dielectric constant is 1, the measurement interval Δy = 7.5 mm, the measurement width is 562.5 mm, and the depth of the point target is 120 mm directly below the origin. Here, the wave number k y of the waveform measured at a position a distance y from the origin is expressed by the following equation (1-11).

Here, θ _ty is the angle between the line connecting the outermost transmitting/receiving antenna and the target and the z-axis. The Nyquist wave number is expressed by the following equation (1-12).

From this, the Nyquist criterion is expressed by the following equation (1-13).

Any waveform of the y coordinate that does not satisfy this condition, that is, data in the aliasing region in FIG. 6, will be alias data (aliasing noise).

From the above equation, we can see that the alias data increases as the wavelength of the electromagnetic wave becomes shorter, the measurement interval becomes longer, or the target becomes closer to the measurement surface. Normally, alias data is also treated as aliasing noise during signal processing. Therefore, an increase in alias data causes a deterioration in image resolution, a decrease in image intensity for shallow targets, and a deterioration in the S/N ratio.

In the following, it is assumed that the measurement interval Δx=10 mm, the measurement interval Δy=19.25 mm, and the wavelength of the maximum frequency in the medium λ _fmax =21.24 mm. In this case, the following equation (1-14) is obtained.

In this case, the Nyquist criterion is not satisfied, and aliasing occurs in both the x-axis and y-axis directions. The Nyquist wave number in this case is expressed by the following equation (1-15).

The maximum wave number obtained from the highest frequency is expressed by the following equation (1-16).

Here, the Fourier transform for the x-axis in step S1-3 is shown in Fig. 7. Also, the Fourier transform for the y-axis in step S1-3 is shown in Fig. 8.
As shown in Figure 7, in the Fourier transform for the x-axis, wavenumbers with ± _kx,nyq as their upper limit are output in range A. Data for wavenumbers exceeding ± _kx,nyq appear as alias data, with a period of _2kx,nyq . Therefore, in the variable substitution process of step S1-4, variable substitution is performed not in range A but in range A' with ± _kxmax as its upper limit. This allows the alias data to be added to the variable substitution process as actually meaningful data rather than discarded as aliasing noise, thereby improving the resolution in the x-axis direction and further increasing the image intensity of shallow targets.
Similarly, for the y-axis, variable substitution processing is performed in range B' instead of range B, as shown in FIG.

Second Embodiment
The data processing method, measurement system, and program of the second embodiment will be described in detail below. The transmitting array antenna 50 and receiving array antenna 52 of the first embodiment are arranged in one direction (the y direction in FIG. 3), but the arrangement of the transmitting array antenna 50 and receiving array antenna 52 of this embodiment is different. The transmitting array antenna 50 and receiving array antenna 52 of this embodiment are arranged in a plane.
In the first embodiment, the coordinates of the transmission points and reception points are both p(x', y', z'), but in this embodiment, the coordinates of the transmission points and reception points are different. In this embodiment, as shown in Fig. 9, a transmission point _p1 ( _x'1 , _y'1 , _z'1 ) and a reception point _p2 ( _x'2 , _y'2 , _z'2 ) are arranged on the xy plane.

Here, let f(x, y, z) be the reflectivity at the reflection point (x, y, z) of the object to be measured. Let s(x' ₁ , x' ₂ , y' ₁ , z' ₁ ) and s(x' 1 , x' ₂ , y' 1 , y' ₂ , z' _{1 , z' 2} , k) be the measurement data at the transmission point p ₁ (x' ₁ , y' ₁ , z' ₁ ) and the reception point p 2 (x' 2 , y' ₂ , z _{' 2 )} . Let λ ₀ be the propagation wavelength of the electromagnetic wave in a vacuum. Let ε _r be the relative dielectric constant of the medium. Let k be the wave number of the propagating electromagnetic wave.

In this case, the measurement data s(x' ₁ , x' ₂ , y' ₁ , y' ₂ , z' ₁ , z' ₂ , k) can be expressed by the following equation.

however,

is.

In equation (2-1), electromagnetic waves are represented as spherical waves, and distance attenuation is omitted. This distance attenuation is omitted because it has little effect on subsequent processing. When the exponential part of the integrand in equation (2-1) is expressed in Fourier transform notation, it becomes the following equation. This is equivalent to decomposing the spherical wave in equation (2-1) into three-dimensional plane waves.

Here, ( _k'x1 , _k'y1 , _k'z1 ) are the components of the wave vector of the spherical wave of the wave propagating from the transmitting point to the reflecting point, and ( _k'x2 , _k'y2 , _k'z2 ) are the components of the wave vector of the spherical wave of the wave propagating from the reflecting point to the receiving point.

Meet the following.

Below, the reflectance f(x, y, z) is derived from s( _x'1 , _x'2 , _y'1 , _y'2 , _z'1 , _z'2 , k) based on equation (2-3). First, a quadruple Fourier transform is performed on both sides of equation (2-3) with respect to _x'1 , _x'2 , _y'1 , _y'2 .

The left side of equation (2-5) is rewritten and rearranged as equation (2-6) below.

Then, equation (2-5) is expressed as equation (2-7).

The following integral is performed on both sides of equation (2-7).

Here, the following variable substitutions are made:

Here, the following equation is obtained from equations (2-9) and (2-10).

The absolute value |J| of the Jacobian in this variable substitution is given by the following equations from equations (2-12) and (2-13).

Substituting equations (2-9), (2-10), (2-14), and (2-15) into equation (2-8) and performing variable substitution, the following equation is obtained.

Here, the integral with respect to (u, v) on the right-hand side of the second line of equation (2-16) is a constant and has been omitted. By performing a triple inverse Fourier transform on both sides of equation (2-16) with respect to (k _x , k _y , k _z ), the reflectance f(x, y, z) is obtained as follows:

To solve equation (2-17), it is necessary to express ( _k'z1 , _k'z2 , k) as ( _k'x1 , _k'x2 , _k'y1 , _k'y2 , _kz ) or ( _kx , u, _ky , v, _kz ).
By rearranging using equations (2-4) and (2-11), and if the transmitting point and receiving point are both located on the xy plane passing through the origin, z' ₁ =z' ₂ =0, so equation (2-17) can be expressed as follows:

As described above, the data processing unit 66 obtains the reflectance f(x, y, z) based on the measurement data s(x' ₁ , x' ₂ , y' ₁ , y' ₂ , z' ₁ , z' ₂ , k).

The data processing method and program of this embodiment will be described below with reference to Fig. 10. Fig. 10 is a flowchart showing the data processing method of this embodiment.
First, the measurement unit 61 acquires the measurement data s ( _x'1 , _x'2 , _y'1 , _y'2 , _z'1 , _z'2 , k) (step S2-1). Then, the data processing unit 66 performs a Hilbert transform on the measurement data s ( _x'1 , _x'2 , _y'1 , _y'2 , _z'1 , _z'2 , k) (step S2-2). This obtains the imaginary component of the frequency data at each measurement point.

Next, the data processing unit 66 performs a quadruple Fourier transform of the measurement data s( _x'1 , _x'2 , _{y'1, y'2, z'1 , z'2 ,} k) with respect to ( _x'1 , _x'2 , _y'1 , _y'2 ) (step S2-3), thereby obtaining S( _k'x1 , k'x2, _k'y1 , _k'y2 , ₀ , 0, k) as shown in equation (2-6).
Next, the data processing unit 66 performs variable substitution on S( _k'x1 , _k'x2 , _k'y1 , _k'y2 , 0,0,k) (step S2-4). Specifically, using equations (2-9), (2-10), (2-14), and (2-15), the function of ( _k'x1 , _k'x2 , _k'y1 , _k'y2 , k) is converted into a function of ( _kx , _ky , _kz ). This results in S( _kx , _ky , _kz ).
Next, the data processing unit 66 performs a triple inverse Fourier transform on S( _kx , _ky , _kz ) (step S2-5), thereby obtaining the reflectance f( _x , _y , _z ) as shown in equation (2-18).

The memory unit 66a stores a program for executing the data processing method of this embodiment. The program stored in the memory unit 66a causes the data processing unit 66 to execute the data processing method of this embodiment.

Here, the frequency of the electromagnetic wave is 4.46484 GHz, the relative dielectric constant is 10, the measurement interval Δx' ₁ = 38.5 mm, the measurement interval Δx' ₂ = 38.5 mm, the measurement interval Δy' ₁ = 38.5 mm, the measurement interval Δy' ₂ = 38.5 mm, the measurement width is 616 mm, and the wavelength of the maximum frequency in the medium λ _fmax = 21.24 mm. In this case, the following equation (2-20) is obtained.

In this case, the Nyquist criterion is not satisfied, and aliasing occurs in each of the x'1 _- axis direction, _x'2- axis direction, _y'1- axis direction, and _y'2- axis direction. The Nyquist wave number in this case is expressed by the following equation (2-21).

The maximum wave number obtained from the highest frequency is expressed by the following equation (2-22).

Here, the Fourier transforms for the _x'1 axis and _x'2 axis in step S2-3 are shown in Fig. 11. Also, the Fourier transforms for the _y'1 axis and _y'2 axis in step S2-3 are shown in Fig. 12.
11, in the Fourier transform for the _x'1 axis and the _x'2 axis, wavenumbers with upper limits of ± _k'x1,nyq and ±k'x2 _,nyq are output in range A. Data with wavenumbers exceeding ± _k'x1,nyq and ±k'x2 _,nyq appear as alias data, repeating in cycles of _2k'x1,nyq and 2k'x2 _,nyq , respectively. Therefore, in the variable substitution process of step S2-4, variable substitution is performed not in range A but in range A' with upper limits of ± _k'x1max and ± _k'x2max , respectively. In this way, by adding the alias data to the variable substitution process as actually meaningful data rather than discarding it as aliasing noise, it is possible to improve the resolution in the _x'1 axis and _x'2 axis directions and further increase the image intensity of shallow targets.
Similarly, for the _y'1 axis and _y'2 axis, as shown in FIG. 12, variable substitution processing is performed not within range B but within range B' having upper limits of ±k' _y1max and ±k' _y2max , respectively.

Third Embodiment
The data processing method, measurement system, and program of the third embodiment will be described in detail below. In the second embodiment, the transmitting array antenna 50 and the receiving array antenna 52 are arranged in a plane, but in this embodiment, the arrangement of the transmitting array antenna 50 and the receiving array antenna 52 is different. In this embodiment, the transmitting array antenna 50 and the receiving array antenna 52 are arranged in a straight line.
Specifically, in this embodiment, the transmitting antennas 10a and receiving antennas 10b are arranged in the y direction as shown in Fig. 13. The movement direction (scanning direction) of the transmitting array antenna 50 and receiving array antenna 52 is defined as the x direction. The direction toward the object to be measured (the transmission direction of electromagnetic waves) as viewed from the transmitting array antenna 50 and receiving array antenna 52 is defined as the z direction.
The movement direction (scanning direction) of the transmitting array antenna 50 and the receiving array antenna 52 may be the y direction.

Therefore, the positional relationship between the object to be measured, the transmitting array antenna 50, and the receiving array antenna 52 can be expressed as shown in FIG.
Here, let the coordinates of the transmitting point be p ₁ (x' ₁ , y' ₁ , z' ₁ ) and the coordinates of the receiving point be p ₂ (x' ₂ , y' ₂ , z' ₂ ). Let the reflectance at the reflection point (x, y, z) of the object to be measured be f(x, y, z). Let the measurement data at _{p 2} (x' ₂ , y' ₂ , z' ₂ ) be s(x' ₁ , x' ₂ , y' ₁ , y' ₂ , z' ₁ , z' ₂ , k). Let the propagation wavelength of the electromagnetic wave in a vacuum be λ _0. Let the relative dielectric constant of the medium be ε _r . Let the wave number of the propagating electromagnetic wave be k.

In this case, the measurement data s(x' ₁ , x' ₂ , y' ₁ , y' ₂ , z' ₁ , z' ₂ , k) can be expressed by the following equation.

however

is.

In equation (3-1), electromagnetic waves are represented as spherical waves, and distance attenuation is omitted. This distance attenuation is omitted because it has little effect on subsequent processing. If the exponential part of the integrand function in the second equation in equation (3-1) is expressed in Fourier transform notation, it becomes the following equation. This is equivalent to decomposing the spherical wave in equation (3-1) into three-dimensional plane waves.

Here, ( _k'x1 , _k'y1 , _k'z1 ) are the components of the wave vector of the spherical wave of the wave propagating from the transmitting point to the reflecting point, and ( _k'x2 , _k'y2 , _k'z2 ) are the components of the wave vector of the spherical wave of the wave propagating from the reflecting point to the receiving point.

Meet the following.

Here, since the x coordinates of transmission point p ₁ (x' ₁ , y' ₁ , z' ₁ ) and reception point p ₂ (x' ₂ , y' ₂ , z' ₂ ) are equal, if we set x' = x' ₁ = x' ₂ , equation (3-3) can be expressed as follows:

Here, the following variable substitutions are made:

The following equation is obtained from equation (3-6):

Calculating the absolute value of the Jacobian from equation (3-7) gives the following equation:

Substituting equations (3-6) and (3-8) into equation (3-5) and performing variable substitution, the following equation is obtained.

Here, the integral with respect to u in the second line of equation (3-9) is omitted because it is a constant. When a triple inverse Fourier transform is performed on both sides of equation (3-9) with respect to (x, y ₁ , y ₂ ), the following equation is obtained.

The left side of equation (3-10) is rewritten and rearranged as follows:

Then, equation (3-10) is expressed as follows:

By performing the following integral on both sides of equation (3-12), the following equation is obtained:

Here, the following variable substitutions are defined for (k' _y1 , k' _y2 ) and (k' _z1 , k' _z2 ).

Here, the following equation is obtained from equation (3-14):

The absolute value of the Jacobian under this variable substitution is given by the following formula:

Substituting equations (3-1), (3-15), and (3-17) into equation (3-13) and performing variable substitution, the following equation is obtained.

Here, the integral with respect to v on the right-hand side of the second line of equation (3-18) is omitted because it is a constant. When a triple inverse Fourier transform is performed on both sides of equation (3-18) with respect to (k _x , k _y , k _z ), the reflectance f(x, y, z) is obtained as follows:

Here, if z'=0 is set in order to align the measurement surface with the x'-y' plane passing through the origin, equation (3-19) can be expressed as follows:

To solve equation (3-20), k must be expressed as (k' _y1 , k' _y2 , k _z ) or (k _y , v, k _z ).
The four simultaneous equations, Equation (3-4), Equation (3-6), Equation (3-15), and Equation (3-21) obtained from the assumption, are solved.

From this, k is expressed by the following formula.

As described above, the data processing unit 66 obtains the reflectance f(x, y, z) based on the measurement data s(x' ₁ , x' ₂ , y' ₁ , y' ₂ , z' ₁ , z' ₂ , k).

The data processing method and program of this embodiment will be described below with reference to Fig. 14. Fig. 14 is a flowchart showing the data processing method of this embodiment.
First, the measurement unit 61 acquires the measurement data s(x', _y'1 , _y'2 , 0, 0, k) (step S3-1). Then, the data processing unit 66 performs a Hilbert transform on the measurement data s(x', _y'1 , _y'2 , 0, 0, k) (step S3-2). This obtains the imaginary component of the frequency data at each measurement point.

Next, the data processing unit 66 performs a triple Fourier transform with respect to (x', y' ₁ , y' ₂ ) on the measurement data s(x', y' ₁ , y' ₂ , 0, 0, k) (step S3-3), thereby obtaining S(k _x , k' _y1 , k' _y2 , 0, 0, k) as shown in equation (3-11).
Next, the data processing unit 66 performs variable substitution on S( _kx , _k'y1 , _k'y2 , 0,0,k) (step S3-4). Specifically, using equations (3-14) and (3-15), the function of ( _kx , _k'y1 , _k'y2 , k) is converted into a function of ( _kx , _ky , v, k). This results in S( _kx , _ky , v, 0,0,k).
Next, the data processing unit 66 performs a triple inverse Fourier transform on S( _kx , _ky , v, 0, 0, k) to ( _kx , _ky , _kz ) (step S3-5), thereby obtaining the reflectance f(x, y, z) as shown in equation (3-20).

The memory unit 66a stores a program for executing the data processing method of this embodiment. The program stored in the memory unit 66a causes the data processing unit 66 to execute the data processing method of this embodiment.

Here, the frequency of the electromagnetic wave is 4.46484 GHz, the relative dielectric constant is 10, the measurement interval Δx = 10 mm, the measurement interval Δy' ₁ = 38.5 mm, the measurement interval Δy' ₂ = 38.5 mm, the measurement width is 616 mm, and the wavelength of the maximum frequency in the medium λ _fmax = 21.24 mm. In this case, the following equation (3-23) is obtained.

In this case, the Nyquist criterion is not satisfied, and aliasing occurs in each of the x-axis direction, _y'1- axis direction, and _y'2- axis direction. The Nyquist wave number in this case is expressed by the following equation (3-24).

The maximum wave number obtained from the highest frequency is expressed by the following equation (3-25).

Here, the Fourier transform for the x-axis in step S3-3 is shown in Fig. 15. Moreover, the Fourier transform for the _y'1 axis and _y'2 axis in step S3-3 is shown in Fig. 16.
15, in the Fourier transform for the x-axis, wavenumbers with upper limits of ± _kx,nyq and ±k'x2 _,nyq are output in range A. Data with wavenumbers exceeding ± _kx, nyq appear as alias data, repeating with a period of _2kx,nyq . Therefore, in the variable substitution process of step S3-4, variable substitution is performed not in range A but in range A' with upper limits of ± _kxmax . In this way, by adding the alias data to the variable substitution process as actually meaningful data rather than discarding it as aliasing noise, it becomes possible to improve the resolution in the x-axis direction and further increase the image intensity of shallow targets.
As for the _y'1 axis and _y'2 axis, similarly to the second embodiment, as shown in FIG. 16, variable substitution processing is performed not within range B but within range B' having upper limits of ±k' _y1max and ±k' _y2max , respectively.

(Simulation results)
The results of a computer simulation of the data processing method of the third embodiment will be described below. The simulation conditions are as follows:

・Frequency band used f _min to f _max : DC to 4.46484 GHz
・Center frequency fc: 2.23242GHz
Wavelength λ _0c of center frequency in vacuum: 134.38 mm
Relative dielectric constant ε _r of the medium: 10
Wavelength λc of the center frequency in the medium: 42.50 mm
・Wavelength λ _fmax of the highest frequency in the medium: 21.24 mm
・Number of measurements Nx: 128pt
Measurement interval Δx: 10 mm
Measurement width: 1280 mm
Number of measurements N _y'1 : 16 pt
Measurement interval Δy'1: 38.5 mm
Measurement width: 616 mm
Number of measurements N _y'2 : 16 pt
Measurement interval Δy'2: 38.5 mm
Measurement width: 616 mm
Coordinates of point target 1 (unit: mm): (640, 307, 20)
Coordinates of point target 2 (unit: mm): (640, 307, 200)

Fig. 17 shows an xz plane image of the point target 1 simulated by the data processing method of Comparative Example 1. Fig. 18 shows an xz plane image of the point target 1 simulated by the data processing method of Example 1. Fig. 19 shows a yz plane image of the point target 1 simulated by the data processing method of Comparative Example 1. Fig. 20 shows a yz plane image of the point target 1 simulated by the data processing method of Example 1.
Fig. 21 shows waveforms in the x-axis direction of the point target 1 of Comparative Example 1 and Example 1. Fig. 22 shows waveforms in the y-axis direction of the point target 1 of Comparative Example 1 and Example 1. Fig. 23 shows waveforms in the z-axis direction of the point target 1 of Comparative Example 1 and Example 1. In all of Figs. 21 to 23, Example 1 is shown by a solid line, and Comparative Example 1 is shown by a dashed line.

Fig. 24 shows an xz plane image of the point target 2 simulated by the data processing method of Comparative Example 1. Fig. 25 shows an xz plane image of the point target 2 simulated by the data processing method of Example 1. Fig. 26 shows a yz plane image of the point target 2 simulated by the data processing method of Comparative Example 1. Fig. 27 shows a yz plane image of the point target 2 simulated by the data processing method of Example 1.
Fig. 28 shows waveforms in the x-axis direction of the point target 2 of Comparative Example 1 and Example 1. Fig. 29 shows waveforms in the y-axis direction of the point target 2 of Comparative Example 1 and Example 1. Fig. 30 shows waveforms in the z-axis direction of the point target 2 of Comparative Example 1 and Example 1. In all of Figs. 28 to 30, Example 1 is shown by a solid line, and Comparative Example 1 is shown by a dashed line.

In Comparative Example 1, in the variable substitution process of step S3-4 in this embodiment, the variable substitution process was performed in range A. In Example 1, in the variable substitution process of step S3-4 in this embodiment, the variable substitution process was performed in range A' with ±k _xmax as the upper limit, rather than range A. Hereinafter, the variable substitution process in Example 1 in range A' with ±k _xmax as the upper limit will be referred to as super-resolution processing.

From Figures 17 to 30, it was confirmed that, compared to Comparison Example 1, Example 1, which performs super-resolution processing, improves the resolution in the x-axis and y-axis directions for both point targets 1 and 2.

Fourth Embodiment
In this embodiment, noise that may occur in the third embodiment is reduced. Measurement data s in the y-axis direction is expressed by the following equation (4-1), where t is the amplitude of the electromagnetic wave propagating from the transmitting antenna to the target, r is the reflectivity of the target, and r is the amplitude of the electromagnetic wave propagating from the target to the receiving antenna.

In this case, if the amplitude components below the Nyquist frequency are t ₀ and r ₀ and the n-times aliased amplitude components (alias data) are t _n and r _n , t and r are respectively expressed by the following equation (4-2).

Substituting equation (4-2) into equation (4-1), the following equation (4-3) is obtained.

In the super-resolution processing shown in the third embodiment, for example, during variable substitution in step S3-4, in the ( _k'y1 , _k'y2 ) space shown in Figure 16, if _k'y1 caused by switching the transmitting antenna and _k'y2 caused by the receiving antenna are both within a range equal to or less than the Nyquist wave number, as shown in the following equation (4-4), all components including alias data (aliasing noise components) become noise components in the signal processing.

In this case, equation (4-3) is expressed as equation (4-5).

Here, r ₀ t ₀ surrounded by a solid line is a signal component, and the part surrounded by a dashed line is a noise component.

Next, when only k' _y1 falls within the aliasing noise region once during variable substitution, the condition of equation (4-6) is met, and equation (4-3) is expressed as equation (4-7).


Here, r ₀ t ₁ surrounded by a solid line is a signal component, and the part surrounded by a dashed line is a noise component.

Next, when k' _y2 falls into the double aliasing noise region during variable substitution, the condition of equation (4-8) is met, and equation (4-3) is expressed as equation (4-9).


Here, r ₂ t ₁ surrounded by a solid line is a signal component, and the part surrounded by a dashed line is a noise component.

As shown in Figure 5, when the _y'1 coordinate of the transmitting antenna measurement point and the _y'2 coordinate of the receiving antenna measurement point are close to the target, i.e., when they are located directly above the target, the spatial frequencies of the measurement data in _{the y'1} and _y'2 axis directions are low. Conversely, as the _y'1 and _y'2 coordinates move farther from the target, the spatial frequencies of the measurement data in the _y'1 and _y'2 axis directions are higher. When the wave number at that point exceeds the Nyquist wave number, it is included in the measurement data as an aliasing noise component. In this way, alias data components (aliasing noise components) tend to occur at points relatively far from the target. Therefore, due to distance attenuation, their amplitudes are smaller than those of data components below the Nyquist wave number. Therefore, the following equation (4-10) generally holds true.

As a result, the magnitude of the second-order terms related to r _n and t _n becomes smaller than the terms including r ₀ or t _0. This is expressed by the following equation (4-11).

Furthermore, unlike air, where the attenuation rate is small, when electromagnetic waves propagate through concrete or underground, where the attenuation rate is large, the magnitudes of r _n and t _n become smaller due to the attenuation rate caused by the medium in addition to the attenuation due to distance. In this case, equation (4-11) can be expressed as the following equation (4-12).

Therefore, the quadratic terms related to t _n and r _n in equations (4-5), (4-7), (4-9), and (4-12) result in data with a relatively poor S/N ratio compared to the other terms. Therefore, in order to improve the S/N ratio, the quadratic terms related to r _n and t _n can be selectively excluded in the variable substitution process of step S3-4. In this case, equation (4-3) can be expressed as the following equation (4-13).

Here, only the components enclosed by solid lines are used in the variable substitution process, and the parts enclosed by dashed lines are not used in the variable substitution process.

In this way, when electromagnetic waves propagate through a medium with a large attenuation rate, the degradation of the S/N ratio of a 3D image can be suppressed by using the filter expressed by equation (4-13) during variable substitution processing. This is expressed as a "cross-shaped _k'y1 - _k'y2 filter" in the ( _k'y1 , _k'y2 ) space. Figure 31 is a conceptual diagram of the cross-shaped _k'y1 - _k'y2 filter.

15 and ₁₆ depend on the position of the target or the directivity of the antenna. If targets in the vicinity of the measurement plane (z=0) are included in the measurement target, the maximum values may be determined by _{the directivities θ bx and θ by} of the antenna, as shown in equation (4-14).

however,

is.

(Simulation result 1)
The results of a computer simulation of the data processing method of the fourth embodiment will be described below. First, the effect of improving the image intensity of a shallow depth target will be verified under the same simulation conditions as in the third embodiment.
Super-resolution processing uses the alias data that was lost. As a result, the target image intensity (image amplitude) increases compared to when super-resolution processing is not performed. This effect is greater for shallow targets with a large amount of alias data than for deep targets. In other words, the image intensity of shallow targets, which has been reduced by aliasing, is largely restored by super-resolution processing.

Figure 32 shows the simulation results for four point targets. The horizontal position of each point target is the center of the measurement surface, and the depths are 20 mm, 200 mm, 400 mm, and 600 mm. In Figure 32, Example 1, in which super-resolution processing was performed, is shown by a solid line, and Comparative Example 1, in which super-resolution processing was not performed, is shown by a dashed line. The magnitude of the 3D image amplitude is normalized such that the amplitude of the 20 mm deep target without super-resolution processing is 1.

The shallower the target depth, the stronger the target image amplitude due to the effect of super-resolution processing. Because distance attenuation is taken into account in the simulation data, the amplitude of deep targets becomes very small, making it difficult to directly compare the recovery effects. Therefore, for the amplitude of targets at the same depth, the recovery rate is defined as the ratio of the 3D image amplitude with super-resolution processing to the 3D image amplitude without super-resolution processing. Figure 33 shows the recovery rate for each point target. Figure 33 shows that the shallower the target, the greater the recovery rate due to super-resolution processing, as there is more alias data. At a depth of 20 mm, the recovery rate reaches 6 times. On the other hand, the deeper the target, the less alias data there is and the smaller the recovery rate. At a depth of 600 mm, the recovery rate is approximately 1.2 times.

From the above, it was confirmed that using alias data through super-resolution processing improves image amplitude intensity. Significant improvements are seen, particularly in images of shallow targets with a large amount of alias data. This is primarily reflected in 3D images as an improvement in the S/N ratio of shallow-depth target images.

(Simulation result 2)
The results of computer simulations of the data processing methods of the third and fourth embodiments will be described below. The simulation conditions are as follows:
・Frequency band used f _min to f _max : DC to 4.46484 GHz
・Center frequency fc: 2.23242GHz
Wavelength λ _0c of center frequency in vacuum: 134.38 mm
Relative dielectric constant ε _r of the medium: 5
Wavelength λc of the center frequency in the medium: 60.1 mm
・Wavelength λ _fmax of the highest frequency in the medium: 30.05 mm
Measurement interval in scanning direction Δx: 10 mm
・Number of measurements Nx: 115pt
Measurement width: 1150 mm
Measurement interval of transmitting antenna Δy1: 38.5 mm
Number of transmitting antennas N _y1 : 12
Measurement interval of receiving antenna Δy2: 38.5 mm
Number of receiving antennas N _y2 : 12pt
Measurement width: 481.25 mm

The results of measurements on a reinforced concrete specimen are shown below. The layout of the transmitting and receiving antennas used for the measurements is shown in Figure 34. The specimen contains multiple rebars with different z-coordinates. Figure 35 is a 3D image of the specimen without super-resolution processing. Figure 36 is a 3D image of the specimen with super-resolution processing. The z-coordinates of multiple rebars are shown in Figures 35 and 36. Comparing Figure 35 and Figure 36, the image of the rebars in Figure 36, which underwent super-resolution processing, is thinner overall, confirming improved image resolution. In addition, the image intensity of the rebars at depths of 35 mm and 45 mm is stronger, confirming the recovery effect of the image amplitude of shallow depth targets. However, it was also confirmed that new unwanted noise (aliasing noise) such as interference fringes was generated between the rebars.

Next, the cross-shaped _k'y1 - _k'y2 filter described in the fourth embodiment and FIG. 31 was applied. FIG _. 37 is a 3D image of a test specimen subjected to super-resolution processing using a cross-shaped _k'y1 - _k'y2 filter. Comparing FIG _. 36, which is a 3D image of a test specimen subjected to super-resolution processing without applying a cross-shaped _k'y1 - _k'y2 filter, with FIG. 37, which is a 3D image of a test specimen subjected to super-resolution processing using a cross-shaped k'y1-k'y2 filter, it was confirmed that the cross-shaped _k'y1 - _k'y2 filter can reduce interference fringe noise (aliasing noise) while suppressing degradation in the resolution of rebar images and the intensity of shallow-depth rebar images.

Next, the influence of the cross-shaped k' _y1 -k' _y2 filter on the resolution was evaluated using simulation data, under the same conditions as in the third embodiment.
FIG. 38 shows waveforms in the x-axis direction of point target 1 in Comparative Example 2 and Example 2. FIG. 39 shows waveforms in the y-axis direction of point target 1 in Comparative Example 2 and Example 2. FIG. 40 shows waveforms in the z-axis direction of point target 1 in Comparative Example 2 and Example 2. FIG. 41 shows waveforms in the x-axis direction of point target 2 in Comparative Example 2 and Example 2. FIG. 42 shows waveforms in the y-axis direction of point target 2 in Comparative Example 2 and Example 2. FIG. 43 shows waveforms in the z-axis direction of point target 2 in Comparative Example 2 and Example 2.
In all of FIGS. 38 to 43, Example 2 is indicated by a solid line, and Comparative Example 2 is indicated by a dashed line.

In Comparative Example 2, the variable substitution process in step S3-4 in the third embodiment was performed in a range A' with ±k _xmax as the upper limit without applying the cross-shaped k' _y1 -k' _y2 filter. In Example 2, the variable substitution process in step S3-4 in the third embodiment was performed in a range A' with ±k _xmax as the upper limit with applying the cross-shaped k' _y1 -k' _y2 filter.
38 to 43, it was confirmed that even when super-resolution processing was performed using the cross-shaped k' _y1 -k' _y2 filter, there was almost no effect on the image resolution in any of the x-axis, y-axis, and z-axis.

Fifth Embodiment
In this embodiment, noise that may occur in the second embodiment is reduced. In the fourth embodiment, super-resolution processing is performed by applying a cross-shaped _k'y1 - _k'y2 filter in the third embodiment. In this embodiment, super-resolution processing is performed by applying a cross-shaped _k'x1 - _k'x2 filter and a cross-shaped _k'y1 - _k'y2 filter in the second embodiment. FIG. 44 is a conceptual diagram of the cross-shaped _k'x1 - _k'x2 filter. That is, in this embodiment, super-resolution processing is performed by applying both the cross-shaped _k'x1 - _k'x2 filter shown in FIG. 44 and the cross-shaped _k'y1 - _k'y2 filter shown in FIG. 31 in the second embodiment.

k' _x1max , k' _x2max , k' _y1max , and k' _y2max depend on the position of the target or the directivity of the antenna. If targets in the vicinity of the measurement plane (z=0) are included in the measurement targets, the maximum values may be determined by the directivities θ _bx and θ _by of the antenna, as shown in equation (5-1).

however,

is.

10a: transmitting antenna 10b: receiving antenna 50: transmitting array antenna 52: receiving array antenna 60: radar device 61: measuring unit 64: system control circuit 66: data processing unit 68: image display units 58 and 59: high frequency switch 62: high frequency circuit 69: encoder

Claims (13)

A data processing method for analyzing scattered waves of waves radiated to an object, comprising:
The wave is emitted to the object from a plurality of transmitting points p ₁ (x' ₁ , y' ₁ , z' ₁ ) arranged two-dimensionally on an xy plane;
The scattered waves reflected with reflectivity f(x, y, z) at a reflection point (x, y, z) on the object are received as measurement values s(x' ₁ , x' 2 , y' ₁ , y' 2 , z' 1 , z' ₂ , k) at a plurality of reception points p ₂ (x' ₂ , y' ₂ , z _{' 2 )} arranged two-dimensionally on _{the xy} plane,
The measurement value s( _x'1 , _x'2 , _y'1 , _y'2 , _z'1 , _z'2 , k) is subjected to a quadruple Fourier transform using equation (1) to obtain S( _k'x1 , _k'x2 , _k'y1 , _k'y2 , _z'1 , _z'2 , k),

Let _k'x1,nyq be the Nyquist wave number in the _x'1 direction determined by the measurement interval in the x direction, and _k'x2,nyq be the Nyquist wave number in the _x'2 direction. Then, perform variable substitution processing in the x direction within the range of _-kx1 ≦ _k'x1 ≦ _kx1 and _-kx2 ≦ _k'x2 ≦ _kx2 (where k'x1 _,nyq ≦ _kx1 and _k'x2,nyq ≦ _kx2 ),
If the Nyquist wave number in the _y'1 direction determined by the measurement interval in the y direction is k'y1 _,nyq and the Nyquist wave number in the _y'2 direction is _{k'y2 ,nyq} , then variable substitution processing in the y direction is performed within the range of _-ky1 ≦ _k'y1 ≦ _ky1 and -ky2 ≦ _k'y2 ≦ _ky2 (where k'y1 _,nyq ≦ _ky1 and _k'y2, nyq ≦ _ky2 ),
The reflectance f(x, y, z) is calculated by performing a triple inverse Fourier transform using equation (2).

Data processing methods.
however,
z' ₁ = z' ₂ = 0
k is the wave number of the propagating wave,
k' _x1 , k' _y1 , k' _z1 are components of the wave vector of the spherical wave of the wave propagating between the transmitting point p ₁ (x' ₁ , y' ₁ , z' ₁ ) and the reflecting point (x, y, z),
k' _x2 , k' _y2 , k' _z2 are components of the wave vector of the spherical wave of the wave propagating between the reflection point (x, y, z) and the reception point 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
is. If the maximum wave number in the _x'1 direction determined by equation (3) is _k'x1max and the maximum wave number in the _x'2 direction is _k'x2max , then variable substitution processing in the x direction is performed within the range of _-k'x1max ≦ _k'x1 ≦ _k'x1max and _-k'x2max ≦ _k'x2 ≦ _k'x2max ,

If the maximum wave number in the _y'1 direction determined by equation (4) is _k'y1max and the maximum wave number in the _y'2 direction is _k'y2max , variable substitution processing in the y direction is performed within the _{ranges -k'y1max ≦ k'y1} ≦ _k'y1max and _-k'y2max ≦k'y2 ≦ _k'y2max .

The data processing method according to claim 1 . In the variable substitution process in the x direction, the ranges of k'x1 _,nyq ≦| _k'x1 |≦ _k'x1max and _k'x2,nyq ≦| _k'x2 |≦ _k'x2max are selectively excluded, and the variable substitution process in the x direction is performed;
In the variable substitution process in the y direction, the ranges of k' _y1,nyq ≦|k' _y1 |≦k' _y1max and k' _y2,nyq ≦|k' _y2 |≦k' _y2max are selectively excluded, and the variable substitution process in the y direction is performed.
The data processing method according to claim 2 . The wave number and wave vector components of the wave satisfy the formula (5).
The data processing method according to any one of claims 1 to 3.
A measurement system for analyzing scattered waves of waves radiated to an object, comprising:
A transceiver unit,
a transmitting unit that radiates the wave toward the object from a plurality of transmitting points p ₁ (x' ₁ , y' ₁ , z' ₁ ) that are two-dimensionally arranged on an xy plane;
a receiving unit that receives the scattered waves reflected with reflectivity f(x, y, z) at a reflection point (x, y, z) on the object as measurement values s(x' ₁ , x' ₂ , y' 1 , y' ₂ , z' 1 , z' ₂ , k) at a plurality of receiving points p _{2 (} x' ₂ , y' ₂ , z' ₂ ) that are two-dimensionally arranged on _{the xy} plane;
a transceiver unit having
A processing device,
a step of performing a quadruple Fourier transform on the measured values s( _x'1 , _x'2 , _y'1 , _y'2 , _z'1 , _z'2 , k) using equation (1) to obtain S( _k'x1 , _k'x2 , _k'y1 , _k'y2 , _z'1 , _z'2 , k);

Let _k'x1,nyq be the Nyquist wave number in the _x'1 direction determined by the measurement interval in the x direction, and _k'x2,nyq be the Nyquist wave number in the _x'2 direction. Then, perform variable substitution processing in the x direction within the range of _-kx1 ≦ _k'x1 ≦ _kx1 and _-kx2 ≦ _k'x2 ≦ _kx2 (where k'x1 _,nyq ≦ _kx1 and _k'x2,nyq ≦ _kx2 );
where the Nyquist wave number in the _y'1 direction determined by the measurement interval in the y direction is k'y1 _,nyq and the Nyquist wave number in the _y'2 direction is _{k'y2 ,nyq} , a procedure for performing variable substitution processing in the _y direction within the range of _-ky1 ≦ _k'y1 ≦ _ky1 and -ky2 ≦k'y2 ≦ _ky2 (where k'y1 _,nyq ≦ _{ky1 and k'y2,} nyq ≦ _ky2 );
A procedure of calculating the reflectance f(x, y, z) by triple inverse Fourier transform using equation (2);

a processing unit for executing the
A measurement system having:
however,
z' ₁ = z' ₂ = 0
k is the wave number of the propagating wave,
k' _x1 , k' _y1 , k' _z1 are components of the wave vector of the spherical wave of the wave propagating between the transmitting point p ₁ (x' ₁ , y' ₁ , z' ₁ ) and the reflecting point (x, y, z),
k' _x2 , k' _y2 , k' _z2 are components of the wave vector of the spherical wave of the wave propagating between the reflection point (x, y, z) and the reception point 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
is. The processing device includes:
If the maximum wave number in the _x'1 direction determined by equation (3) is _k'x1max and the maximum wave number in the _x'2 direction is _k'x2max , then variable substitution processing in the x direction is performed within the range of _-k'x1max ≦ _k'x1 ≦ _k'x1max and _-k'x2max ≦ _k'x2 ≦ _k'x2max ,

If the maximum wave number in the _y'1 direction determined by equation (4) is _k'y1max and the maximum wave number in the _y'2 direction is _k'y2max , variable substitution processing in the y direction is performed within the _{ranges -k'y1max ≦ k'y1} ≦ _k'y1max and _-k'y2max ≦k'y2 ≦ _k'y2max .

The measurement system according to claim 5 . The processing device includes:
a step of selectively excluding a range of k'x1 _,nyq ≦| _k'x1 |≦ _k'x1max and k'x2 _,nyq ≦| _k'x2 |≦ _k'x2max in the variable substitution process in the x direction, and performing the variable substitution process in the x direction;
a step of selectively excluding a range of k' _y1,nyq ≦|k' _y1 |≦k' _y1max and k' _y2,nyq ≦|k' _y2 |≦k' _y2max in the variable substitution process in the y direction;
To execute
The measurement system of claim 6 . The wave number and wave vector components of the wave satisfy the formula (5).
The measurement system according to any one of claims 5 to 7.
A program for analyzing scattered waves of waves radiated to an object,
The procedure involves performing a quadruple Fourier transform on the measured values s( _x'1 , _x'2 , _y'1 , _y'2 , _z'1 , _z'2 , k) using equation (1) to obtain _S ( _k'x1 , k'x2, _k'y1 , _k'y2 , _z'1 , _z'2 , k), and


Let _k'x1,nyq be the Nyquist wave number in the _x'1 direction determined by the measurement interval in the x direction, and _k'x2,nyq be the Nyquist wave number in the _x'2 direction. Then, perform variable substitution processing in the x direction within the range of _-kx1 ≦ _k'x1 ≦ _kx1 and _-kx2 ≦ _k'x2 ≦ _kx2 (where k'x1 _,nyq ≦ _kx1 and _k'x2,nyq ≦ _kx2 );
where the Nyquist wave number in the _y'1 direction determined by the measurement interval in the y direction is k'y1 _,nyq and the Nyquist wave number in the _y'2 direction is _{k'y2 ,nyq} , a procedure for performing variable substitution processing in the _y direction within the range of _-ky1 ≦ _k'y1 ≦ _ky1 and -ky2 ≦k'y2 ≦ _ky2 (where k'y1 _,nyq ≦ _{ky1 and k'y2,} nyq ≦ _ky2 );
A procedure for calculating the reflectance f(x, y, z) by triple inverse Fourier transform using equation (2),

A program that causes a computer to execute the following.
however,
z' ₁ = z' ₂ = 0
k is the wave number of the propagating wave,
k' _x1 , k' _y1 , k' _z1 are the components of the wave vector of the spherical wave of the wave propagating between the transmitting point p ₁ (x' ₁ , y' ₁ , z' ₁ ) and the reflecting point (x, y, z),
k' _x2 , k' _y2 , k' _z2 are components of the wave vector of the spherical wave of the wave propagating between the reflection point (x, y, z) and the reception point 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
is. where the maximum wave number in the _x'1 direction determined by equation (3) is _k'x1max and the maximum wave number in the _x'2 direction is _k'x2max , a procedure for performing variable substitution processing in the x direction within the range of _-k'x1max ≦ _k'x1 ≦ _k'x1max and _-k'x2max ≦ _k'x2 ≦ _k'x2max ;

When the maximum wave number in the _y'1 direction determined by equation (4) is _k'y1max and the maximum wave number in the _y'2 direction is _k'y2max , a procedure of performing variable substitution processing in the y direction within the ranges _-k'y1max ≦ _k'y1 ≦ _k'y1max and _-k'y2max ≦ _k'y2 ≦k'y2max _;

The program according to claim 9, further causing a computer to execute the following: a step of selectively excluding a range of k'x1 _,nyq ≦| _k'x1 |≦ _k'x1max and k'x2 _,nyq ≦| _k'x2 |≦ _k'x2max in the variable substitution process in the x direction, and performing the variable substitution process in the x direction;
a step of selectively excluding a range of k' _y1,nyq ≦|k' _y1 |≦k' _y1max and k' _y2,nyq ≦|k' _y2 |≦k' _y2max in the variable substitution process in the y direction;
The program according to claim 10, further causing a computer to execute the following: The measurement value s( _x'1 , _x'2 , _y'1 , _y'2 , _z'1 , _z'2 , k) is a value of the scattered wave received at a plurality of receiving points _p2 (x'2, _y'2 , _z'2 ) arranged two-dimensionally on the xy plane, which is a wave radiated to the object from a plurality of transmitting points p1 ( _x'1 , y'1, z'1) arranged two-dimensionally on the xy plane and reflected at the reflecting point (x, y, _z ) on the object with the reflectivity _f ( _x , _y , z).
The program according to any one of claims 9 to 11. The wave number and wave vector components of the wave satisfy the formula (5).
The program according to any one of claims 9 to 11 .