Как было сказано выше, для получения полной информации о трех компонентах вектора скорости ветра в декартовой системе координат с помощью КДЛ необходимо провести измерения, как минимум, при трех различных положениях зондирующего пучка в пространстве. В соответствии с рисунком 1, поясняющим принцип конического сканирования атмосферы, на заданной высоте h, определяемой необходимостью измерения, радиальную составляющую скорости ветра можно записать следующим образом:

,

где и - угол азимута лазерного луча, который изменяется в интервале [0; 2р), а ц - угол места лазерного луча, изменяющийся в диапазоне [0; р/2].

Таким образом, при измерении радиальной скорости ветра в одной точке записывается только одно уравнение при трех неизвестных. При измерении радиальной составляющей скорости ветра в трех независимых точках в пространстве на заданной высоте для трех неизвестных параметров записывается система из трех независимых уравнений. Решение данной системы уравнений позволяет определить полный вектор скорости ветра. При непрерывном вращении лазерного луча по углу азимута путем отбора его дискретных положений получается переопределенная система уравнений, которая решается методом наименьших квадратов.

Отметим, что измерение полного вектора скорости ветра при коническом сканировании воздушного пространства лазерным лучом на различных высотах достигается путем соответствующего изменения фокусного расстояния оптической системы КДЛ. Ключевым моментом принципа конического сканирования является предположение о “замороженности" ветровых характеристик атмосферы за время сканирования на заданной высоте. Поэтому в разработанном алгоритме измерение на одной высоте происходит в течение одной секунды.

Рисунок 1. Принцип конического сканирования атмосферы лазерным лучом.

В случае, когда ветровые характеристики атмосферы за период полного сканирования лазерным лучом по углу азимута на заданной высоте остаются постоянными, радиальная составляющая вектора скорости ветра изменяется по синусоидальному закону, как это показано на рисунке 2.

Рисунок 2. Закон изменения радиальной составляющей скорости ветра при постоянстве ветровых характеристик атмосферы при коническом сканировании.

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

,

где a, b - аппроксимирующие коэффициенты.

Направление ветра по азимуту определяется как угол , соответствующий максимальной положительной скорости. Величина горизонтальной составляющей скорости определяется амплитудой гармонической функции. Вертикальным смещением полученного графика определяется величина вертикальной составляющей вектора скорости.

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

Основными задачами программно-алгоритмического обеспечения КДЛ являются:

1. Определение спектральных характеристик гетеродинного сигнала на выходе приемного устройства.

2. Обнаружение полезного сигнала в спектральных данных и формирование набора частот, характеризующего скоростные параметры воздушных масс.

3. Подготовка спектральных данных для аппроксимации путем удаления выбросов из набора частот с целью уменьшения ошибочных вычислений параметров ветра.

4. Проведение операции аппроксимации с целью определения искомых параметров функции приближения.

5. Окончательный расчёт параметров ветра.

Блок-схема алгоритма вычисления параметров ветра приведена на рисунке 3.

Рисунок 3. - Блок-схема алгоритма вычисления параметров ветра.

Алгоритм вычисления параметров ветра рассмотрим на примере обработки сигнала с профилометра лидарного ветрового ПЛВ-300 [3], разработанного ООО “Лазерные системы" в Санкт-Петербурге.

Сигнал с выхода приемника обратно рассеянного излучения содержит переменную составляющую с частотой до 50 МГц. Обрабатываемая полоса частот непосредственно задается диапазоном скоростей, определяемых лидаром. Аналого-цифровое преобразование осуществляется на частоте 100 МГц. Объем одиночной выборки содержит 512 отсчетов временного сигнала. Длительность выборки по времени - 5.12 мкс. Полученный временной сигнал поступает на блок быстрого преобразования Фурье (БПФ). На выходе блока БПФ вычисляется квадрат модуля комплексного Фурье-спектра.

Квадраты модулей спектров усредняются. Количество спектров в серии - 4096. Общая продолжительность одной серии съема данных - ~21 мс. На рисунке 4 представлена схема обработки сигнала.

Рисунок 4. - Схема обработки сигнала в единичном измерении.

На рисунке 5 показан вид результирующего спектрального сигнала.

Рисунок 5. - Результирующий спектральный сигнал.

Спектр содержит искомую дискретную компоненту в виде колоколообразного импульса. Отношение сигнал/шум находится в диапазоне 3 - 10 для нормальных атмосферных условий с МДВ в диапазоне 7 - 20 км. Ширина сигнала, определяемая по уровню 0.7 от максимального значения, составляет 1 - 5 МГц. Задача алгоритма заключается в поиске дискретной составляющей полезного спектра, определении ее максимального значения и вычислении соответствующей ему частоты доплеровского сдвига:

,

где N₀ - положение, соответствующее нулевой скорости ветра.

После этого вычисляется абсолютное значение радиальной составляющей скорости ветра

,

где м - рабочая длина волны КДЛ.

Для выделения полезного сигнала несущего информацию о доплеровском сдвиге частоты, необходимо провести обработку спектра, усредненного за 4096 реализаций. Основные этапы обработки спектра представлены ниже:

1. Нормировка спектра S.

Данная операция позволяет выделить изолинию спектра. Характерный вид спектра с неравномерной фоновой характеристикой показан на рисунке 6. Неравномерность полученного спектра обусловлена шумовыми характеристиками лазера, приемника и трака АЦП.

a) С помощью медианного фильтра с фиксированным окном спектр сглаживается, из него удаляются резко выделяющиеся значения (в том числе полезный колоколообразный сигнал). На выходе фильтра получается спектр . На рисунке 7 показан сглаженный спектр, исходный вид которого приведен на рисунке 6.

b) Производится нормировка: .

Рисунок 6. - Усреднённый спектральный сигнал.

Рисунок 7. - Изолиния спектрального сигнала.

2. Определение частоты доплеровского сдвига в спектре полезного сигнала.

a) Выделение полезного колоколообразного сигнала на нормированном спектре. Для этого производится сглаживание “скользящим" или “бегущим” средним с маленьким окном, чтобы удалить случайные выбросы. В результате на выходе получается спектр , показанный на рисунке 8.

b) Расчет математического ожидания m и среднеквадратического отклонения s спектра .

Рисунок 8. - Спектральный сигнал после нормировки и сглаживания.

Для точек нормированного фильтрованного спектра, удовлетворяющих условию спектра >, рассчитывается центральная частота колоколообразного импульса , множитель k - задает порог чувствительности и лежит в пределах от 3 до 4, s - стандартное отклонение, m - математическое ожидание. Положение центральной частоты рассчитывается в соответствии с известным центроидным методом по формуле центра масс:

,

где - частота i-й точки спектра,

- масса i-й точки (значение по амплитуде).

Кроме приведенного выше центроидного метода могут быть использованы такие методы, как метод гауссовой аппроксимации [4] либо метод с использованием порядковых статистик [5].

На рисунке 9 показан уровень максимума спектрального сигнала и рассчитанный уровень порога .

3. Пункты 1 и 2 повторяются для окна медианного фильтра шириной в 2 раза большей. Форма сигнала (ширина колоколообразного импульса) зависит от текущих метеоусловий и может меняться. Алгоритм повторяется с большим окном медианного фильтра, чтобы не пропустить импульсы полезного сигнала различной ширины. Это позволяет использовать алгоритм независимо от текущей формы сигнала.

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

Рисунок 9. - Результат сканирования на фиксированной высоте.

На рисунке 9 представлен результат сканирования на фиксированной высоте 90 метров в течении нескольких минут. На графике по горизонтальной оси отложены номера точек на кругу сканирования, по вертикальной оси - значения радиальных скоростей, соответствующих центральной частоте найденных на каждом кругу сканирования полезных спектральных сигналов.

4. Алгоритм удаления выбросов из набора частот

В рассматриваемом КДЛ линия зондирования за время ~1 секунду описывает коническую поверхность с углом раствора 22є по отношению к вертикальной оси. Для измерения скорости ветра на заданной высоте осуществляется за один цикл три оборота сканирования. Четвертый оборот используется для определения направления ветра.

Наборы частот с каждого оборота объединяются в общий массив по принципу:

,

где n - число радиальных проекций за один оборот, в качестве столбцов выступает номер оборота.

Спектры радиальных проекций помимо сигнальной составляющей также содержат шум, который в зависимости от атмосферных условий, может в меньшей или большей степени мешать выделению полезного сигнала. Таким образом, в результирующий набор частот попадают нежелательные шумовые частоты, которые должны быть удалены до вычисления параметров аппроксимирующей функции. Это повышает точность аппроксимации, и, соответственно, точность определения параметров скорости ветра.

На рисунке 10 приведена блок-схема алгоритма удаления выбросов из набора частот. Для удаления выбросов из набора частот рассчитываются среднее значение (m) и среднеквадратичное отклонение (s) набора частот . Далее:

1. По порогу отсекаются все точки, амплитуда которых превышает заданный порог , где .

2. Оставшиеся точки (подмножество ) передаются в алгоритм аппроксимации.

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

Рисунок 10. - Блок-схема алгоритма удаления выбросов из набора частот.

Рисунок 11. - Аппроксимация полученных данных гармонической зависимостью.

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

5. Алгоритм определения параметров аппроксимирующей функции

Аппроксимирующая функция задается выражением:

.

Без потери общности ее общий вид показан на рисунке 12.

Рисунок 12. - Вид аппроксимирующей функции.

Очевидно, выражение для вектора скорости имеет вид:

.

Горизонтальная составляющая вектора скорости:

.

Вертикальная составляющая вектора скорости:

.

Как указывалось выше, на практике имеем:

.

Функция зависит от переменной и трёх параметров a, b и .

Тогда для n измерений за один оборот сканирования можно записать следующее выражение:

где , а .

Очевидно, для определения параметров a, b, с в соответствии с методом наименьших квадратов необходимо решить систему уравнений:

А именно:

Решение системы уравнений содержит искомые параметры a, b и c, по которым вычисляется горизонтальная и вертикальная составляющие вектора скорости ветра, а также его направление. Отметим, что методы, алгоритмы и сама процедура определения параметров a, b и c в рассмотренной выше постановке детально проработаны в диссертации Смалихо И.Н. [6].

Описанная серия измерений производится при фиксированной фокусировке приемо-передающего телескопа, соответствующей высоте, на которой производятся измерения. Для измерения вертикального профиля скорости ветра эта серия повторяется при нескольких различных фокусировках. Основным конечным результатом измерений является вертикальный профиль скорости ветра, измеренный в нескольких заранее заданных точках по высоте, привязанный к топографическим координатам. Кроме того, должно быть известно, как сориентированы оси x, y относительно направления на север (с точностью 0.2° - 0.3°). Отображение вертикального профиля скорости ветра в зависимости от высоты на мониторе оператора ПЛВ-300 для примера показано на рисунке 13.

Рисунок 13. - Отображение на мониторе оператора данных о профиле ветра.

6. Выводы