Abstract:
A method for measuring wind velocity including using a Fourier Transform (FT) and a correlation between a wave vector number and temporal frequency of a wind to calculate wind velocity. Contour line deviations in a series of images of a far object are evaluated by performing a spatial discrete fourier transform (DFT) on deviations within each image, and subsequently a time discrete fourier transform (DFT) on the Fourier coefficients obtained by the spatial DFT, to get the frequency dependence of each Fourier coefficient.

Description:
RELATED U.S. APPLICATION DATA 
     This application is a 371 of PCT/IL2005/001370 Dec. 22, 2005, which claims benefit of 60/639,307 Dec. 28, 2004. 
    
    
     FIELD OF THE INVENTION 
     This invention relates generally to methods for wind measurement. 
     BACKGROUND OF THE INVENTION 
     The ground temperature is different from the air temperature above it. Heat transfer causes the lower air layer to change its temperature and to differ from layers above it. 
     Mixing of air layers due to several processes (wind, elevation of hot air, etc.) creates streams of air that circulate within the atmosphere in a closed form. These whirlwinds or eddies, first created at large size (L 0 —outer scale,) are not stable and decay to smaller sizes in a cascade process, until they reach a smallest size (l 0 —inner scale) and finally dissipate. The scales L 0  and l 0  depend on local conditions and are hard to predict. The whole phenomenon is called turbulence. 
     The turbulence eddies have a different density than its surround air and therefore somewhat different index of refraction (Hot air is less dense than cold air—at a fixed pressure). For example, during day time, when the ground is hotter than the air, the eddies density is less than that of the surrounding air. The effect of these eddies, with different indices of refraction, on light rays passing through them, is that the eddies act as weak lenses. 
     The appearance of these “lenses” along a line of sight is random and the effect changes with time. Typically many such weak “lenses” are present along a line of sight. 
     Passive methods for wind measurement based on the phenomenon of the wind acting as a refractive lens are known. For example, U.S. Pat. No. 5,469,250 to Holmes, describes a passive wind profilometer that measures wind transverse to a line of sight as a function of range. Schlieren which are natural occurrences in wind act as refracting lenses which move along with the wind. The patches of refracted light are detected by use of a light-intensity-sensitive television. camera and a microprocessor. The range of the schlieren are determined from the scale size of the refracted patches of light. The microprocessor, using a specific algorithm for processing the light intensity data, then statistically compares the positions of the schlieren at different times to determine the wind velocity. 
     However, the Holms patent tends to be inaccurate, inter alia, because it tends to pick up artifacts that are not related to the wind velocity. 
     SUMMARY OF THE INVENTION 
     The present invention seeks to provide a novel passive method for wind measurement, as described hereinbelow. Unlike the prior art, the present invention picks up only frequencies that are related to wind velocity, and eliminates changes due to internal structure changes of the turbulence eddies. In accordance with one non-limiting embodiment of the invention, the method measures wind velocity by using a Fourier Transform (FT) and a correlation between wave vector number and temporal frequency. 
    
    
     DETAILED DESCRIPTION OF EMBODIMENTS 
     Freezed Turbulence Model 
     The turbulent eddies are dragged with the wind, so that their center moves in the wind velocity. These eddies crosses a given line of sight (L.O.S), in a speed equal to the perpendicular wind vector component. Ignoring the internal dynamics of the turbulence, as described above, the typical time a given eddy disturbs the light along the L.O.S is proportional to the eddy size divided by its velocity component perpendicular to the L.O.S. 
     This model, named Freezed Turbulence, is an approximation which usually holds for large eddies (that have long life time) but may fail for the smallest ones. 
     Turbulence Effects on Images 
     Taking a picture of a far object through a turbulent atmosphere, the image is distorted in several ways: 
     1) Intensity fluctuations (Stars blinks due to this effect). 
     2) Position change due to light rays shifting by the turbulence eddies (“lenses”). 
     3) Glare. 
     The second effect (position change) is used here to measure the perpendicular wind velocity component. 
     As a turbulence eddy crosses a camera field of view, it actually crosses many (nearby) lines of sight. We will diagnose two orthogonal cases: 
     1) A horizontal line target (object). 
     2) A vertical line target. 
     It is noted that the eddy is dragged horizontally with the wind. 
     1. The perturbation to the picture moves along the horizontal line when the eddy passes, with a speed proportional to the wind velocity component perpendicular to the L.O.S. 
     2. The perturbation created by the eddy, on a horizontal line, exists at a time period inversely proportional to that wind velocity component. 
     The present invention takes advantage of the relation between the perturbation size and the perturbation speed of change. 
     In order to correlate size (as seen by the camera) to period of time, we use the Fourier transformation (spatial and temporal). As our space is divided into quantized quantities by the hardware (pixels) and the time is quantized by the movie frame rate (time between pictures) we are using a discrete version of this technique. 
     Discrete Fourier Transform (DFT) 
     For length N input vector x, the DFT is a length N vector with elements 
     
       
         
           
             
               
                 j 
                 = 
                 
                   √ 
                   
                     - 
                     1 
                   
                 
               
               ; 
               
                 
                   X 
                   ⁡ 
                   
                     ( 
                     k 
                     ) 
                   
                 
                 = 
                 
                   
                     ∑ 
                     
                       n 
                       = 
                       1 
                     
                     N 
                   
                   ⁢ 
                   
                     
                       x 
                       ⁡ 
                       
                         ( 
                         n 
                         ) 
                       
                     
                     * 
                     
                       exp 
                       ⁡ 
                       
                         ( 
                         
                           
                             - 
                             j 
                           
                           ⁢ 
                           
                               
                           
                           * 
                           2 
                           * 
                           pi 
                           * 
                           
                             ( 
                             
                               k 
                               - 
                               1 
                             
                             ) 
                           
                           * 
                           
                             
                               ( 
                               
                                 n 
                                 - 
                                 1 
                               
                               ) 
                             
                             / 
                             N 
                           
                         
                         ) 
                       
                     
                   
                 
               
             
             , 
             
               
 
             
             ⁢ 
             
               1 
               &lt;= 
               k 
               &lt;= 
               
                 N 
                 . 
               
             
           
         
       
     
     If the x vector is in time space, that is, the indices (n) of the elements stand for discrete points on the time axis, then the indices of the resulted X (k) are points along a so-called frequency axis (usually using f instead of k). 
     If the x vector is in real space, that is, the indices (n) of the elements stand for discrete points along an axis, then the indices of the resulted X (k) are points along a so-called wave-vector axis. 
     The transformed vector elements (Fourier coefficients) are complex numbers with amplitude, A, and phase, φ:
 
X k =A k e jφ     k    
 
     The squares of the amplitudes are called Power Spectrum. 
     In real space, an object draw of length D centered at x 0  will have significant DFT element X(k) at k=2π/D with phase φ k =2 πkx 0 . Other elements for higher k values may appear according to the internal structure of the draw. 
     The number of significant Fourier elements is equal to the number of original data points. If the data points are real numbers, however, the second half of the Fourier coefficient is equal to the first half. 
     Taking a movie of a far object, the contour lines seem to move due to the atmospheric turbulence. On the average, however, they maintain the right position. 
     By registration of the contour line positions, pixel by pixel, in each picture and measure the deviation with respect to the average contour, one may obtain a matrix of numbers wherein each element represents a tangent (to the contour line) deviation at a given contour point at a given time (picture of movie). In other words, each row contains the deviations along the contour in a given picture (time), and each column contains the changes history of a given point on the contour. (Columns and rows can be interchanged of course). 
     For example if the contour is a straight line along the x axis, a row of the matrix includes the deviations in the y direction (caused by turbulence) in one picture of the movie. 
     For the explanation to be clear we can look at any portion of the contour as a straight line (the tangent) segment. This straight segment can be viewed as a composition of horizontal and vertical segments (vector decomposition). 
     The treatment and explanation for horizontal (in the plane containing the wind velocity vector) and for vertical (perpendicular to the wind velocity vector) lines are somewhat different. 
     In both cases the main point is that turbulence is a phenomenon of many closed form air eddies. Each eddy projects a disturbance on the recording imaging system of size proportional to its own size. 
     Before going into details, a mathematical explanation is presented. 
     Spatial DFT of Recorded Data: 
     Taking a row m (t=τm is the time where X is the leg between pictures) in the data matrix of deviations and making a Discrete Fourier Transform operation on it, one may get some significant Fourier elements, {X(t,k)} in accordance with the size of the disturbances. Each element can be separated into amplitude and phase:
 
 X ( m,k )= A   k ( t ) e   f2πkx     0     (l)   ; t=mτ 
 
     For a single eddy disturbance, we can clearly separate two cases: 
     1) A vertical line—the disturbance starts and ends at same height so that x 0  is fixed. What is changed during time is the amplitude, A, that grows up and than grows down until it disappear when the eddy goes away. The most significant amplitude is for k=2π/D where D is proportional to the eddy size. 
     2) A horizontal line—the disturbance stays for long time changing its position along the line: x 0 =V x t 
     wherein V x  is the velocity of deformation motion (proportional to the velocity that the eddy crosses the field of view of the camera). So that:
 
φ k =2πkVt
 
The most significant (with highest amplitude) k is again the one with k=2π/D.
 
     Time DFT 
     After operating the DFT process for all rows, we obtain a matrix with elements X(m,k), where m=t/τ is the time index. 
     We now perform another DFT along columns. This gives use the frequency dependence of each k Fourier coefficient. For each coefficient there is a typical frequency that measures its rate of change. One way to calculate the average frequency, which is used here, is according to a weighted average: 
     
       
         
           
             fk 
             = 
             
               
                 
                   ∑ 
                   
                     f 
                     = 
                     0 
                   
                   
                     f 
                     
                       ma 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       x 
                     
                   
                 
                 ⁢ 
                 
                   fP 
                   ⁡ 
                   
                     ( 
                     
                       f 
                       , 
                       k 
                     
                     ) 
                   
                 
               
               
                 
                   ∑ 
                   
                     f 
                     = 
                     0 
                   
                   
                     f 
                     
                       ma 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       x 
                     
                   
                 
                 ⁢ 
                 
                   P 
                   ⁡ 
                   
                     ( 
                     
                       f 
                       , 
                       k 
                     
                     ) 
                   
                 
               
             
           
         
       
     
     Here P=|X| 2  is the power spectrum elements. 
     1) Horizontal line—for the simple case studied above for this kind of target the time dependence is only in the phase and the DFT yields a Kroneker delta function:
 
δ(2πf−2πkV)
 
     so that only f k =kV Fourier element should present for each k. In real life there are some other contributions but the averaging process will give the right value. 
     2) For the horizontal line case the typical frequency is inversely proportional to the time this element exists, that is the disturbance with size of length D y =2π/k 
     This time depends on eddy size and velocity so that: 
     
       
         
           
             
               
                 T 
                 
                   l 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   n 
                 
               
               = 
               
                 
                   D 
                   x 
                 
                 / 
                 
                   V 
                   x 
                 
               
             
             ; 
             
               f 
               = 
               
                 
                   
                     2 
                     ⁢ 
                     π 
                   
                   
                     T 
                     
                       i 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       n 
                     
                   
                 
                 = 
                 
                   
                     2 
                     ⁢ 
                     π 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       V 
                       x 
                     
                   
                   
                     D 
                     x 
                   
                 
               
             
           
         
       
     
     The fact that we are dealing with a circulating air stream enforces the relation:
 
D x ≈D y  
 
     So putting D y  instead of D x  and further express it using k one gets:
 
f k =kV x  
 
     Again this is an approximated way to estimate the expected average frequency while the exact proof is more complicated. 
     So for both vertical and horizontal lines the velocity is given by: 
     
       
         
           
             
               For 
               ⁢ 
               
                   
               
               ⁢ 
               any 
               ⁢ 
               
                   
               
               ⁢ 
               
                 k 
                 · 
                 
                   V 
                   x 
                 
               
             
             ∝ 
             
               
                 〈 
                 fk 
                 〉 
               
               k 
             
           
         
       
     
     The algorithm further takes the average of the results for V x  from several values of k. 
     The k values in use of the algorithm were picked according to a “life time” of turbulent eddies, according to their size and other factors, such as target line length, compared to expected outer scale (L 0 ) of the eddies. 
     Some applications of the invention include, without limitation: 
     Fire control system—to correct for wind drag of bullets, shells and rockets. 
     Movies correction—using the k−ω relations to correct movies pictures that are distorted by turbulence. 
     Simulation of turbulence effects on movies—Taking a “clean” movie and adding turbulence effects to it for using in simulators.