JP2014182032A - Rotation estimation system of rotary flying body, rotation estimation method, and rotation estimation program - Google Patents

Abstract

A rotation parameter is accurately estimated.
A rotation estimation system (1) irradiates a sphere (10) flying with rotation with electromagnetic waves from three directions, receives reflected electromagnetic waves from the sphere (10), and receives reception strength of the three-axis reflected electromagnetic waves. A measuring device 5 that measures time series information about the frequency, a frequency domain conversion unit 61 of a calculation device 6 that converts three-axis time series information into frequency information, and peak intervals Δf _x , Δf _y in the three-axis frequency information, a threshold processing operation section 62 for specifying a Delta] f _z, respectively, based on the peak interval triaxial specified, the angular velocity omega _x of the three axes of the spheres _10, omega _y, and angular separation unit 64 for estimating the omega _z, three An angular velocity vector deriving unit 66 that outputs the magnitude of the angular velocity vector ω having the angular velocity of the axis as the rotational speed of the sphere 10 and outputs the direction of the angular velocity vector ω as the rotational axis position of the sphere 10.
[Selection] Figure 2

Description

  The present invention relates to a rotation estimation system, a rotation estimation method, and a rotation estimation program for a rotating projectile.

  2. Description of the Related Art Conventionally, for example, a technique for measuring a rotation parameter of a rotating flying object that flies with rotation, such as a trajectory, a rotation speed, and a rotation direction of a baseball ball at the time of pitching or a golf ball after launch, is known. For example, Patent Document 1 discloses a technique for obtaining the rotational speed and trajectory of a sphere from a Doppler shift.

Japanese Patent No. 4865735

  However, the conventional method for measuring the rotation parameters of the rotating projectile has room for further improvement in the measurement accuracy of the rotation parameters.

  The present invention has been made in view of the above, and an object of the present invention is to provide a rotation projectile rotation estimation system, a rotation estimation method, and a rotation estimation program capable of accurately estimating the rotation parameters of the rotation projectile. .

  In order to solve the above problems, a rotation estimation system for a rotating projectile according to the present invention is a rotation estimation system that estimates a rotation speed and a rotation axis position of a rotating projectile that flies with rotation, and the rotation flight Measuring means for irradiating the body with electromagnetic waves from three axial directions, receiving reflected electromagnetic waves from the rotating flying object with the three axes, and measuring time-series information regarding the reception intensity of the three-axis reflected electromagnetic waves; Conversion means for converting the three-axis time series information measured by the measurement means into frequency information, and peak interval specifying means for specifying peak intervals in the three-axis frequency information converted by the conversion means, respectively. And angular velocity estimation means for estimating the three-axis angular velocity of the rotating projectile based on the three-axis peak interval specified by the peak interval specifying means, The magnitude of the angular velocity vector whose component is the angular velocity of the three axes estimated by the degree estimating means is output as the rotational speed of the rotating projectile, and the direction of the angular velocity vector is output as the rotational axis position of the rotating projectile. And an output means.

  Further, in the above rotation estimation system, the angular velocity estimation means is configured so that each of the peak intervals of the three axes is around the other two axes excluding an axis serving as a viewpoint when the rotating projectile is viewed from each axial direction. It is preferable to estimate the three-axis angular velocity of the rotating projectile based on the characteristic that the information corresponds to the rotational velocity obtained by combining the angular velocity of the rotation.

In the above rotation estimation system, when the angular velocity estimation means represents the peak intervals of the three axes as Δf _x , Δf _y , Δf _z , and the angular velocity of the three axes as ω _x , ω _y , and ω _z , respectively. It is preferable to estimate the angular velocity of the three axes using the following mathematical formula.

  In addition, the rotation estimation system obtains the orbital change of the rotating projectile from the three-axis directions and compares it with the expected orbital change from the flight direction of the rotating projectile. It is preferable to include a determination unit that determines a positive / negative direction of the estimated angular velocity.

In the above rotation estimation system, the determination means sets the measurement start point of the reflected electromagnetic wave from the rotating projectile as P ₀ (x ₀ , y ₀ , z ₀ ) and the arbitrary measurement point before the trajectory change as P 0. _h (x _h , y _h , z _h ), p _i (x _i , y _i , z _i ) as the measurement end point after the trajectory change, and the three points from the measurement start point p ₀ to the measurement end point p _i When the axis trajectory change amount is expressed as Δx, Δy, Δz, gravitational acceleration is expressed as g, initial velocity is expressed as v ₀ , and wind speeds of the three axes are expressed as w _x , w _y , w _z It is preferable to determine the direction.

  Further, in the above rotation estimation system, the peak interval specifying unit extracts a peak having a predetermined threshold value or more from the triaxial frequency information converted by the conversion unit, and calculates a peak interval between adjacent peaks. Preferably, the calculated peak interval is extracted with a high frequency, and the peak interval is specified based on the extracted peak interval.

  The rotation estimation system includes a display unit, the rotating projectile is a ball game ball, and based on the rotation speed and the rotation axis position of the ball game ball output by the output unit, It is preferable to display information on the type or trajectory of the ball on the display means.

  Similarly, in order to solve the above-mentioned problem, the rotation estimation method of the rotating projectile according to the present invention is a rotation estimation method for estimating the rotation speed and the rotation axis position of the rotating projectile that flies with rotation, Irradiate the rotating projectile with electromagnetic waves from three directions, receive the reflected electromagnetic waves from the rotating projectile with the three axes, and measure time-series information regarding the reception intensity of the three-axis reflected electromagnetic waves. A measurement step; a conversion step for converting the three-axis time-series information measured in the measurement step into frequency information; and a peak for identifying a peak interval in the three-axis frequency information converted in the conversion step. Estimating the three-axis angular velocity of the rotating projectile based on the interval specifying step and the three-axis peak interval specified by the peak interval specifying step A velocity estimation step, the magnitude of an angular velocity vector having the three-axis angular velocity estimated in the angular velocity estimation step as a rotational velocity of the rotating projectile, and the direction of the angular velocity vector of the rotating projectile And an output step of outputting as the rotation axis position.

  Similarly, in order to solve the above-mentioned problem, a rotation estimation program for a rotating projectile according to the present invention is a rotation estimation program for estimating the rotation speed and rotation axis position of a rotating projectile that flies with rotation, Irradiate the rotating projectile with electromagnetic waves from three directions, receive the reflected electromagnetic waves from the rotating projectile with the three axes, and measure time-series information regarding the reception intensity of the three-axis reflected electromagnetic waves. A measurement function, a conversion function for converting the time-series information of the three axes measured by the measurement function into frequency information, and a peak for identifying a peak interval in the frequency information of the three axes converted by the conversion function, respectively An angular velocity estimator that estimates the three-axis angular velocity of the rotating projectile based on the interval identifying function and the three-axis peak interval identified by the peak interval identifying function The angular velocity vector having the three-axis angular velocity estimated by the angular velocity estimating function as the rotational velocity of the rotating projectile, and the direction of the angular velocity vector as the rotational axis position of the rotating projectile And an output function for outputting as a computer.

  The rotation estimation system, the rotation estimation method, and the rotation estimation program according to the present invention have an effect that the rotation parameters of the rotation projectile can be accurately estimated.

FIG. 1 is a schematic diagram showing a schematic configuration of a rotation estimation system for a rotating projectile according to an embodiment of the present invention. FIG. 2 is a functional block diagram of the control device in FIG. FIG. 3 is a flowchart showing the derivation process of the angular velocity vector of the sphere executed by the rotation estimation system of the present embodiment. FIG. 4 is a diagram showing the change in the angle of the radar cross-sectional area per rotation of the soft ball for baseball. FIG. 5 is a diagram showing the change in the angle of the radar cross-section per rotation of a sphere with a built-in corner reflector. FIG. 6 is a diagram for explaining a radar sphere tracking method. FIG. 7 is a diagram illustrating an example of a temporal change in the reception intensity on the x-axis. FIG. 8 is a diagram illustrating an example of a temporal change in reception intensity on the y-axis. FIG. 9 is a diagram illustrating an example of a temporal change in reception intensity on the z axis. FIG. 10 is a diagram illustrating an example of frequency information of x-axis reception intensity. FIG. 11 is a diagram illustrating an example of frequency information of reception intensity on the y-axis. FIG. 12 is a diagram illustrating an example of frequency information of reception intensity on the z axis. FIG. 13 is a diagram illustrating a case where a spectrum peak that originally occurs cannot be obtained in the frequency information of reception intensity. FIG. 14 is a diagram illustrating a case where a spectrum of a component other than the rotation of the sphere is obtained in the frequency information of the reception intensity. FIG. 15 is a flowchart showing a subroutine of threshold value processing and mode value processing in step S103 of the flowchart of FIG. FIG. 16 is a diagram illustrating an example of a threshold setting method in threshold processing. FIG. 17 is a schematic diagram showing a change in the trajectory of a sphere in the xyz coordinate system. FIG. 18 is a diagram in which the trajectory of the sphere in FIG. 17 is projected onto the zx plane coordinates. FIG. 19 is a diagram in which the trajectory of the sphere in FIG. 17 is projected onto the xy plane coordinates. FIG. 20 is a schematic diagram showing an outline of an angular velocity vector. FIG. 21 is a schematic diagram showing a configuration of an example in which a baseball curve ball is measured from three directions. FIG. 22 is a diagram illustrating a temporal change in the x-axis reception intensity acquired in the embodiment. FIG. 23 is a diagram illustrating a temporal change in the reception intensity on the y-axis acquired in the example. FIG. 24 is a diagram illustrating a temporal change in the reception intensity on the z-axis acquired in the example. FIG. 25 is a diagram illustrating a frequency spectrum of x-axis reception intensity. FIG. 26 is a diagram illustrating a frequency spectrum of the reception intensity on the y-axis. FIG. 27 is a diagram illustrating a frequency spectrum of the received intensity on the z-axis. FIG. 28 is a diagram illustrating a frequency spectrum of the reception intensity on the y-axis obtained by reacquisition of the peak interval. FIG. 29 is a diagram illustrating a frequency spectrum of z-axis reception intensity obtained by reacquiring the peak interval. FIG. 30 is a diagram illustrating changes in the trajectory of the sphere in the x-axis in the example. FIG. 31 is a diagram illustrating a change in the y-axis trajectory of a sphere in the example. FIG. 32 is a diagram showing a change in the z-axis trajectory of a sphere in the example.

  DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS Embodiments of a rotational projectile rotation estimation system according to the present invention will be described below with reference to the drawings. In the following drawings, the same or corresponding parts are denoted by the same reference numerals, and the description thereof will not be repeated.

[Embodiment]
First, the configuration of the rotation estimation system 1 according to the present embodiment will be described with reference to FIGS. FIG. 1 is a schematic diagram showing a schematic configuration of a rotation estimation system according to an embodiment of the present invention, and FIG. 2 is a functional block diagram of a control device in FIG.

(Configuration of rotation estimation system)
As shown in FIG. 1, the rotation estimation system 1 irradiates a rotating flying object 10 that rotates and flies with electromagnetic waves from three axial directions, and uses the received reflected electromagnetic waves to rotate the rotational speed of the rotating flying object. And the direction of the rotation axis. In the present embodiment, a sphere 10 such as a baseball is illustrated as an example of the rotating flying object 10. In the present embodiment, the x axis, the y axis, and the z axis that are orthogonal to each other with the sphere 10 as the origin are set as the three axes that irradiate the sphere 10 with electromagnetic waves and receive the reflected electromagnetic waves. Here, the z-axis is defined as a positive direction vertically upward, and the x-axis and y-axis directions are horizontal directions orthogonal to the z-axis. It is assumed that the sphere 10 starts to fly in the positive direction of the x axis (the direction indicated as “sphere movement direction” in FIG. 1) at the initial velocity along the x axis.

  As shown in FIG. 1, the rotation estimation system 1 includes an x-axis radar 2, a y-axis radar 3, a z-axis radar 4, a measurement device 5 (measurement means), a calculation device 6, a recording device 7, and an output device 8 (display means). ).

  The x-axis radar 2, the y-axis radar 3, and the z-axis radar 4 are arranged on the x, y, and z axes, respectively, and irradiate the sphere 10 with electromagnetic waves 11, 13, and 15 from the x, y, and z axis directions. Thus, the reflected electromagnetic waves 12, 14, and 16 reflected by the sphere 10 are received. As the electromagnetic wave, a continuous wave or a pulse wave is used.

  The measuring device 5 measures the time variation (time series information) of the reception intensity of the reflected electromagnetic waves 12, 14, and 16 obtained by the x-axis radar 2, the y-axis radar 3, and the z-axis radar 4, the radio wave arrival time, the reception frequency, and the like. taking measurement.

  The calculation device 6 performs frequency domain transformation, threshold processing, mode processing, angular velocity separation, angular velocity vector derivation, Doppler shift derivation, Doppler shift, or initial velocity derivation using radio wave arrival time of the received signal obtained by the measurement device 5. , Orbital change derivation, angular velocity ambiguity removal, etc.

  As shown in FIG. 2, the calculation device 6 has a frequency domain conversion unit 61 (conversion unit), a threshold processing calculation unit 62 (peak interval specifying unit), and a mode value processing calculation unit 63 (peak interval specifying) with respect to the above functions. Means), an angular velocity separation unit 64 (angular velocity estimation unit), an angular velocity direction determination unit 65 (determination unit), and an angular velocity vector derivation unit 66 (output unit).

  The frequency domain conversion unit 61 converts the time series information of the x, y, and z axes acquired by the measurement device 5 into frequency information.

  The threshold processing calculator 62 extracts peaks that are equal to or greater than a predetermined threshold from the x, y, and z axis frequency information converted by the frequency domain converter 61, and calculates a peak interval between adjacent peaks.

The mode value processing calculation unit 63 extracts a high frequency among the peak intervals calculated by the threshold value processing calculation unit 62, and specifies the peak intervals Δf _x , Δf _y , Δf _z based on the extracted peak intervals. To do.

The angular velocity separation unit 64 calculates the angular velocities ω _x , ω _y , and ω _z of the sphere 10 based on the identified peak intervals Δf _x , Δf _y , and Δf _z of the x, y, and z axes. presume.

The angular velocity direction determination unit 65 determines the positive / negative direction of the angular velocities ω _x , ω _y , and ω _z estimated by the angular velocity separation unit 64.

The angular velocity vector deriving unit 66 derives an angular velocity vector ω having the components of the x, y, and z-axis angular velocities ω _x , ω _y , and ω _z estimated by the angular velocity separating unit 64. Details of processing of each functional block of the computing device 6 will be described later with reference to FIG.

  The recording device 7 is a device that stores the value obtained by the calculation device 6, and more specifically, information on the angular velocity vector derived by the angular velocity vector deriving unit 66, and rotation of the sphere 10 that can be derived from the angular velocity vector. Saves information such as speed and rotation axis direction.

  The output device 8 outputs the value obtained by the calculation device 6 and performs angular velocity vector display on the display, Internet distribution, and the like.

  The measuring device 5, the computing device 6, the recording device 7, and the output device 8 physically include a CPU (Central Processing Unit), a RAM (Random Access Memory), a ROM (Read Only Memory), and an EEPROM (Electrically Erasable Programmable). A microcomputer having hardware such as a storage device such as an only memory (HDD) and a hard disk drive (HDD), an interface for communicating with each part inside and outside the device, an input device such as a switch, a keyboard and a mouse, and a display device such as a display. It is. All or some of the functions of the measurement device 5, the calculation device 6, the recording device 7, and the output device 8 described in the present embodiment allow a predetermined application program to be read on hardware such as a CPU, RAM, and ROM. Thus, the interface, input device, display device, and the like are operated under the control of the CPU, and data is read and written in the RAM, ROM, and storage device.

  In addition, application programs (particularly the measurement device 5, the frequency domain conversion unit 61, the threshold value processing calculation unit 62, the mode value processing calculation of the measurement device 5, the calculation device 6, and the output device 8). The function of the unit 63, the angular velocity separation unit 64, the angular velocity direction determination unit 65, and the angular velocity vector derivation unit 66) may be stored in a computer-readable recording medium, or may be configured as a program product. Here, the “recording medium” is any memory card, USB memory, SD card, flexible disk, magneto-optical disk, ROM, EPROM, EEPROM, CD-ROM, MO, DVD, Blu-ray Disc, etc. Of “portable physical media”. The application program may be stored in an application program server connected to the rotation estimation system 1 via an arbitrary network, and may be downloaded in whole or in part as necessary. .

(Operation of the rotation estimation system)
Next, with reference to FIGS. 3-20, operation | movement of the rotation estimation system 1 of this embodiment is demonstrated. FIG. 3 is a flowchart showing the derivation process of the angular velocity vector of the sphere executed by the rotation estimation system of the present embodiment. The process of the flowchart of FIG. 3 is performed, for example, every predetermined time while the sphere 10 to be measured is flying.

(Measurement of reflected electromagnetic wave)
In step S101, the measuring device 5 uses the x-axis radar 2, the y-axis radar 3, and the z-axis radar 4 to irradiate the sphere 10 with the electromagnetic waves 11, 13, and 15 from the x, y, and z axis directions. The reflected electromagnetic waves 12, 14, and 16 reflected by the sphere 10 are received for a predetermined time. The measuring device 5 uses the reflected electromagnetic waves 12, 14, and 16 for a predetermined time received by the x-axis radar 2, the y-axis radar 3, and the z-axis radar 4, and reflects the reflected electromagnetic waves 12, 14 on the x, y, and z axes. , 16, the time change (time series information) of the received intensity is measured.

Here, with reference to FIGS. 4 and 5, the reason for measuring the temporal change of the reception intensity of the reflected electromagnetic waves 12, 14, and 16 in this embodiment will be described. FIG. 4 is a diagram showing the change in the angle of the radar cross-sectional area per rotation of the soft ball for baseball. The horizontal axis of FIG. 4 shows the rotation angle (deg) of a soft ball for baseball as an example of the sphere 10, and the vertical axis shows the radar cross section (dB · m ² ). “Radar cross section (RCS)” is an amount representing the degree to which an object scatters electromagnetic waves, and is a parameter that varies according to the received intensity of received electromagnetic waves. As shown in FIG. 4, the radar cross-sectional area changes non-uniformly according to the rotation angle of the soft ball for baseball. Therefore, it can be considered that there is a strong correlation between the phase change of the rotation of the sphere 10 and the time change of the received intensity of the reflected electromagnetic waves 12, 14, 16. The rotation estimation system 1 of this embodiment finally estimates the rotation axis direction and rotation speed of the sphere 10. In order to perform these estimations with high accuracy, it is preferable to include information relating to the rotational phase of the sphere 10 as input information for estimation. Therefore, in the present embodiment, the time variation of the reception intensity of the reflected electromagnetic waves 12, 14, 16 that has a strong correlation with the rotational phase of the sphere 10 is used.

  If the RCS to be measured is small and difficult to distinguish from other objects, for example, a “game ball” disclosed in Japanese Patent No. 5142366 may be used. This game ball is a spherical reflector configured such that eight corner reflectors are orthogonal to each other. FIG. 5 is a diagram showing a change in the angle of the radar cross-sectional area per rotation of the game ball. By applying the configuration of the game ball to the measurement target of the rotation estimation system 1, it is possible to achieve a significant improvement in RCS and stability of angle change. Further, the reflected electromagnetic wave can be easily detected by the radar.

  When the distance between the x-axis radar 2, the y-axis radar 3, the z-axis radar 4 and the sphere 10 is short, the x-axis radar 2, the y-axis radar 3, and the z-axis radar according to the movement of the sphere 10. 4 may be moved, and the sphere 10 may be tracked so that the sphere 10 is always arranged at the three-axis origin.

ここで、ｃは光速度である。 As a tracking method of the sphere, for example, the time variation of the arrival time from the irradiation of the electromagnetic waves 11, 13, 15 to the reception of the reflected electromagnetic waves 12, 14, 16 can be used. FIG. 6 is a diagram for explaining a radar sphere tracking method. As shown in FIG. 6, x-axis radar 2, y-axis radar 3, the time variation of the arrival time in the z-axis radar 4 respectively Delta] t _x, Delta] t _y, When Delta] t _z, x-axis radar 2, y-axis radar 3 The time variation ΔR _x , ΔR _y , ΔR _z of the distance between the z-axis radar 4 and the sphere 10 can be calculated by the following equation (1).

Here, c is the speed of light.

In the present embodiment, the sphere 10 is set so as to be always located at the origin during the flight in the xyz coordinate system including the x, y, and z axes. ΔR _x , ΔR _y , ΔR _z calculated by the above equation (1) are movement amounts in the respective axis directions from the origin of the sphere 10 in the xyz coordinate system. Therefore, by tracking the position of the x-axis radar 2, y-axis radar 3, and z-axis radar 4 according to the amount of movement, the sphere 10 can be tracked while maintaining the origin of the xyz coordinate system. it can.

For example, in the case of the x-axis radar 2, when the movement amount ΔR _{y in the} y-axis direction is positive, it is in the positive direction of the y-axis (the direction indicated as “+ R _y ” in FIG. 6), and ΔR _y is negative. , The position of the x-axis radar 2 is moved by | ΔR _y | in the negative direction of the y-axis (the direction indicated by “−R _y ” in FIG. 6). Further, when the movement amount ΔR _z in the z-axis direction is positive, the z-axis is in the positive direction (the direction indicated as “+ R _z ” in FIG. 6), and when ΔR _z is negative, the z-axis is in the negative direction. In the direction indicated by “−R _z ” in FIG. 6, the position of the x-axis radar 2 is moved by | ΔR _z |. Similarly, the y-axis radar 3 moves in the + R _x direction or −R _x direction by | ΔR _x |, and moves in the + R _z direction or −R _z direction by | ΔR _z |. The z-axis radar 4 moves by | ΔR _x | in the + R _x direction or −R _x direction, and moves by | ΔR _y | in the + R _y direction or −R _y direction.

  Further, noise generated from other than the sphere 10 included in the reflected electromagnetic waves 12, 14, and 16 received by the x-axis radar 2, the y-axis radar 3, and the z-axis radar 4 is, for example, a moving target indicator (MTI). ) Or the like.

  An example of the time change of the received intensity on the x, y, and z axes measured in step S101 is shown in FIGS. FIG. 7 is a diagram illustrating an example of a temporal change in the reception intensity on the x axis, FIG. 8 is a diagram illustrating an example of a temporal change in the reception intensity on the y axis, and FIG. It is a figure which shows an example of a time change. The vertical axis in FIGS. 7 to 9 represents the relative reception intensity (dB), and the horizontal axis represents the measurement time (s (seconds)). The graphs shown in FIGS. 7 to 9 are the changes over time in the relative received intensity on the x-, y-, and z-axes of the softball for baseball rotated once per second around the z-axis. Here, the “relative reception strength” is a relative value of the reception strength based on the minimum value of the reception strength. When the process of step S101 is completed, the process proceeds to step S102.

(Frequency conversion processing)
In step S102, the time change (time-series information) of the received intensity on the x, y, and z axes is converted into frequency domain information (frequency information) by the frequency domain conversion unit 61 of the calculation device 6. The frequency domain transforming unit 61 can use a known frequency transforming method such as Fourier transform, wavelet transform, Wigner distribution or the like.

  In addition, the measurement values to be converted to the frequency domain are restricted in that the values of the measurement start point and the end point must be the same and the period within the measurement range must be an integer. However, in actual measurement, it is impossible to always satisfy this condition, and if this is ignored and converted into the frequency domain, problems such as a decrease in frequency resolution and occurrence of extra values around the peak occur. Therefore, a window function is used to solve this problem. A window function refers to a function that becomes 0 outside a certain finite interval. By applying this to a signal, the values of the measurement start point and end point become the same and the period becomes an integer, so that the frequency resolution is improved. There are various types of window functions, such as a Hamming window, a Hanning window, and a rectangular window. The frequency domain conversion unit 61 can use these known window functions.

  An example of the frequency information of the received intensity on the x, y, and z axes converted in step S102 is shown in FIGS. FIG. 10 is a diagram illustrating an example of frequency information of x-axis reception strength, FIG. 11 is a diagram of an example of frequency information of reception strength of the y-axis, and FIG. It is a figure which shows an example of frequency information. 10 to 12, the vertical axis represents the power spectrum (dB), and the horizontal axis represents the frequency (Hz / rps). The frequency information shown in FIGS. 10 to 12 is obtained by frequency-converting the temporal change in the received intensity on the x, y, and z axes shown in FIGS. In this example, Fourier transform is used as the frequency conversion method, and Hanning window is used as the window function. As shown in FIGS. 10 to 12, in the frequency domain, the reception intensity on the x, y, and z axes can obtain a periodic peak. Note that the frequency interval between adjacent peaks at this time coincides with the rotational speed of each axis. When the process of step S102 is completed, the process proceeds to step S103.

(Threshold processing and mode processing)
In step S <b> 103, threshold processing is performed by the threshold processing unit 62 of the computing device 6, and mode processing is performed by the mode processing unit 63. As shown in FIGS. 10 to 12, the frequency spectrum of the received intensity on the x, y, and z axes has a periodic peak such that the peak interval coincides with the rotational speed on each axis. However, for example, when the peak of the frequency that is originally generated cannot be obtained as in the region indicated by the dotted line in FIG. 13, or is not related to the rotation of the sphere 10 as in the region indicated by the dotted line in FIG. 14. There may be cases where peaks occur. In this case, the correct peak interval cannot be obtained from the frequency spectrum, which may adversely affect the estimation of the rotational speed. The threshold value processing and the mode value processing in this step are performed in order to remove elements that are unnecessary for obtaining a correct peak interval, as indicated by dotted lines in FIGS.

  In step S103, the subroutine shown in FIG. 15 is executed. FIG. 15 is a flowchart showing a subroutine of threshold value processing and mode value processing in step S103. In FIG. 15 and the following description, the x, y, and z axis symbols are omitted for convenience, but the series of processing shown in the flowchart of FIG. 15 is performed in the x, y, and z axis acquired in step S102. This is performed individually for each frequency information of the reception intensity.

  The process in step S103 will be described in detail with reference to the flowchart of FIG. 15. First, an arbitrary threshold value is set by the threshold value processing calculation unit 62, and threshold value processing for extracting a peak equal to or higher than this threshold value from the frequency information of the received intensity is performed. (Step S201).

Here, an example of the threshold setting method will be described with reference to FIG. FIG. 16 is a diagram illustrating an example of a threshold setting method in threshold processing. The vertical axis in FIG. 16 indicates the power spectrum (dB), and the horizontal axis indicates the frequency (Hz / rps). In general, the frequency spectrum of the received intensity of the radar includes a direct current component in a region around the frequency 0. As the received power of the radar increases, the direct current component increases. Due to the influence of such a direct current component, for example, as shown in FIG. 16, the peak value of the power spectrum at the frequency 0 may become extremely larger than other components. In order to extract peaks with high accuracy in threshold processing, it is desirable to set the threshold by excluding the influence of such DC components. Therefore, for example, as shown in FIG. 16, to exclude the range of up to 0 to F _a that contains a DC component, to identify the frequency f _max of the maximum value of the peak of the above-domain frequency f _a. Then, by calculating the average value of the power spectrum in the interval of f _max from the frequency f _a, it sets the average value as a threshold value. The frequency fa at the boundary of the region including the DC component is, for example, _a frequency at which the differential value of the power spectrum switches from negative to positive, that is, a frequency at which the power spectrum decreases from the maximum value at frequency 0 and starts to increase. Note that the threshold setting method described with reference to FIG. 16 is an example, and other methods can be used as appropriate.

Returning to FIG. 15, the threshold value processing calculation unit 62 obtains peaks f ₁ , f ₂ , f ₃ ,... F _i from the frequency spectrum of the received intensity of each axis (step S 202). For example, as illustrated in FIGS. 10 to 12, the threshold processing calculation unit 62 sets a predetermined threshold for each frequency spectrum of the x, y, and z axes, and sets a frequency that peaks from the frequency spectrum equal to or higher than the threshold to f ₁ , _{f 2,} f 3, to get as is ... _{f i.} The acquired sign of each peak is given, for example, in order from the smaller frequency to the larger frequency. Then, peak intervals Δf ₁ , Δf ₂ , Δf ₃ ,... Δf _i between adjacent peaks are acquired (step S 203). For example, the first peak interval Δf ₁ can be calculated by Δf ₁ = f ₂ −f ₁ .

Subsequently, the mode value processing is performed by the mode value processing calculation unit 63. First, among the peak intervals Δf ₁ , Δf ₂ , Δf ₃ ,... Δf _i acquired in step S203, a minimum value Δf _min1 is acquired, and a peak interval in the range of Δf _min1 ≦ Δf _min1 + 2Δf _r is acquired. Let them be group 1. Here, Delta] f _r is the frequency resolution of the Fourier transform.

Next, the minimum value Δf _min2 is acquired from the peak intervals Δf ₁ , Δf ₂ , Δf ₃ ,... Δf _i obtained in step S203 by removing the group 1, and a range of Δf _min2 ≦ Δf _min2 + 2Δf _r is obtained. The peak intervals at are acquired and set as group 2.

After performing this calculation until all of the peak intervals Δf ₁ , Δf ₂ , Δf ₃ ,... Δf _i acquired in step S203 belong to any group, the group with the largest number of elements is the mode. Obtained as a value group Δf _m1 (step S204).

Next, it is determined whether or not the number of elements of the mode value group Δf _m1 is less than 4 (step S205). When the number of elements of the mode value group Δf _m1 is less than 4 (Yes in step S205), the number of peaks is insufficient and the frequency information does not have a tendency of periodic rotation. Assuming that no occurrence has occurred, the peak interval Δf of the frequency spectrum of the corresponding axis is set to 0 (step S206), and the process returns to the main flow.

On the other hand, if the number of elements of the mode value group Δf _m1 is 4 or more (No in step S205), then the number of elements of the mode value group Δf _m1 is the peak interval Δf acquired in step S203. ₁ , Δf ₂ , Δf ₃ ,... Δf _i is checked to see if it is half or more of the total number (step S 207). When the number of elements of the mode value group Δf _m1 is more than half of the whole (Yes in step S207), the average value of the peak interval values of each element of the mode value group Δf _m1 (in FIG. 15, above Δf _m1 Is set as the peak interval Δf of the frequency spectrum of the corresponding axis (step S208), and the process returns to the main flow.

When the number of elements of the mode value group Δf _m1 is less than half of the whole (No in step S207), the group having the next largest number of elements after the mode value group Δf _m1 is acquired as the second mode value group Δf _m2. (Step S209).

The mode value group Δf _m1 and the second mode value group Δf _m2 are left, and peaks relating to peak intervals belonging to other groups are deleted (step S210). For example, when the peak interval Δf ₃ (= f ₄ −f ₃ ) does not belong to any of the groups Δf _m1 and Δf _m2 , the peak f ₃ is deleted.

  When the process of step S210 is completed, the process returns to step S202, the peak interval is acquired again, and the above process is repeated.

(Angular velocity estimation)
Returning to FIG. 3, in step S <b> 104, the angular velocity separation unit 64 uses the x, y, z axis frequency spectrum peak intervals Δf _x , Δf _y , Δf _z calculated in step S 103, x, y, The z-axis angular velocities ω _x , ω _y , and ω _z are calculated. Here, the peak intervals Δf _x , Δf _y , and Δf _z of the frequency spectrum of the x, y, and z axes can be considered as information corresponding to the rotational speed when the sphere 10 is viewed from each axial direction. At this time, the rotation speed is obtained by combining the angular velocities of rotation around the other two axes excluding the viewpoint axis. For example, the peak interval Delta] f _x of the frequency spectrum of the x-axis, when projecting the rotation of the spherical body 10 in the yz plane, which is information corresponding to the rotational speed of the rotation on the yz plane. The rotational speed at this time is a combination of the angular speed ω _y due to the rotation around the y axis and the angular speed ω _z due to the rotation around the z axis. Therefore, the peak intervals Δf _x , Δf _y , Δf _z of the frequency spectrum of the x, y, z axes and the angular velocities ω _x , ω _y , ω _{z of} the x, y, z axes are expressed by the following (2) to ( 4) There is a relationship of the formula. By using this parameter, the angular velocities can be separated.

When the above equations (2) to (4) are expanded, equations (5) to (7) relating to the angular velocities ω _x , ω _y , and ω _z can be derived.

The angular velocity separation unit 64 substitutes the peak intervals Δf _x , Δf _y , Δf _z of the frequency spectra of the x, y, and z axes calculated in step S103 into the above expressions (5) to (7), The angular velocities ω _x , ω _y , and ω _z of the x, y, and z axes are calculated. When the process of step S104 is completed, the process proceeds to step S105.

(Angular velocity direction judgment)
In step S105, the angular velocity direction determination unit 65 determines the positive / negative direction of the angular velocities ω _y and ω _z calculated in step S104. Positive and negative ambiguities remain in the above equations (5) to (7) for deriving the angular velocity. The angular velocity direction determination unit 65 specifies the direction of the angular velocity using the change distance of the sphere 10. Since the Magnus effect and gravity act on the sphere 10 flying in the air, the direction of change of the flight trajectory is determined by the rotation direction. The trajectory changes at this time are simultaneously acquired from the three directions of the x, y, and z axes (“trajectory change amounts Δx, Δy, Δz” described below), and the trajectory changes expected from the flight direction of the sphere 10 (below, By comparing with “reference trajectory variation x _m , y _m , z _m ” to be described, the sign of the angular velocity is determined and the ambiguity is removed.

With reference to FIGS. 17-19, the angular velocity direction determination method in the angular velocity direction determination part 65 is demonstrated in detail. FIG. 17 is a schematic diagram showing a trajectory change of the sphere 10 in the xyz coordinate system. As described above, in the present embodiment, it is assumed that the sphere 10 starts to fly along the positive direction of the x axis. The x-axis radar 2, the y-axis radar 3, and the z-axis radar 4 acquire the spatial coordinates of the sphere 10 simultaneously with the measurement of the reception intensity of the sphere 10. As shown in FIG. 17, the measurement start point of the reflected electromagnetic wave from the sphere 10 at this time is p ₀ (x ₀ , y ₀ , z ₀ ), any before the trajectory of the sphere 10 changes due to the influence of the Magnus effect, etc. Is the measurement start point p on the x, y, and z axes, where p _h (x _h , y _h , z _h ) and the measurement end point after the trajectory change are p _i (x _i , y _i , z _i ). distance from ₀ to the measurement ending point _{p i} (the orbit change amount) respectively [Delta] x, [Delta] y, and Delta] z.

FIG. 18 is a diagram in which the trajectory of the sphere 10 in FIG. 17 is projected onto the zx plane coordinates. Based on the change in the trajectory of the sphere 10 in the zx plane coordinates shown in FIG. 18, the sign of the angular velocity ω _y of the y axis is determined. The trajectory change z _m (hereinafter also referred to as “reference trajectory change amount”) of the sphere 10 assuming that no trajectory change due to the rotation of the sphere 10 occurs is the trajectory change amount z _s according to the flight direction. orbit change amount of gravity z _f, the sum of orbital variation z _w of wind. By comparing the reference trajectory change amount z _m with the actually measured trajectory change amount Δz, the sign of the angular velocity ω _y of the y axis is determined.

The trajectory change amount z _s according to the flight direction of the z axis is calculated as follows. First, a measurement starting point p _0, acquires an arbitrary measurement point P _h before track changes, to obtain a change in the z-direction distance z _s' at the measurement point P _h than their difference. When the sphere 10 is to perform the constant-velocity rectilinear motion passing through two points of measurement points p _h and the measurement starting point p _0, track variation z _s by flying direction of the spherical body 10, it is calculated by the following equation (8) Can do.

In addition, when the sphere 10 flies in the air, a drop due to gravity occurs. The trajectory change amount z _f due to gravity can be calculated by the following equation (9), where g is the acceleration of gravity and t is the flight time of the sphere 10.

ここで、ｖ_０は初速度、ｗ_ｘはｘ軸方向の風速を表す。なお、初速度ｖ_０の導出には、ドップラーレーダを用いてドップラーシフトから算出する手法、または球体１０の軌道変化とその所要時間から算出する手法を用いることができる。 Further, the orbit change amount Δx in the x-axis direction can be calculated by the following equation (10).

Here, v ₀ represents the initial speed, and w _x represents the wind speed in the x-axis direction. In order to derive the initial velocity v ₀ , a method of calculating from the Doppler shift using Doppler radar, or a method of calculating from the trajectory change of the sphere 10 and its required time can be used.

From the above equation (9) and (10), the orbit change amount _{Z f} by gravity, can be calculated by the following equation (11).

Next, determine the orbit change the amount of z _w by the wind. In the actual environment, especially when measuring outdoors, the trajectory changes due to wind. If the wind velocity in the z-axis direction at this time is w _z , the trajectory change amount z _w due to this wind force can be calculated by the following equation (12).

Further, the flight time t of the sphere 10 is expressed by using the initial velocity v ₀ , the wind force w _{x in} the x-axis direction, and the orbit change Δx in the x-axis direction based on the above equation (10), and is substituted into the equation (12). Then, the trajectory change amount z _w from wind can be represented by the following equation (13).

Since the reference trajectory change amount z _{m in} the z-axis direction of the sphere 10 is the sum of the trajectory change amount z _s according to the flight direction, the trajectory change amount z _f due to gravity drop, and the trajectory change amount z _w due to wind force, as described above. Based on the above equations (8), (11), and (13), it can be expressed by the following equation (14).

The reference trajectory change amount z _{m in} the z-axis direction of the sphere 10 thus calculated serves as a threshold value for determining whether the y-axis angular velocity ω _y of the sphere 10 is positive or negative, and this z _m and the actually measured z-axis The direction of the angular velocity ω _{y on} the y axis can be determined by comparing the amount of change in the direction of trajectory Δz.

When the trajectory change amount Δz is larger than the reference trajectory change amount z _m (Δz> z _m ), that is, when the sphere 10 is above the lower end of z _{m in} FIG. 18, the sphere 10 is lifted by the Magnus force. It is in. The condition for receiving this force is backspin, and the sphere 10 rotates clockwise when viewed from the positive direction side of the y-axis (the back side in FIG. 18). Accordingly, when Δz> z _m , the angular velocity ω _{y of the} y axis is in the negative direction.

On the other hand, according to if there if trajectory variation Delta] z is the reference trajectory change amount _{z m} is smaller than (Δz _<z m), that is, the spheres 10 as shown in FIG. 18 from the lower end of the _{z m} downward, Delta] z and _{z m} Rotation It is the sum of downward changes. The condition for the downward change is top spin, and the sphere 10 rotates counterclockwise when viewed from the positive side of the y-axis. Therefore, when Δz <z _m , the y-axis angular velocity ω _y is in the positive direction. From the above, the sign of the y-axis angular velocity ω _y can be determined by the following equation (15).

FIG. 19 is a diagram in which the trajectory of the sphere 10 in FIG. 17 is projected onto the xy plane coordinates. Based on the trajectory change of the sphere 10 in the xy plane coordinates shown in FIG. 19, the sign of the z-axis angular velocity ω _z is determined. Y-axis direction of the orbit change amount of the spherical body 10 when the track changes due to the rotation of the sphere 10 is assumed not to occur y _{m (hereinafter} also referred to as "reference trajectory change amount") is a trajectory change amount y _s by flight direction , it is the sum of the orbit change amount y _w by the wind. And the reference trajectory change amount y _m, by the comparison between the measured trajectory variation [Delta] y, to determine the sign of the angular velocity omega _z of the z-axis.

The trajectory change amount _s according to the flight direction of the y-axis is calculated as follows. First, a measurement starting point p _0, obtains the arbitrary measurement points p _h before track changes, to obtain a change in the y-direction distance y _s' at the measurement point p _h than their difference. At this time, when an extension of the trajectory of the sphere 10 passes through the two points of measurement points p _h and the measurement starting point p _0, track variation y _s by flying directions of the sphere 10 can be calculated by (16) below .

Next, determine the orbit change amount y _w by the wind. When the wind velocity in the y-axis direction and w _y, track variation y _w by the wind can be calculated by (17) below.

Further, the flight time t of the sphere 10 is expressed by using the initial velocity v ₀ , the wind force w _{x in} the x-axis direction, and the orbit change amount Δx in the x-axis direction. _w can be expressed by the following equation (18).

Reference trajectory change amount y _m in the y-axis direction of the spherical body, as described above, the trajectory change amount of the flying direction y _s, and since the sum of orbital variation y _w of wind, of the (16), (18) Based on the formula, it can be expressed by the following formula (19).

Reference trajectory change amount y _m in the y-axis direction of the spherical body 10 that this is calculated as becomes a threshold value for determining the sign of the angular velocity omega _z of z-axis of the sphere 10, and the y _m, actually measured y-axis by comparing the direction of the trajectory change amount [Delta] y, it is possible to determine the direction of the angular velocity omega _z of the z-axis.

If track variation [Delta] y is larger than the reference trajectory change amount _{y m (Δy> y m)} , i.e. if the spheres 10 as shown in FIG. 19 on the left side of the left end of the _{y m,} the positive side in the z-axis (Fig. 19 When viewed from the front side of the sheet, the sphere 10 rotates counterclockwise. Therefore, in the case of [Delta] y> _{y m,} the angular velocity omega _z of the z-axis has a positive direction.

On the other hand, if the trajectory change amount [Delta] y is the reference trajectory change amount y _m is smaller than ([Delta] y <y _m), i.e. if the sphere 10 in FIG 19 is to the right from the left end of the y _m, when viewed from the positive side in the z-axis The sphere 10 is rotating clockwise. Therefore, in the case of [Delta] y _{<y m,} the angular velocity omega _z of the z-axis is negative direction. Thus, the positive and negative angular velocity omega _z of the z-axis can be determined by (20) below.

Note that the sphere 10 tends to change the trajectory in the same way regardless of whether the angular velocity ω _x of the x axis, which is the axis of travel, is positive or negative. Therefore, in this embodiment, the angular velocity direction determination unit 65 selects an arbitrary direction for the x-axis angular velocity ω _x . When the process of step S105 is completed, the process proceeds to step S106.

ここで、ｉ，ｊ，ｋは、それぞれｘ，ｙ，ｚ軸の単位ベクトルである。角速度ベクトルの概略を図２０に示す。図２０に示す角速度ベクトルωのベクトル量（長さ）が球体１０の回転速度、角速度ベクトルωの向きが球体１０の回転軸方向となる。ステップＳ１０６の処理が完了すると、本制御フローを終了する。 (Derivation of angular velocity vector)
In step S106, the angular velocity vector deriving unit 66 outputs the angular velocity vector ω of the sphere 10 using the angular velocities ω _x , ω _y , and ω _z of the respective axes calculated in steps S104 and S105. The angular velocity vector ω can be expressed by the following equation (21).

Here, i, j, and k are unit vectors of the x, y, and z axes, respectively. An outline of the angular velocity vector is shown in FIG. The vector amount (length) of the angular velocity vector ω shown in FIG. 20 is the rotational speed of the sphere 10, and the direction of the angular velocity vector ω is the rotational axis direction of the sphere 10. When the process of step S106 is completed, this control flow ends.

(Use of angular velocity vector)
In addition, after deriving the angular velocity vector ω, the calculation device 6 can output the information of the derived angular velocity vector ω to the storage device 7 and the output device 8, store it in the storage device 7, and output it to the output device 8. Further, the rotation speed and the rotation axis direction of the sphere 10 can be displayed on the output device 8 as character information or graphic information.

  Next, the effect of the rotation estimation system 1 according to the present embodiment will be described.

The rotation estimation system 1 of the present embodiment estimates the rotation speed and rotation axis position of the sphere 10 flying with rotation. The rotation estimation system 1 irradiates the sphere 10 with the electromagnetic waves 11, 13, and 15 from the x, y, and z directions, and receives the reflected electromagnetic waves 12, 14, and 16 from the sphere 10 with the x, y, and z axes. Then, the measuring device 5 that measures the time series information regarding the reception intensity of the reflected electromagnetic waves of the x, y, and z axes, and the time series information of the x, y, and z axes measured by the measuring device 5 are converted into frequency information, respectively. A frequency domain conversion unit 61 of the measuring device 6 that performs measurement, and a threshold processing calculation unit 62 that specifies peak intervals Δf _x , Δf _y , Δf _z in the x, y, and z axis frequency information converted by the frequency domain conversion unit 61. And the angular velocity separation for estimating the angular velocities ω _x , ω _y , and ω _z of the sphere 10 based on the identified peak intervals Δf _x , Δf _y , and Δf _z of the x, y, and z axes. 64 and the angular velocity separation unit 64 And x, y, the angular velocity omega _x of _z-axis, as the rotational speed of the omega _y, omega sphere 10 the magnitude of the angular velocity vector omega _{of z} to a component, also the direction of the angular velocity vector omega as a rotation axis position of the sphere 10 And an angular velocity vector deriving unit 66 for outputting.

  With this configuration, information regarding the rotational phase change of the sphere 10 flying with rotation (time-series information regarding the received intensity of the reflected electromagnetic wave) can be acquired from three directions in a multifaceted manner. It can be grasped with high accuracy. In addition, since the rotational speed and the rotational axis position of the sphere 10 are estimated based on these various pieces of information, it is possible to estimate a rotation parameter that preferably conforms to the actual rotational behavior of the sphere 10 and to rotate the sphere 10. Parameters can be estimated with high accuracy.

In the rotation estimation system 1 of the present embodiment, the angular velocity separation unit 64 serves as a viewpoint when the sphere 10 is viewed from each axial direction with respect to the triaxial peak intervals Δf _x , Δf _y , and Δf _z. The three-axis angular velocities ω _x , ω _y , and ω _z of the sphere 10 are estimated based on the characteristic that the information corresponds to the rotational speed obtained by synthesizing the angular velocities of rotation around the other two axes excluding the axes. More specifically, the angular velocity separation unit 64 sets the peak intervals of the x, y, and z axes to Δf _x , Δf _y , Δf _z , x, y, and z axis respectively, and the angular velocities to ω _x , ω _y , and ω _z , respectively. When expressed, the three-axis angular velocities are estimated using the above equations (5) to (7).

As described above, the frequency spectrum of the received intensity on the x, y, and z axes of the rotating sphere 10 has a characteristic of having periodic peaks such that the peak interval coincides with the rotation speed on each axis. Therefore, by using the above equations (5) to (7) using the peak intervals Δf _x , Δf _y , Δf _z of the frequency spectrum of the received intensity on the x, y, z axes, the x, y, The z-axis angular velocity can be accurately estimated.

Further, the rotation estimation system 1 of the present embodiment acquires the trajectory changes Δx, Δy, Δz of the sphere 10 from the x, y, and z axis directions, and the reference trajectory change x expected from the flight direction of the sphere 10. An angular velocity direction determination unit 65 that determines the positive / negative direction of the angular velocities ω _x , ω _y , ω _z estimated by the angular velocity separation unit 64 by comparing with _{m 1} , y _m , z _m is provided.

As described above, the angular velocities ω _x , ω _y , and ω _z calculated by the angular velocity separation unit 64 using the above expressions (5) to (7) include both positive and negative and remain ambiguous. It is. Therefore, the ambiguity of the angular velocity calculated by the angular velocity separation unit 64 can be resolved by specifying the direction of the angular velocity based on the trajectory change of the sphere 10 with the above configuration. Thereby, the estimation accuracy of the rotation parameter of the sphere 10 finally derived can be improved.

Further, in the rotation estimation system 1 of the present embodiment, the angular velocity direction determination unit 65 uses P ₀ (x ₀ , y ₀ , z ₀ ) as the measurement start point of the reflected electromagnetic wave from the sphere 10 and performs arbitrary measurement before the trajectory change. The point is P _h (x _h , y _h , z _h ), the measurement end point after the trajectory change is p _i (x _i , y _i , z _i ), x from the measurement start point p ₀ to the measurement end point p _i , Y, z-axis trajectory change amounts are represented by Δx, Δy, Δz, gravitational acceleration g, initial velocity v ₀ , x, y, z-axis wind speeds w _x , w _y , w _z The positive / negative direction of the angular velocity is determined by equations (15) and (20).

When the rotation direction of the sphere 10 is positive or negative, the trajectory change of the sphere 10 due to this rotation is clearly different. With the above configuration, the trajectory change amount generated by the influence of the rotation of the sphere 10 can be cut out by offsetting the influence of the trajectory change due to the flying direction of the sphere, gravity drop, and wind force. Therefore, the positive and negative directions of the angular velocities ω _x , ω _y , and ω _z of the sphere 10 can be accurately determined.

Further, in the rotation estimation system 1 of the present embodiment, the threshold processing calculation unit 62 extracts peaks that are equal to or greater than a predetermined threshold from the frequency information of the x, y, and z axes converted by the frequency domain conversion unit 61, and adjacent The peak interval between the peaks to be calculated is calculated, the mode value processing calculation unit 63 extracts a high frequency among the peak intervals calculated by the threshold value processing calculation unit 62, and the peak interval is calculated based on the extracted peak interval. Δf _x , Δf _y , and Δf _z are specified.

The frequency spectrum of the received intensity on the x, y, and z axes for the sphere 10 that rotates and has a periodic peak whose peak interval matches the rotational speed on each axis. However, in this frequency spectrum, for example, when a peak unrelated to the rotation of the sphere 10 occurs due to the influence of noise or the like when measuring the reflected electromagnetic wave, or when the peak of the originally generated frequency cannot be obtained, May contain uncertainties. Such uncertainties may adversely affect the final rotation parameter estimation. Therefore, with the above-described configuration, it is possible to exclude uncertain elements from peak interval information by performing mode value processing and specifying peak intervals Δf _x , Δf _y , Δf _z . This makes it possible to specify the correct peak interval related to the rotation of the sphere 10 and improve the estimation accuracy of the rotation parameter.

  The rotation estimation system 1 of this embodiment can be applied to ball games. For example, if the sphere 10 to be measured is a ball game ball, the ball type of the ball is estimated from the rotation axis position based on the rotation speed and rotation axis position of the ball game ball output by the angular velocity vector deriving unit 66. It is possible to generate information regarding the ball type and trajectory of the ball game ball that rotates and fly, such as estimating the degree of change of the changing ball from the rotation speed and predicting the start of the ball. And the information produced | generated in this way can also be displayed on the display apparatus of the output device 8 or a ball game hall, and can also be shown to a player or an audience.

  The effect of applying the rotation estimation system 1 of the present embodiment to a ball game will be further described. The effect when this embodiment is applied to a ball game will be described on both the player side and the audience side.

  First, the effect on the audience side will be explained. The ball used for the game has a predetermined standard, and it is difficult to measure the rotation parameter in real time during the game because the surface processing is particularly difficult. However, the rotation estimation system 1 of the present embodiment is an electromagnetic wave. Is used to measure the angular velocity vector from the change in received intensity, so that it is possible to measure the rotation parameter in real time even with an unprocessed sphere.

  In this way, if the information such as the ball rotation speed and rotation axis position measured in real time can be displayed on the screen of the venue and provided information to the audience, what kind of ball the player threw or hit This makes it possible to provide services such as improving the sense of reality and estimating tactics.

  As an effect on the player side, since the rotation parameter of the ball thrown or hit during the game can be acquired in real time, the technical correction on the spot can be performed. Also, since it is possible to leave the parameters of the ball thrown during the game as data, analysis of the vague definitions such as “heavy ball” and “ball with elongation”, improvement of weak points using data, It is possible to improve skills such as training proposals and confirming the condition by comparing parameters during good and bad conditions, which can lead to further human resource development.

  A conventional method (see Patent Document 1) for deriving the rotational speed of a sphere using a Doppler shift generated by the rotation of the sphere is known. The Doppler shift is a phenomenon in which the frequency of the reflected wave changes due to the relative velocity between the wave generation source (electromagnetic wave transmitter) and the measurement object (sphere). In this conventional method, the cause of Doppler shift is the rotational speed of the sphere. On the other hand, the rotation estimation system 1 of the present embodiment calculates the rotation speed and the rotation axis of the sphere 10 based on the temporal change of the reception intensity generated by the change in the radar cross-sectional area of the sphere 10. The radar cross-sectional area of the sphere 10 used as input information in the present embodiment does not affect the Doppler shift, and is information that cannot be acquired by the conventional method described above. As described above, the rotation estimation system 1 according to the present embodiment is different in the input information used for derivation of the sphere parameter from the conventional method using the Doppler shift. The analysis method for deriving the parameters is also different.

  Next, the present invention will be described more specifically based on examples. FIG. 21 shows a configuration of an example in which a baseball curve ball is measured from three directions. As the sphere 10 to be measured, a “game ball” described in Japanese Patent No. 5142366 is used. The transmission power is 0 (dBm), the transmission frequency is 10 (GHz), the measurement time is 0.4096 (seconds), the measurement point is 8192 (points), the frequency resolution is 2.44 (Hz), and the x-axis radar 2 is The back net, y-axis radar 3 was placed on the 1st base, and z-axis radar 4 was placed on the ceiling. The distances between the initial position (pitcher mound) of the sphere 10 and the x-axis radar 2, the y-axis radar 3, and the z-axis radar were 36.4 (m), 40 (m), and 50 (m), respectively. . The rotational speed of each axis of the sphere 10 was all 20 (rps). The wind is not blowing.

  FIGS. 22, 23, and 24 show temporal changes in the received intensity related to the reflected electromagnetic waves from the sphere 10 obtained under the above conditions. 22, 23, and 24 are diagrams showing temporal changes in the received strengths of the acquired x, y, and z axes, respectively. 22-24, the vertical axis | shaft shows receiving strength (pW) and a horizontal axis shows measurement time (s (second)).

  When the time change of the received intensity of each axis obtained at this time is Fourier transformed, the frequency spectrum shown in FIGS. 25, 26 and 27 can be obtained. 25, 26, and 27 are diagrams showing frequency spectra of received intensities on the x, y, and z axes, respectively. 25 to 27, the vertical axis represents the power spectrum (dB), and the horizontal axis represents the frequency (Hz / rps). In these frequency spectra, the maximum value excluding the DC component peak is when f = 80.57 (Hz / rps) for each axis. Therefore, when the threshold value is obtained from the average value in this range, the threshold value is −123.03 (dB) in the frequency spectrum of the x axis shown in FIG. In the case of the y-axis frequency spectrum shown in FIG. 26, the threshold value was −128.98 (dB). In the case of the z-axis frequency spectrum shown in FIG. 27, the threshold value was −132.72 (dB). A peak and its interval are obtained by performing threshold processing using the threshold derived in this way. The number of peak intervals obtained at this time is 18 in the x-axis frequency spectrum shown in FIG. 25, 21 in the y-axis frequency spectrum shown in FIG. 26, and 20 in the z-axis frequency spectrum shown in FIG. It was.

  Next, the peak interval was grouped. The minimum peak intervals in the frequency spectrum of each axis shown in FIGS. 25, 26, and 27 were all 19.53 (Hz / rps). A value up to 24.41 (Hz / rps) obtained by adding a value obtained by doubling the frequency resolution 2.44 (Hz) from this value is group 1. Next, the minimum peak interval excluding group 1 was determined, and the same processing was performed. This process was performed until all values belonged to any group.

  As a result of this processing, in the case of the x-axis and z-axis frequency spectra shown in FIGS. 25 and 27, the group is 19.53 (Hz / rps) ≦ 24.41 (Hz / rps) group 1, 58.59 ( Hz / rps) ≦ 63.47 (Hz / rps) group 2 and 78.13 (Hz / rps) ≦ 83.01 (Hz / rps) group 3. In the case of the y-axis frequency spectrum shown in FIG. 26, the group is 19.53 (Hz / rps) ≦ 24.41 (Hz / rps) group 1 and 41.50 (Hz / rps) ≦ 46.38. (Hz / rps) group 2, 61.04 (Hz / rps) ≦ 65.92 (Hz / rps) group 3, 78.13 (Hz / rps) ≦ 83.01 (Hz / rps) group 4 It became.

Since grouping is completed, the number of elements in the group is counted, and the group with the largest number and the ratio in the whole are confirmed. In the case of the x-axis frequency spectrum shown in FIG. 25, there are nine groups 3 including the most elements. Since the total number of elements is 18, the ratio of the number of elements in the group to the total number of elements is 50%. Since elements in the group accounts for more than half of the total, and acquires the average value of the elements in this group, make it a peak interval Delta] f _x of x-axis. Since the average value of Group 3 is 80.02 (Hz / rps9), Δf _x = 80.02 (Hz / rps).

In the case of the y-axis frequency spectrum shown in FIG. 26, the groups having the largest number of elements are seven groups 1 and 4, and the ratio of the number of elements in the group to the total number of elements is 47.62%. In the case of the z-axis frequency spectrum shown in FIG. 27, there are seven groups 2 and 3 that contain the most elements, and the ratio of the number of elements in the group to the total number of elements is 35%. Since the ratio of the number of elements in the group does not exceed 50%, the peak f _i corresponding to the number i of the peak interval Δf _i other than these groups is deleted. In the case of the y-axis frequency spectrum shown in FIG. 26, the peak intervals belonging to groups 2 and 3 are Δf ₂ , Δf ₄ , Δf ₁₀ , Δf ₁₂ , Δf ₁₅ , Δf ₁₇ , Δf ₁₉ , Corresponding peaks f ₂ , f ₄ , f ₁₀ , f ₁₂ , f ₁₅ , f ₁₇ , f ₁₉ were deleted. Similarly, in the case of the z-axis frequency spectrum shown in FIG. 27, the peak intervals belonging to group 1 are Δf ₂ , Δf ₁₀ , Δf ₁₂ , Δf ₁₃ , Δf ₁₅ , Δf ₁₇ , and therefore correspond to these numbers. deleting the peak _{f 2, f 10, f 12} , f 13, f 15, f 17.

Next, the peak number i and the peak interval Δf _i were obtained again in the y-axis and z-axis frequency spectra from which the peaks were deleted. FIG. 28 and FIG. 29 are diagrams illustrating frequency spectra of received intensity on the y-axis and the z-axis obtained by reacquiring the peak interval. 28 and 29, the vertical axis represents the power spectrum (dB), and the horizontal axis represents the frequency (Hz / rps).

  In the frequency spectrum of the received intensity on the y-axis and z-axis shown in FIGS. 28 and 29, the peak intervals were grouped again. In the case of the y-axis frequency spectrum shown in FIG. 28, the total number of elements of the peak interval is 15, and the group is 19.53 (Hz / rps) ≦ 24.41 (Hz / rps) group 1, 61.04 ( (Hz / rps) ≦ 65.92 (Hz / rps) group 2, 78.13 (Hz / rps) ≦ 83.01 (Hz / rps) group 3, 139.16 (Hz / rps) ≦ 144.01 It became group 4 of (Hz / rps). In the case of the z-axis frequency spectrum shown in FIG. 29, the total number of elements of the peak interval is 15, and the group is 58.59 (Hz / rps) ≦ 63.47 (Hz / rps) group 1, 78.13 ( Hz / rps) ≦ 83.01 (Hz / rps) group 2 and 100.10 (Hz / rps) ≦ 104.98 (Hz / rps) group 3.

Since grouping is completed, the number of elements in the group is counted, and the group with the largest number and the ratio in the whole are confirmed. In the case of the frequency spectrum of the y-axis shown in FIG. 28, the group containing the largest number of elements is 11 in group 3, whereas the total number of elements is 15. Therefore, the ratio of the number of elements in the group to the total number of elements is 73. %. In the z-axis frequency spectrum shown in FIG. 29, the group containing the largest number of elements is 11 in group 2, whereas the total number of elements is 15. Therefore, the ratio of the number of elements in the group to the total number of elements is 73. %. Since both of these groups occupy more than half of the total number of elements, the average value of the elements in this group is obtained, and these are used as the peak intervals Δf _y and Δf _{z of} the y-axis and z-axis. The average value of group 3 in FIG. 28 is 79.59 (Hz / rps), and the average value of group 2 in FIG. 29 is 79.59 (Hz / rps), so Δf _y = 79.59 (Hz / rps) Δf _z = 79.59 (Hz / rps).

Since the sphere 10 has four periods when it rotates once, it is necessary to divide Δf _x , Δf _y , and Δf _z by four. As a result, the peak intervals of the frequency spectrum of each axis were Δf _x = 19.90 (Hz / rps), Δf _y = 19.90 (Hz / rps), and Δf _z = 19.90 (Hz / rps).

As the measurement conditions of this example, the rotational speeds of the respective axes of the sphere 10 were all set to 20 (rps) as described above. When the peak intervals Δf _x , Δf _y , Δf _z calculated by the above calculation are compared with a set value (= 20 rps), the estimation error with respect to the set value is suppressed to a very small value of 0.5%.

となった。 When the peak intervals Δf _x , Δf _y , Δf _z obtained as described above are separated into angular velocities ω _x , ω _y , ω _z of the respective axes using the equations (5) to (7), ω _x = 39.79. (Rad / s), ω _y = 39.79π (rad / s), ω _z = 39.79π (rad / s), and the angular velocity vector ω is

It became.

Next, ambiguity was removed based on the trajectory change of the sphere 10. 30, 31, and 32 are diagrams respectively showing changes in trajectories of the sphere 10 in the x, y, and z axes in the embodiment. 30 to 32, the vertical axis represents the moving distance (m) of the sphere 10, and the horizontal axis represents the measurement time (s (seconds)). As shown in FIGS. 30, 31, and 32, when the coordinates of the measurement start point in this measurement are p ₀ (0, _0, 0), the amount of change in the trajectory of each axis is Δx = 13.56 (m), Δy = 0.25 (m) and Δz = −0.5 (m). The measurement point at this time is _{p 0 = 0 (s),} p h = 0.12 (s), and a p i = 0.4096 (s). Since the trajectory change of the x-axis in the _{p h} is 4 (m), the initial velocity of the sphere becomes 120 (km / h).

z _s , z _f , and z _w in the zx plane coordinates are 0.4068 (m), −0.8125 (m), and 0 (m) from the equations (8), (11), and (13), respectively. It became. From these, the change in the z-axis trajectory (reference trajectory change) when the change due to the rotation of the sphere 10 does not occur is −0.4057 (m) from the equation (14). On the other hand, since Δz = −0.5 (m), ω _y > 0 from the formula (15).

となった。 y _s and y _w in the xy plane coordinates are −0.3465 (m) and 0 (m) from the equations (16) and (18). From these, the change in the trajectory of the y-axis when the change due to the rotation of the sphere 10 does not occur (reference trajectory change amount) is −0.4057 (m) from the equation (19). On the other hand, since Δy = 0.25 (m), ω _z > 0 from Equation (20). From this, the angular velocity vector ω is

It became.

  The above calculation was actually operated 10 times with a computer program, and the maximum value, minimum value, and average value of the calculation time were measured. The CPU of the computer used for this calculation was Core (registered trademark) 2 Duo E7300 (2.66 (GHz)) manufactured by Intel (registered trademark), and the memory was 2 (GB). As a programming language, Matlab (registered trademark) was used to create a computer program for this calculation. The calculation time using the computer and the program language is a maximum value of 0.687 seconds, a minimum value of 0.578 seconds, and an average value of 0.631 seconds. In any case, the angular velocity vector can be calculated in 1 second or less. It was. Since these calculation times depend on the machine specifications, the program language, and how to write the program, the calculation time can be further shortened by changing the above conditions.

  As mentioned above, although embodiment of this invention was described, the said embodiment was shown as an example and is not intending limiting the range of invention. The above-described embodiment can be implemented in various other forms, and various omissions, replacements, and changes can be made without departing from the spirit of the invention. The above-described embodiments and modifications thereof are included in the invention described in the claims and equivalents thereof, as long as they are included in the scope and gist of the invention.

  In the above-described embodiment, the sphere 10 has been described as an example of the rotating flying object that is the target for which the rotation estimation system 1 according to the present invention estimates the rotation parameter. However, the rotating flying object only needs to be an object that rotates and flies. For example, an object having a shape other than a sphere, such as an ellipsoid, can be included in the rotating flying object. In the above embodiment, a ball used mainly for a ball game such as a baseball is assumed as the rotating flying object. However, for example, an object other than the ball game such as a cannonball can be used as a target.

1 rotation estimation system 5 measuring device (measuring means)
6 Calculator 61 Frequency domain converter (conversion means)
62 Threshold processing calculation unit (peak interval specifying means)
63 Mode value processing calculation section (peak interval specifying means)
64 Angular velocity separation unit (angular velocity estimation means)
65 Angular velocity direction determination unit (determination means)
66 Angular velocity vector deriving section (output means)
8 Output device (display means)
10 Sphere (Rotating flying object)

Claims (9)

A rotation estimation system for estimating a rotation speed and a rotation axis position of a rotating projectile that flies with rotation,
Irradiate the rotating projectile with electromagnetic waves from three directions, receive the reflected electromagnetic waves from the rotating projectile with the three axes, and measure time-series information regarding the reception intensity of the three-axis reflected electromagnetic waves. Measuring means;
Conversion means for converting the time series information of the three axes measured by the measurement means into frequency information, and
Peak interval specifying means for specifying each of the peak intervals in the triaxial frequency information converted by the converting means;
Angular velocity estimating means for estimating an angular velocity of the three axes of the rotating projectile based on the peak distance of the three axes specified by the peak interval specifying means;
The magnitude of the angular velocity vector whose component is the angular velocity of the three axes estimated by the angular velocity estimating means is output as the rotational speed of the rotating projectile, and the direction of the angular velocity vector is output as the rotational axis position of the rotating projectile. Output means for
A rotation estimation system for a rotating projectile, comprising: The angular velocity estimation means includes
Each of the three-axis peak intervals is information corresponding to a rotational speed obtained by synthesizing angular speeds of rotation around the other two axes excluding the viewpoint axis when the rotating projectile is viewed from each axial direction. The rotational estimation system for a rotating projectile according to claim 1, wherein the three-axis angular velocity of the rotating projectile is estimated based on the above. The angular velocity estimating means represents the three-axis peak intervals by Δf _x , Δf _y , Δf _z , and the three-axis angular velocities by ω _x , ω _y , ω _z , respectively. The rotation estimation system for a rotating projectile according to claim 2, wherein an angular velocity of the shaft is estimated.

The trajectory change of the rotating projectile is obtained from the three axis directions, and compared with the trajectory change expected from the flight direction of the rotating projectile, the positive / negative direction of the angular velocity estimated by the angular velocity estimating means is obtained. The rotation estimation system for a rotating projectile according to any one of claims 1 to 3, further comprising determination means for determining. The determination means includes
The measurement start point of the reflected electromagnetic wave from the rotating projectile is P ₀ (x ₀ , y ₀ , z ₀ ), the arbitrary measurement point before the trajectory change is P _h (x _h , y _h , z _h ), and the trajectory change. The subsequent measurement end point is p _i (x _i , y _i , z _i ), the three-axis trajectory variation from the measurement start point p ₀ to the measurement end point p _i is Δx, Δy, Δz, and gravitational acceleration. ₅ , wherein the initial velocity is represented by v ₀ , and the three-axis wind velocity is represented by w _x , w _y , w _z , the positive / negative direction of the angular velocity is determined by the following discriminant: The rotation estimation system of the described rotational projectile.

The peak interval specifying means is
Extracting a peak above a predetermined threshold from the triaxial frequency information converted by the conversion means, calculating a peak interval between adjacent peaks,
The rotation flight according to any one of claims 1 to 5, wherein a high frequency of the calculated peak intervals is extracted, and the peak interval is specified based on the extracted peak intervals. Body rotation estimation system. A display means,
The rotating projectile is a ball for ball games,
The information on the ball type or trajectory of the ball game ball is displayed on the display unit based on the rotation speed and the rotation axis position of the ball game ball output by the output unit. The rotation estimation system for a rotating flying object according to any one of claims 6 to 6. A rotation estimation method for estimating a rotation speed and a rotation axis position of a rotating projectile that flies with rotation,
Irradiate the rotating projectile with electromagnetic waves from three directions, receive the reflected electromagnetic waves from the rotating projectile with the three axes, and measure time-series information regarding the reception intensity of the three-axis reflected electromagnetic waves. A measurement step;
A conversion step for converting the three-axis time series information measured by the measurement step into frequency information, respectively.
A peak interval specifying step for specifying each of the peak intervals in the triaxial frequency information converted by the converting step;
An angular velocity estimating step for estimating an angular velocity of the three axes of the rotating projectile based on the peak interval of the three axes specified by the peak interval specifying step;
The magnitude of the angular velocity vector whose component is the angular velocity of the three axes estimated by the angular velocity estimating step is output as the rotational speed of the rotating projectile, and the direction of the angular velocity vector is output as the rotational axis position of the rotating projectile. An output step to
A rotation estimation method for a rotating projectile, comprising: A rotation estimation program for estimating the rotation speed and rotation axis position of a rotating projectile that flies with rotation,
Irradiate the rotating projectile with electromagnetic waves from three directions, receive the reflected electromagnetic waves from the rotating projectile with the three axes, and measure time-series information regarding the reception intensity of the three-axis reflected electromagnetic waves. Measurement function,
A conversion function for converting the three-axis time series information measured by the measurement function into frequency information, and
A peak interval specifying function for specifying each of the peak intervals in the triaxial frequency information converted by the conversion function;
An angular velocity estimation function for estimating an angular velocity of the three axes of the rotating projectile based on the peak interval of the three axes specified by the peak interval specifying function;
The magnitude of the angular velocity vector whose component is the angular velocity of the three axes estimated by the angular velocity estimation function is output as the rotational speed of the rotating projectile, and the direction of the angular velocity vector is output as the rotational axis position of the rotating projectile. Output function to
Rotational projectile rotation estimation program for realizing a computer.

Images