Friday 7 September 2012

Calculating wind speed from the GPS track

Coding the android app is progressing well and I have most hardware parts for the first ten prototypes. In coding the Android app I got a little fixated on developing an algorithm for estimating wind speed and direction from the gps location flight path. I think what I have done is new. [EDIT: See the comments below from Edgar, I guess what I did was not new]. If it is not new, then as far as I can tell from pretty extensive Google and literature searches, it is undocumented in the public arena. This post gets pretty technical, so first..

The bottom line

Using information from only the GPS, this algorithm provides a pretty good estimate of windspeed, wind direction, and airspeed. The algorithm does not rely on turning circles or airspeed measurements. It does assume that the aircraft is flying at a reasonably constant airspeed over a time period, or at least that the airspeed is not overlay biased towards any heading. It also assumes that the wind speed is pretty constant over the same time period. It does rely on some changes in heading over the time period. If the aircraft heading is distributed over a 60 to 90 degree arc then the windspeed, direction and airspeed will be calculated pretty accurately. By pretty well, I mean less than about 20% error, which is kind of like the variation of windspeed and wind direction caused by gusts in any case.

For a typical paraglider boating around in thermal or on a ridge the assumptions above hold true. Even with a little speed bar our airspeed range is reasonably confined, and most people do not use it while thermalling. A time period of 90 seconds works well from my testing so far. I think that this algorithm would work less well for hang glider and sailplane pilots as they vary airspeed much more often. It will be interesting to get some practical feedback.

My main motivation for working out the windspeed and wind direction is to ensure I can predict where a packet of lift that you have just flown through will have moved to later in time. That magic is for a later post. Now the technical stuff.

What the GPS gives us

Most GPS spit out position information every second. In addition to position we get heading and ground speed. The GPS is not merely taking two location points and calculating the vector between them. Rather, it measures the Doppler shift in the signal from each satellite and converts that to a heading and speed. This is quite nifty and I imagine particularly more accurate than just taking the difference between two points.

The wind triangle

The picture below tells a thousand words. On the left is a 2D flight path in pink. I have marked three points on the path (p1, p2 and p3). The green arrows represent the ground speed and heading vectors at each point. You can also see the blue arrow representing the windspeed and wind direction. Notice that groundspeed at point p2 and p3 is much greater (longer arrow) when the wind is from the side at point p1. On the right side are these three groundspeed vectors plotted on the XY velocity plane, and the windspeed plotted. Note that each of the red airspeed vectors have the same magnitude |a| (our assumption of constant airspeed), and the the windspeed vector is assumed to be the same at each point. The greyed out triangle is the 'wind triangle' for a particular point in time. The most important thing in this picture is the circle. For any three points drawn on a plane (not on a straight line), there is only one circle that fits them all. To work that circle out we only need our three groundspeeds. The centre of the circle gives us the windspeed vector, and its radius is the airspeed.

Practical data

In theory, with just three gps heading and speed readings, provided we were not travelling in a straight line, we could determine windspeed and direction. In practice there are all kinds of errors, our assumptions do not hold exactly, and so on. But with more data we can find the circle that best fits our points. The image below is a graph of all speed and headings from a ~35 min scratchy paragliding flight on a ridge. The units are in m/s.  When I first saw this I was amazed at how consistent both my airspeed and the windspeed was over the 30 minutes. In flight I experienced the range of gusts and direction changes that we are used to when flying a ridge. Yet, it is still easy to see where the circle fits if you were to draw it by hand.

Below is 90 seconds of data from the middle of the flight. This time I have hand drawn a circle on it. The errors of the points are much less.

Circular Regression

To fit a circle to a bunch of points requires the use of advanced math. It is called circular regression. Nikolai Chernov from the University of Alabama had an excellent summary of his work presented here.  He has provided a bunch of algorithms. I chose the Taubin-based Newton method for a number of reasons: Speed (it runs fast enough on Android), I am not worried about numerical stability for the kind of data we will get, it does not need matrix functions, and it is characterised as 'perhaps the best algebraic fit'. The Kasa fit worked well enough but the Taubin was much better at small arcs. In addition to the 'fit' produced I then calculate the mean squared error of the predicted centre and radius of the circle. Some of the methods in JBone code provided inspiration. I will share the actual Java code from the Android app with the initial release in a few weeks time.

Below is an animation of successive circle fits over time. The red dot is the origin. The fits start at 90s. The center pink dot size represents the mean squared error of the circle location. The two circles represent the mean squared error associated with radius.

The graph below is the calculated windspeed and direction over time associated with the above flight. It matches what I experienced.

I have added a simple IIR filter on the windspeed (the yellow line). This smooths out second to second fluctuations and is what will be implemented in code.

Further work

This can be better. The circular regression algorithm tries to minimize the sum of the square errors from all of the points. If it could be altered to weight the points over time, with the most recent having the highest weight, then it might be possible to get a more time sensitive fit. In addition, the current algorithm calculates a full fit with each new point. It might be possible to alter it to reduce the computational intensity by saving some of the computed data for older points, and by using the previous fit as a starting point for the next fit.