Метод векторного осреднения при измерении направления ветра

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

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

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

    Метеостанция в Кирове

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

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

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

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

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

    Примечание для особо въедливых
    Для полностью корректного осреднения таким способом надо, чтобы величина 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 = round ((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;
    
    Поделиться публикацией
    Комментарии 19
      +1

      Вы забыли отнормировать результат осреднения (wCos /= hypot(wCos, wSin)).
      Но лучше просто использовать atan2(wSin, wCos) — тогда не надо ни нормировать, ни разбирать варианты по знаку wSin.

        0
        Спасибо, однако. Про atan2 я без вас бы и не вспомнил, конечно, его удобнее применять.
          0
          Было бы интересно увидеть графически как было «случайные выборки из непрерывного процесса болтания флюгера туда-сюда с частотой намного большей, чем раз в 8 и тем более 16 секунд» и как стало.
            0
            Да не было это, конечно, зафиксировано, нет такой нужды что-то тут точно мерить, так что оценка на глаз. Разброс показаний (любых) вообще очень хорошо на глаз оценивается — если разница очень заметна и необъяснима, то все плохо, надо сглаживать. На примере направления — был север, через 8 секунд стал юго-запад, потом северо-восток и сразу запад — уже видно, что никуда не годится. Потому что на этих данных заключения о том, куда же дует ветер, сделать нельзя. После исправления подобные метания сократились до пределов примерно одного квадранта, что уже гораздо лучше.
              0
              Обычные флюгеры достаточно массивны, чтобы не колыхаться как флаг. Может надо было увеличить инерцию конструкции, добавив «маховик» на ось.
                0
                Можно, конечно. Было бы у меня время и ресурсы, я бы поэкспериментировал, чтобы знать, какова оптимальная масса флюгера. Потому что слишком большая — понятно, что это тоже плохо, увеличится трение и он перестанет реагировать на хоть и постоянный, но слабый ветерок. Значит, нужно будет еще и дорабатывать опоры и в конце концов полностью менять всю конструкцию. И, возможно, не один раз до получения нужного соотношения масса/трение. А пока возможности приобретения такого опыта нет, проще все-таки вносить поправки электронным методом, а не делать это трудноисправляемой механической доработкой.
                ЗЫ: между прочим, есть еще одна теоретическая возможность сделать оптимальный по инерции флюгер. Для этого нужна смазка, обладающая эффектом, противоположным тиксотропии: то есть, имеющая малую вязкость при слабом сдвигающем усилии и резко увеличивающем ее при увеличении усилия. С такой смазкой флюгер перестал бы болтаться при резких порывах и стал бы двигаться более плавно. Но про такую смазку я только слышал, и не уверен, что она работает в нужном мне диапазоне сдвигающих усилий. Так что это тоже требует экспериментов, и как бы даже не поболее чистой механики.
                  0
                  Без этой смазки этот эффект был давно реализован для стрелочных приборов — применяли торможение токами Фуко. У стрелки был металлический флажок, движущийся в разрезе магнита. При резком ускорении возникают большие тормозящие токи. При медленном движении стрелки торможение отсутствует.

                  Т.е. в данной конструкции надо диск (край диска) поместить в разрез подковообразного магнита.

                  P.S. для читателей комментария, не знакомых с токами Фуко www.youtube.com/watch?v=01JCZgpy7tw
                    0
                    Вот видите, вы еще и третий путь мне напомнили (я о нем просто забыл, а ведь когда-то в институте проходил устройство стрелочных приборов. Посыпаю голову пеплом). И все требуют механической доработки без гарантии, что это поможет (магниты тоже бывают очень разные и я лично вот не на раз соображу, как тут и чего посчитать заранее).
                    Так что проще все-таки с программой поэкспериментировать.
                      0
                      В случае магнитов нужна механическая точность изделия. Я бы строил на основе блока головок неисправного жесткого диска.

                      В общем случае действительно программа проще, если сам флюгер не разворачивается в произвольном направлении, генерируя заметное количество ложных данных.
                        0
                        Отсутствие ложных данных обеспечивается при трех условиях: а)флюгер хорошо отбалансирован и самостоятельно не разворачивается в преимущественную сторону при отсутствии ветра; б)его не тормозит существенно трение в опорах и в)не заносит при резких порывах. Первое и второе обеспечить несложно продуманной конструкцией, аккуратным изготовлением и сборкой (у меня для балансировки имеется передвижной груз в передней части), а третьего не происходит при достаточной массе (флюгер из латуни, а груз свинцовый, общим весом около 200-300 грамм), и тем, что порывы ветра не имеют формы строго прямоугольного импульса.
                  0
                  Вдогонку: электронный метод в теории позволяет сымитировать любое поведение флюгера. Не поможет простое осреднение за цикл — пойдем дальше, будем подбирать период осреднения, выдумывать скользящие средние. В принципе нужно просто привлечь теорию: построить динамическую модель поведения флюгера, и наложить на нее желаемый частотный фильтр. Беда тут только в том, что не очень хорошо известно, какой именно фильтр «желаемый». Потому и проще поэкспериментировать, подобрав параметры фильтра вручную.
              0
              Какой результат вернёт алгоритм, если набор данных для усреднения состоит из равного количества значений с азимутами X и X+180° и одинаковыми амлитудами?
                0
                Беру свой легкомысленный ответ (см. ниже) назад. Ошибка действительно теоретически возможна. Написал по этому поводу UPD в статье.
                  0
                  atan2 не поможет — в таких условиях результат будет определяться ошибками округления и может быть любым. А математически тут и вовсе неопределенность.
                    0
                    Я не о неопределенности, с ней все ясно, в изначальном ответе я все сказал. Ошибки округления тут роли не играют, ибо код (16 градаций) грубее любых ошибок. Речь идет о том, что два значения вблизи 90 и 270 градусов, например cos(89°)=0.017 и cos(269°)=-0.017, в среднем дадут ноль, т.е. arccos даст 90° (или 270°, без разницы), а должно быть 179. atan2 эту проблему снимает, специально проверял. Другое дело, что у меня arccos потому и работал без видимых сбоев, что в реальности разброс на 180 градусов (на четыре значения притом) не встречается.
                      0
                      atan2 от двух нолей даст NaN
                        0
                        Двух нулей там быть не может в принципе. Все-так меряются реальные углы.
                          0
                          Все-таки в приведенном вырожденном случае там именно два ноля и будет, если не учитывать ошибки округления. Потому случай и вырожденный.
                0
                Чисто умозрительное построение, конечно, в реальности вероятность такого совпадения практически равна нулю (особенно, если учесть амплитуды, которые у меня не учитываются). Но если все таки случится, то какая разница, какое? 90 или 270 — оба ответа будут правильными. Можно просчитать по алгоритму, какое именно (или смоделировать на контроллере), но мне лень.

                Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                Самое читаемое