Год назад мной была описана метеостанция с измерителем направления и скорости ветра. По опыту двух сезонов эксплуатации в нее были внесены некоторые изменения и улучшения, которые частично описаны там же в дополнениях к основному тексту.

Одно из таких изменений касается расчета среднего значения направления скорости ветра. К моему изумлению, в Сети я ничего по этой теме толкового не нашел, только на одном из форумов народ самостоятельно почти додумался до метода векторного осреднения (но только почти, задачу в общем виде и там не решили). Потому я посчитал полезным вынести эту тему в отдельную публикацию — вдруг кому еще пригодится. Метод может использоваться для осреднения любых векторных величин, не только ветра или течения.

Заметим, что каноническая метода проведения метеорологических измерений ветра такова: усреднять вектор (т.е. и скорость и направление) в течение 10 минут (в России и большинстве стран мира принято именно так, говорят, в США и некоторых других странах иначе). При этом измерения проводятся на высоте 10 м от поверхности Земли. Обеспечить в наколеночных условиях все это довольно затруднительно: и высоты 10 м не достичь на открытом пространстве (не строить же специальную вышку, причем вдали от домов и деревьев, искажающих поток ветра), и температуру с влажностью нужно, наоборот, мерить в тени и вблизи от поверхности. Из компактного прибора выносной датчик превращается в целую измерительную систему (см. фото территории метеостанции в городе Кирове).



И полученный результат — среднее за 10 минут — будет мало информативен. Вблизи поверхности ветер гораздо более порывистый, чем на высоте, за 10 минут может раз двадцать сменить скорость и направление, и информация об этих порывах для обывателя куда более информативна, чем средняя величина. Напомню, что датчик скорости у нас демонстрирует максимальную величину из четырех измерений за цикл 8 с, и это оказалось правильным выбором (по сути мы вместо датчика средней скорости получили датчик пульсаций).

Но вот флюгер оказался более капризным, чем датчик скорости. По изначальному алгоритму моей метеостанции (который выбирался исходя из максимально возможного энергосбережения), направление измерялось раз в цикл, то есть даже пульсаций не получалось: были случайные выборки из непрерывного процесса болтания флюгера туда-сюда с частотой намного большей, чем раз в 8 и тем более 16 секунд.

Поэтому было принято решение осреднить направление вектора скорости за цикл, проводя измерения раз в две секунды и вычисляя среднее. А дело это не вот так, чтобы с полпинка решается — значения направления не образуют равномерного массива чисел, которые можно напрямую сложить и поделить (одно слово — вектор, а не фуфло скалярное). Обычно приводится пример со значениями 1 градус и 359 градусов: легко сообразить, что в среднем это будет ровно 360 (или 0, без разницы), но обычная арифметика выдаст число 180 градусов.

Здесь не надо ничего изобретать — все уже изобретено до вас. Задача решается методом векторного осреднения, хорошо известным тем, кто имел дело с измерениями ветров или течений. Метод по сути очень простой: коли мы не можем напрямую усреднять углы, то давайте усреднять проекции вектора на оси координат, которые, по определению, есть величина скалярная, то есть поддаются обычной арифметике без вопросов.

Проекции текущего вектора ветра W’ (апостроф здесь играет роль надстрочной черты) на оси Х и Y есть wx=Wa•cos(?) и wy=Wa•sin(?), где Wa — модуль вектора (значение скорости), а ? — значение угла между вектором и нулевой осью координат. Если усреднить эти величины проекций, а потом средние преобразовать обратно в вектор, то мы получим истинное значение средней скорости и направления ветра.


Этот замечательный метод (назовем его полным векторным осреднением) имеет один кардинальный недостаток с практической точки зрения: при отсутствии предмета измерений (то есть когда стоит полный штиль, что вполне бытовой случай) он выдает математически некорректный результат: так как скорость ветра равна нулю, то и обе проекции равны нулю, чего быть не может (так как sin и cos взаимодополняющие функции). Точнее, быть-то может, но вот извлечь информацию из такой ситуации принципиально нельзя. А что прикажете демонстрировать на дисплее? Если честно, я так и не знаю корректного способа обойти эту ситуацию (в измерителях течений, которые я когда-то конструировал, циклы осреднения составляли часы, и считалось, что за это время хоть какое-то шевеление произойдет).

Но в нашем случае задача, к счастью, облегчается — скорость нам осреднять не надо, и мы можем обойтись единичными проекциями вектора, без учета величины его модуля. Иными словами, можно оперировать чистыми синусами и косинусами, которые никогда не принимают нулевое значение оба сразу: даже когда ветра нет, застывший в недвижимости флюгер какое-то направление ведь показывает. Назовем такой метод упрощенным векторным осреднением направления (может, у него есть какое-то официальное название, но я не в курсе).

Сложность теперь осталась только одна: превратить вычисленные средние значения проекций обратно в значение угла. Для этого обычно используют функцию ? = arctan(sin(?)/cos(?)), но если уж вычислять через обратные тригонометрические функции, то проще взять просто arcos(cos(?)) (или arcsin(sin(?)), все равно), а дополнять этот результат для получения полного круга (т.е. от 0 до 359 градусов), анализируя знаки проекций, все равно придется в любом случае: все обратные функции выдают результат в пределах полукруга (от 0 до 180 или от -90 до +90). (См. по этому поводу UPD в конце статьи).

Формализуем все сказанное в реальный алгоритм (применительно к Arduino). Для начала будем считывать показания направления не каждый цикл, а каждое измерение (после значения частоты анемометра). Полученный результат в коде Грея (у нас он обозначался, как wind_Gray типа byte, см. ту публикацию) мы преобразуем в обычный двоичный код, и, как и частоту анемометра, поместим его в глобальный массив, который объявим, как wDirAvr [4], где 4 — число измерений в цикле. Преобразование четырехзначного кода Грея в двоичный код расписывать не будем — это можно сделать несколькими способами на усмотрение программиста и описано в любом справочнике.

Этот двоичный код будет принимать значения от 0 до 15, причем договоримся, что углы мы будем отсчитывать, не как сдвинутые географы/топографы/навигаторы, а как нормальные люди, учившие тригонометрию в школе — против часовой стрелки. То есть, если нулевому значению соответствует север, то 90 градусов — не восток, как у «сдвинутых», а запад. Так как градаций направления у нас 16, то для получения направления в обычных градусах дуги значение кода надо умножить на 22,5 (360/16).

Теперь собственно алгоритм упрощенного векторного осреднения направления из 4-х значений кода:

. . . . .
 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 = (wind_Rad*180/M_PI)/22.5; //среднее направление в коде 0-15
. . . . .

Последней строкой мы преобразуем среднее, выраженное в радианах, в среднее, выраженное в нашем коде от 0 до 15. Можно потом преобразовать обратно в код Грея, тогда не придется даже менять программу в основном модуле для отображения направления.

Вот, собственно, и весь алгоритм. Я боялся, что вычисление косинусов-арккосинусов тормознет слабенький (по нынешним 32-разрядным временам) контроллер Arduino, но ничего не произошло: он проглотил код, даже не моргнув… светодиодом, наверное.

UPD: У автора алгоритм работает в данном виде в течение нескольких месяцев без сбоев. Однако комментаторы навели меня на возможную, хотя и крайне маловероятную на практике ошибку: если данные содержат две группы измерений, отстоящих друг от друга примерно на 180 градусов вблизи нулевого значения косинуса (т.е. около 90 и 270 градусов), то алгоритм выдаст ошибочное значение. Для ее избежания вместо acos() следует использовать функцию atan2(wSin,wCos), которая выдает сразу корректный результат с учетом знаков синуса и косинуса (спасибо aamonster за подсказку). Строку, где производится вычисление значения среднего wSin, при этом следует раскомментировать, а строка с поправкой на знак wSin не нужна. Вместо нее надо вставить приведение к положительным значениям угла (так как atan2 выдает значения от pi до -pi):

if (wind_Rad<0) wind_Rad=2*M_PI+wind_Rad;