Vector averaging method for measuring wind direction
A year ago, I described a weather station with a wind speed and direction meter. According to the experience of two seasons of operation, some changes and improvements were made to it, which are partially described there in the supplements to the main text.
One of these changes concerns the calculation of the average value of the wind speed direction. To my amazement, I didn’t find anything sensible on the web, only on one of the forums the people themselves almost thought of the vector averaging method (but only almost, they didn’t solve the problem in general). Therefore, I considered it useful to bring this topic into a separate publication - suddenly, who else would come in handy. The method can be used for averaging any vector quantities, not just wind or current.
Note that the canonical method of conducting meteorological measurements of wind is as follows: to average the vector (ie, speed and direction) for 10 minutes (in Russia and most countries of the world, this is usually the case, they say, in the USA and some other countries otherwise). In this case, measurements are performed at a height of 10 m from the surface of the Earth. It is rather difficult to provide all this in knee-down conditions: you cannot reach 10 m in open space (you should not build a special tower, and far from houses and trees distorting the wind flow), and the temperature with humidity should be measured in the shade and near from the surface. From a compact device, the remote sensor turns into a whole measuring system (see the photo of the territory of the weather station in the city of Kirov).
And the result obtained - an average of 10 minutes - will be of little informative. Near the surface, the wind is much more gusty than at a height, in 10 minutes it can change the speed and direction twenty times, and the information about these gusts for the average person is much more informative than the average value. Let me remind you that the speed sensor we have demonstrates the maximum value of four measurements per cycle of 8 s, and this turned out to be the right choice (in fact, instead of the average speed sensor, we got a pulsation sensor).
But the weather vane turned out to be more capricious than the speed sensor. According to the initial algorithm of my meteorological station (which was chosen on the basis of the maximum possible energy saving), the direction was measured once per cycle, that is, even the pulsations did not work: there were random samples from the continuous process of wind vane rolling back and forth with a frequency much higher than once every 8 and more than 16 seconds.
Therefore, it was decided to average the direction of the velocity vector per cycle, taking measurements every two seconds and calculating the average. But this is not the case in such a way that it resolves with half a ping - the direction values do not form a uniform array of numbers that can be directly folded and divided (one word is a vector, not a garbage). An example is usually given with values of 1 degree and 359 degrees: it is easy to figure out that on average it will be exactly 360 (or 0, no difference), but ordinary arithmetic will produce a number of 180 degrees.
There is no need to invent anything - everything has already been invented before you. The problem is solved by the vector averaging method, well known to those who dealt with measurements of winds or currents. The method is essentially very simple: if we cannot directly average the angles, then let's average the projections of the vector on the coordinate axes, which, by definition, are scalar, that is, amenable to ordinary arithmetic without question.
The projections of the current wind vector W '(the apostrophe here plays the role of a superscript line) on the X and Y axes are wx = Wa • cos (α) and wy = Wa • sin (α), where Wa is the vector modulus (velocity value), and α - the value of the angle between the vector and the zero axis of coordinates. If we average these values of the projections, and then convert the averages back into a vector, we get the true value of the average speed and direction of the wind.
This remarkable method (let's call it full vector averaging) has one fundamental drawback from a practical point of view: in the absence of a measurement object (that is, when there is complete calm, which is quite a household case), it produces a mathematically incorrect result: since the wind speed is zero, and both projections are equal to zero, which cannot be (since sin and cos are complementary functions). More precisely, it can be, but it is fundamentally impossible to extract information from such a situation. And what would you order to show on the display? To be honest, I still do not know the correct way to get around this situation (in the flow meters that I once designed, the averaging cycles made up the clock, and it was believed that during this time at least some kind of movement would occur).
But in our case, the task, fortunately, is easier - we don’t need to average the speed, and we can do with single projections of the vector, without taking into account the magnitude of its module. In other words, it is possible to operate with pure sines and cosines, which never take a zero value both at once: even when there is no wind, a weather vane frozen in real estate shows a certain direction. Let's call such a method a simplified vectorial averaging of a direction (maybe it has some kind of official name, but I am not up to date).
The difficulty now remains only one: turn the calculated average values of the projections back into the angle value. For this, the function α = arctan (sin (α) / cos (α)) is usually used, but if we calculate it via inverse trigonometric functions, then it is easier to take simply arcos (cos (α)) (or arcsin (sin (α)), all the same), and to supplement this result to get a full circle (ie, from 0 to 359 degrees), analyzing the signs of the projections, you still have to in any case: all inverse functions produce a result within a semicircle (from 0 to 180 or from -90 to +90). (See about this UPD at the end of the article).
We formalize everything said in the real algorithm (with reference to Arduino). To begin with, we will read the direction indications not every cycle, but every measurement (after the value of the anemometer frequency). The result obtained in the Gray code (we denoted it as wind_Gray of the byte type, see that publication ) we convert to ordinary binary code, and, like the frequency of the anemometer, put it into the global array, which we declare as wDirAvr [4], where 4 is the number of measurements in the cycle. We will not paint the conversion of the four-digit Gray code to binary code - this can be done in several ways at the discretion of the programmer and is described in any reference book.
This binary code will take values from 0 to 15, and we will agree that we will count the angles, not as shifted geographers / topographers / navigators, but as normal people who have studied trigonometry in school - counterclockwise. That is, if the zero value corresponds to the north, then 90 degrees is not east, as in “shifted” ones, but west. Since we have 16 gradations of direction, then to get directions in normal degrees of arc, the code value must be multiplied by 22.5 (360/16).
Now the algorithm itself of the simplified vectorial averaging of direction from 4 code values:
In the last line, we convert the average, expressed in radians, into the average, expressed in our code from 0 to 15. You can then convert back to Gray code, then you don't even have to change the program in the main module to display the direction.
That's the whole algorithm. I was afraid that the calculation of cosine-arccosinuses would be slowed down by the weak (at the current 32-bit times) Arduino controller, but nothing happened: he swallowed the code without even blinking ... with a LED, I guess.
UPD: The author's algorithm works in this form for several months without a glitch. However, the commentators put me on a possible, though extremely unlikely error in practice: if the data contain two groups of measurements that are approximately 180 degrees apart from each other (ie, about 90 and 270 degrees), then the algorithm will produce an erroneous value. To avoid it, instead of acos (), you should use the atan2 (wSin, wCos) function, which immediately gives the correct result, taking into account the sine and cosine signs (thanks to aamonster for the hint). The line where the mean value of wSin is calculated, should be uncommented, and the line adjusted for the wSin sign is not needed. Instead, you need to insert a cast to positive values of the angle (since atan2 gives values from pi to -pi):
One of these changes concerns the calculation of the average value of the wind speed direction. To my amazement, I didn’t find anything sensible on the web, only on one of the forums the people themselves almost thought of the vector averaging method (but only almost, they didn’t solve the problem in general). Therefore, I considered it useful to bring this topic into a separate publication - suddenly, who else would come in handy. The method can be used for averaging any vector quantities, not just wind or current.
Note that the canonical method of conducting meteorological measurements of wind is as follows: to average the vector (ie, speed and direction) for 10 minutes (in Russia and most countries of the world, this is usually the case, they say, in the USA and some other countries otherwise). In this case, measurements are performed at a height of 10 m from the surface of the Earth. It is rather difficult to provide all this in knee-down conditions: you cannot reach 10 m in open space (you should not build a special tower, and far from houses and trees distorting the wind flow), and the temperature with humidity should be measured in the shade and near from the surface. From a compact device, the remote sensor turns into a whole measuring system (see the photo of the territory of the weather station in the city of Kirov).
And the result obtained - an average of 10 minutes - will be of little informative. Near the surface, the wind is much more gusty than at a height, in 10 minutes it can change the speed and direction twenty times, and the information about these gusts for the average person is much more informative than the average value. Let me remind you that the speed sensor we have demonstrates the maximum value of four measurements per cycle of 8 s, and this turned out to be the right choice (in fact, instead of the average speed sensor, we got a pulsation sensor).
But the weather vane turned out to be more capricious than the speed sensor. According to the initial algorithm of my meteorological station (which was chosen on the basis of the maximum possible energy saving), the direction was measured once per cycle, that is, even the pulsations did not work: there were random samples from the continuous process of wind vane rolling back and forth with a frequency much higher than once every 8 and more than 16 seconds.
Therefore, it was decided to average the direction of the velocity vector per cycle, taking measurements every two seconds and calculating the average. But this is not the case in such a way that it resolves with half a ping - the direction values do not form a uniform array of numbers that can be directly folded and divided (one word is a vector, not a garbage). An example is usually given with values of 1 degree and 359 degrees: it is easy to figure out that on average it will be exactly 360 (or 0, no difference), but ordinary arithmetic will produce a number of 180 degrees.
There is no need to invent anything - everything has already been invented before you. The problem is solved by the vector averaging method, well known to those who dealt with measurements of winds or currents. The method is essentially very simple: if we cannot directly average the angles, then let's average the projections of the vector on the coordinate axes, which, by definition, are scalar, that is, amenable to ordinary arithmetic without question.
The projections of the current wind vector W '(the apostrophe here plays the role of a superscript line) on the X and Y axes are wx = Wa • cos (α) and wy = Wa • sin (α), where Wa is the vector modulus (velocity value), and α - the value of the angle between the vector and the zero axis of coordinates. If we average these values of the projections, and then convert the averages back into a vector, we get the true value of the average speed and direction of the wind.
Note for especially corrosive
Для полностью корректного осреднения таким способом надо, чтобы величина Wa (значение скорости) измерялось строго одновременно с величиной угла. На практике за этим нужно следить только в случае, если периоды существенных по амплитуде пульсаций потока меньше или сравнимы с временем проведения измерений. Для ветра (и для практически всех случаев естественного течения воды) обычно следить не требуется, так как время измерений у нас составляет максимум доли секунды, а значимые пульсации ветра, конечно, длятся дольше. Высокочастотными неоднородностями потока мы может пренебречь, потому что они ни на что не влияют: инерция реальных физических тел (в том числе и датчиков) много больше этих неоднородностей, и мы их просто не почувствуем — клочок бумаги будет трепетать, но не более того. В крайнем случае они вылезут, как небольшой случайный шум, не влияющий существенно на качество измерений.
This remarkable method (let's call it full vector averaging) has one fundamental drawback from a practical point of view: in the absence of a measurement object (that is, when there is complete calm, which is quite a household case), it produces a mathematically incorrect result: since the wind speed is zero, and both projections are equal to zero, which cannot be (since sin and cos are complementary functions). More precisely, it can be, but it is fundamentally impossible to extract information from such a situation. And what would you order to show on the display? To be honest, I still do not know the correct way to get around this situation (in the flow meters that I once designed, the averaging cycles made up the clock, and it was believed that during this time at least some kind of movement would occur).
But in our case, the task, fortunately, is easier - we don’t need to average the speed, and we can do with single projections of the vector, without taking into account the magnitude of its module. In other words, it is possible to operate with pure sines and cosines, which never take a zero value both at once: even when there is no wind, a weather vane frozen in real estate shows a certain direction. Let's call such a method a simplified vectorial averaging of a direction (maybe it has some kind of official name, but I am not up to date).
The difficulty now remains only one: turn the calculated average values of the projections back into the angle value. For this, the function α = arctan (sin (α) / cos (α)) is usually used, but if we calculate it via inverse trigonometric functions, then it is easier to take simply arcos (cos (α)) (or arcsin (sin (α)), all the same), and to supplement this result to get a full circle (ie, from 0 to 359 degrees), analyzing the signs of the projections, you still have to in any case: all inverse functions produce a result within a semicircle (from 0 to 180 or from -90 to +90). (See about this UPD at the end of the article).
We formalize everything said in the real algorithm (with reference to Arduino). To begin with, we will read the direction indications not every cycle, but every measurement (after the value of the anemometer frequency). The result obtained in the Gray code (we denoted it as wind_Gray of the byte type, see that publication ) we convert to ordinary binary code, and, like the frequency of the anemometer, put it into the global array, which we declare as wDirAvr [4], where 4 is the number of measurements in the cycle. We will not paint the conversion of the four-digit Gray code to binary code - this can be done in several ways at the discretion of the programmer and is described in any reference book.
This binary code will take values from 0 to 15, and we will agree that we will count the angles, not as shifted geographers / topographers / navigators, but as normal people who have studied trigonometry in school - counterclockwise. That is, if the zero value corresponds to the north, then 90 degrees is not east, as in “shifted” ones, but west. Since we have 16 gradations of direction, then to get directions in normal degrees of arc, the code value must be multiplied by 22.5 (360/16).
Now the algorithm itself of the simplified vectorial averaging of direction from 4 code values:
. . . . .
float wSin=0; //проекция sin
float wCos=0; //проекция cos
float wind_Rad; //направление в радианах
for (byte i=0; i<4; i++){
wind_Rad= ((float(wDirAvr[i])*22.5)*M_PI/180); //направление в радианах
wSin=wSin+sin(wind_Rad);//сумма нормированных проекций вектора sin
wCos=wCos+cos(wind_Rad);//сумма нормированных проекций вектора cos
}
// wSin=wSin/4;//средняя sin – она нам не нужна, потому закомментирована
wCos=wCos/4; //средняя cos
wind_Rad = acos(wCos); //среднее направление в радианах через arccos
if (wSin<0) wind_Rad=2*M_PI-wind_Rad; //поправка на знак sin
int wind_G = round ((wind_Rad*180/M_PI)/22.5); //среднее направление в коде 0-15
. . . . .
In the last line, we convert the average, expressed in radians, into the average, expressed in our code from 0 to 15. You can then convert back to Gray code, then you don't even have to change the program in the main module to display the direction.
That's the whole algorithm. I was afraid that the calculation of cosine-arccosinuses would be slowed down by the weak (at the current 32-bit times) Arduino controller, but nothing happened: he swallowed the code without even blinking ... with a LED, I guess.
UPD: The author's algorithm works in this form for several months without a glitch. However, the commentators put me on a possible, though extremely unlikely error in practice: if the data contain two groups of measurements that are approximately 180 degrees apart from each other (ie, about 90 and 270 degrees), then the algorithm will produce an erroneous value. To avoid it, instead of acos (), you should use the atan2 (wSin, wCos) function, which immediately gives the correct result, taking into account the sine and cosine signs (thanks to aamonster for the hint). The line where the mean value of wSin is calculated, should be uncommented, and the line adjusted for the wSin sign is not needed. Instead, you need to insert a cast to positive values of the angle (since atan2 gives values from pi to -pi):
if (wind_Rad<0) wind_Rad=2*M_PI+wind_Rad;