Метод расчёта семейства функций Бесселя с помощью быстрого преобразования Фурье при реализации теоретических основ активно-пассивной низкочастотной томографии Текст научной статьи по специальности «Математика»

УДК 621.391.26

МЕТОД РАСЧЁТА СЕМЕЙСТВА ФУНКЦИЙ БЕССЕЛЯ С ПОМОЩЬЮ БЫСТРОГО ПРЕОБРАЗОВАНИЯ ФУРЬЕ ПРИ РЕАЛИЗАЦИИ ТЕОРЕТИЧЕСКИХ ОСНОВ АКТИВНО-ПАССИВНОЙ НИЗКОЧАСТОТНОЙ ТОМОГРАФИИ

© 2005 г. П.А. Стародубцев, А.Н.Иванов

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

В 70-80-е гг. прошлого столетия специалистами Акустического института имени академика Н.Н. Андреева, института прикладной физики РАН было высказано предположение о возможности применения параметрических просветных акустических методов для обнаружения ПО по изменению спектральных, амплитудно-фазовых характеристик комбинационных волн разностной и суммарной частоты как результата параметрического взаимодействия высокочастотных просветных сигналов (ПС) и преобразования их в низкочастотные на гидроакустической барьерной линии.

Но как показала практика, основными недостатками параметрических просветных методов, ограничивающими их широкое применение для обнаружения ПО, являются:

- низкая эффективность параметрического преобразования при нелинейном взаимодействии исходных высокочастотных волн в низкочастотные комбинационные;

- интенсивность волн комбинационных частот зависит от нелинейных характеристик среды (солено -сти, температуры, насыщенности микроорганизмами или газовыми пузырьками) и изменяется в широких пределах от 4 -10-7 до 19 -10-5 относительно интенсивности высокочастотных волн накачки;

- накапливающееся с расстоянием комбинационное взаимодействие сферических волн, рассеянных отдельными рассеивателями на частоте ПС с квазиплоской, инфранизкочастотной (ИНЧ) волной шумо-излучения ПО наблюдается достаточно слабо.

В классических параметрических системах данное взаимодействие составляет 0,01 ... 0,1 %. В системах с применением дополнительных мер и устройств (ультразвуковая пузырьковая область) эффективность преобразования энергии данных волн может составлять 0,1 %.

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

- низкочастотном просветном методе широкомасштабного интенсивного нелинейного взаимодействия и параметрического преобразования ПС возмущенной областью, созданной движущимся ПО во всех направлениях их совместного распространения;

- «квазидифракционном» методе реконструкции процесса взаимодействия низкочастотного ПС с возмущенной областью, созданной движущимся ПО на горизонтально разнесенных приемниках;

- гидроакустической системе как широкомасштабном параметрическом излучателе-приемнике с низкочастотной подсветкой (накачкой) контролируемой трассы, являющейся одновременно и объемной гидроакустической базой.

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

1. Фазовая модуляция с появлением в спектре низкочастотного ПС дополнительных гармоник суммарной и разностной частоты с двумя, четырьмя и более знаками как результат нелинейного взаимодействия ПС с ИНЧ составляющими возмущенной области, созданной движущимся ПО.

2. Нелинейное взаимодействие несущей низкочастотного ПС с наиболее интенсивными ИНЧ составляющими возмущенной области, созданной движущимся ПО, и параметрическое преобразование его в комбинационные волны суммарной и разностной частоты.

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

никновения внутренних корабельных волн. На глубинах меньших 100 м отношение амплитуд внутренних волн к глубине (параметр нелинейности) в этом случае достигает значения 0,1 и более. Возмущения гидрологических полей, в том числе скорости звука, таким волновым движением будут качественно отличаться от аналогичных возмущений, создаваемых акустическими с бесконечно малыми амплитудами, т.е. линейными волнами от шума корпуса ПО. При этом различные части корпуса ПО регенерируют гармонические и квазигармонические полнопериод-ные, полупериодные и четвертьпериодные ИНЧ волны (^) в вертикальной и горизонтальной плоскости при неустойчивом движении ПО. При распространении такой суммарной внутренней волны с конечными амплитудами происходит постепенное изменение формы волны вследствие разности в скоростях движения различных участков его профиля. Периодические волны в процессе распространения эволюционируют в возмущения с «борообразным» (подобным профилю ударных волн) профилем.

Увеличение крутизны внутренней волны соответствует нарастанию высокочастотных гармоник, т.е. обогащению ее спектра высокочастотными составляющими волнового поля, которые вызывают интенсивное изменение спектра ПС. Спектр такого ПС в общем случае содержит бесконечное число гармонических составляющих, частоты которых равны ю 0 ± kQ., а интенсивность этих составляющих пропорциональна значениям функции Бесселя Jk (mk).

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

(L-, 0 г ) = [Po(t)а сп (±Jk (mk )cos ю 0

с 0 k =1

+k Q i (1 - cos 0 i) timk)] ± P±x,

f

t; --

L

(1)

где ti - суммарное время распространения ПС по разным акустическим лучам в пределах г лучевой трубки; асп - коэффициент пространственного затухания на частоте ПС; Li|c0 - время фазового модуляционного (параметрического) взаимодействия низкочастотного ПС с ИНЧ гидрофизическими и другими полями возмущенной области, созданной движущимся ПО, в пределах 1-й лучевой трубки; k - любое вещественное число от 1 до г ; mк - индекс фазовой модуляции

(Y-1 + cos 0i )Po (t)

sin

п — (1 - cos 0 i )

Q, i

2pо(сo)3P(Q i,Li)

п—(1 - cos 0 i )

Qi i

Аю (2)

Jk (mk) - функция Бесселя k индекса от аргумента mk, показывающая зависимость коэффициента модуляции от амплитуды модулирующего сигнала; ю 0 - несущая частота ПС; Аю - девиация частоты ПС; ^ г - ИНЧ колебания возмущенной области ПО с частотами ^ 1, ^ 2,..., ^ i соответственно; Po (t) -исходное давление ПС на излучателе; Li - проекция плоскости сечения области взаимодействия (в том числе и параметрического с появлением комбинационных волн разностной и суммарной частоты) ПС с ИНЧ составляющими возмущенной области на линию, соединяющую излучатель и приемник; P(Q. i, Li) - суммарное давление г-й гармоники ИНЧ в возмущенной области, созданной движущимся ПО; р 0 - плотность водной среды; 0 г - угол между взаимодействующими ПС и гармониками ИНЧ возмущенной области; у - параметр нелинейности среды (морская вода - 3,5 ед., возмущенная область - 10 ... 15 ед.); ± Р±£ - давление волн комбинационных частот суммарной и разностной частот.

Наиболее сложным при проведении расчетов суммарного давления в точке приема является точность определения функция Бесселя. В научных и инженерных расчётах широко используется, в отличие от выражения (1), функция Бесселя Jn(x), формула которой для целого и нецелого значения порядка п имеет следующий вид [1]:

Í „2\k

Jn (Х) ' 2 | k?o k!Г(n + k +1)

(3)

где Г(п) = |е ttn ldt - гамма-функция, которая при

0

целочисленном значении аргумента п может быть интерполирована с помощью выражения Г(п+1)=п!. Данная функция является решением дифференциального уравнения Бесселя для целого и нецелого значения п:

.2 d2y

dx

2 + x— +(x 2 + n 2 2 dx

(x2 + n 2 )y = 0 .

В литературе [1] встречается другое выражение для расчёта функции Бесселя Jn(x):

1 2п

Jn (x) = — J exp [j'x sin 0 - jn0] d0 =

2п

1 2п

= — í exp [ jx sin0]exp [-jn0] d0, 2п 0

которое по существу представляет собой преобразование Фурье от функции ехр/хяшб], где х - постоянная величина, 6 - переменная, эквивалентная времени, п - переменная, эквивалентная частоте.

2п

Выражение (4) с помощью замены ё6^— и

2п

N

N

может быть преобразовано в формулу ДПФ

следующего вида:

1 N-1

Jn (x) = N 2 exp

N к=0

2п, jx sin I — к

1 N

exp

1 N-1

= N 2 exp

N к=0

. , 2п jx sinl — к

1 N

W

2n , - j—пк

N

пк

N

(5)

где W

пк = exp

2п , - j—пк

N

- комплексная величина,

эквивалентная фазовому коэффициенту в формуле ДПФ; к=0, ..., N - 1 - номер выборки по 6.

В ряде задач, где используется функция Бесселя Jn(x) [3], требуется выполнять обработку информации в реальном масштабе времени. Для расчета функции с помощью выражения (5) в данной статье предлагается использовать один из известных алгоритмов БПФ, например, алгоритм БПФ по основанию 2 (алгоритм Кули - Тьюки). В этом случае для одного из N отсчётов переменной x будут рассчитываться сразу N значений функций Jn(x) порядков n=0,..., N - 1.

Результаты подобных расчетов в системе инженерных и научных расчётов MATLAB 6.5 отображены на рис. 1. В вычислительном эксперименте выполнялись расчёты функций Бесселя для различных порядков n = 0, ...,15 (с шагом 1) для области определения x = 0, ... ,16 (с шагом 0,01). В самом начале выполнялись расчёты функций Бесселя по традиционным соотношениям (1), результат вычислений с помощью функции MATLAB besselj показан на рис. 1 а. Затем выполнялся расчёт с помощью алгоритма БПФ (функция MATLAB fft), рис. 1 б. На указанных рисунках показаны графики функции Бесселя только для порядка n = 0, ..., 7. После этого результаты вычислений сравнивались между собой. На рис. 1 в и г показаны графики абсолютной ошибки вычислений для различных диапазонов x = 0, ..., 16 и x = 0, ..., 4, которая рассчитывалась как разница AJn(x) между результатами двух различных методов.

Как видно из графиков на рис. 1 в и г, разница растёт быстрее с увеличением порядка n, а при n = 7 наблюдается самый быстрый её рост. Но даже для больших порядков на рис. 1 г ошибка вычислений не превышает 10 - 3 для x = 0, ..., 4. Такую особенность данного метода вычисления можно объяснить двумя причинами: периодичностью спектра при вычислении ДПФ и явлением наложения (aliasing), наблюдаемым при неправильной дискретизации непрерывных величин [3].

Jn(x) 0,5 0 -0,5

Jn(x) 0,5 0 -0,5

Jx)

0,2 0 -0,2 -0,4

AJn(x)-10-0,8 0,4 0 -0,4 -0,8

0 2 4 6 8 10 а

0

12 14 x

0 2 4 6 8 10 12 14 x

б

0 2 4 6 8 10 12 14 x

0 0,5

1

2,5 3

3,5

1,5 2

г

Рис. 1. Результаты расчётов в системе МЛТЬЛЕ: а - по выражению (1); б - расчёт с помощью БПФ; в и г ошибка вычисления для х = 0, ..., N и х = 0, ..., N/4

Известно, что результат дискретного преобразования Фурье имеет периодическую и симметричную структуру, поэтому из него надо выбрать только часть составляющих для п = 0, ..., N/2 - 1. В выражении (3)

. , 2п

гармоническая функция ехр

jx sinl —к

1 N

зависит от

двух переменных: x, которая играет роль частоты, и sinl к I, которое представляет собой периодически

N

изменяющееся в диапазоне от -1 до 1 дискретное время. Поэтому при росте переменной х выше N/2 происходит наложение периодических составляющих преобразования Фурье. Для иллюстрации приведём трёхмерные графики функции Бесселя (рис. 2), рассчитанные в МЛ ТЬЛЕ рассмотренными выше методами для значений порядка п = 0, ..., 63 (с шагом 1) и независимой величины х = 0, ..., 31 (с шагом 0,1).

0

2

в

7

1

. 0,5

0 -

-0,5 80

я

40

1

0,5 0 -0,5

40

40

20

и 1 0

1

0,5

Иг

-0,5 80

40

Из рис. 2 в видно, что различие между результатами расчётов равно нулю, если порядок n = 0, ..., 31, а независимая переменная х = 0, ..., 31. На основании вышесказанного следует, что при вычислении функций Бесселя для порядков n = 0, ..., N/2 следует рассчитывать с помощью алгоритма БПФ сразу N значений, но учитывать как результат только первые N/2 значения, потому что остальные составляющие будут периодически повторяться. Кроме этого, следует помнить, что эти вычисления более точны методом БПФ только при x<N/2 (в данном случае играет роль частоты Найквиста). При x>N/2 абсолютное значение ошибки будет резко увеличиваться из-за эффекта наложения, наблюдаемого в ДПФ. Но при проведении расчетов для низкочастотной акустической томографии данных точностных характеристик вполне достаточно.

Для сравнительной оценки быстродействия предлагаемого метода расчёта в системе MATLAB 6.5 проводился вычислительный эксперимент. С помощью функций программного таймера tic и toc, описанных в [2], определялось время, затраченное на вычисления. Полученные при этом результаты для традиционного [1] и предлагаемого методов расчёта приведены в таблице и на рис. 3. Они показывают, что предлагаемый метод имеет в 3.8 раз меньшее быстродействие, чем расчёт функции Бесселя традиционным способом по выражению (3).

160 120 80 40 0

л

m

6 7 8 9 2r

Рис. 2. К описанию особенностей вычисления функций Бесселя: а - традиционным методом, б - с помощью БПФ, в - разница между ними

Рис. 3. К сравнительному анализу методов расчёта: сплошная линия - расчёт традиционным способом, пунктирная линия - расчёт с помощью БПФ

Время расчёта семейства функций Бесселя для n = 0, ..., 2Г - 1

Таблица

Метод расчёта Быстродействие r

2 3 4 5 6 7 8 9 10

Время расчёта, с

Традиционный 0,12 0,151 0,29 0,521 0,861 1,462 3,675 8,122 25,397

С помощью БПФ 0,37 0,821 1,783 5,047 10,766 21,761 43,342 86,405 172,39

Примечание: следует учесть, что предлагаемый метод расчёта (с помощью БПФ) требует в два раза больше значений по п = 0, ..., 2Г - 1, т.е. г = 3, ..., 11 вместо г = 2, ..., 10.

0

и

х

а

х

б

3

4

5

и

х

в

Но так как при реализации активно-пассивной томографии первым анализируется и уточняется энергетический спектр совокупности ПС, зафиксированного приемником, а в случае, когда в энергетическом спектре отсутствуют достаточно интенсивные спектральные компоненты исходного низкочастотного сигнала, однозначно осуществляется операция компрессии через аддитивный компенсатор помех, которая позволяет путем набора наиболее оптимальных весовых окон выявить и усилить спектральные компоненты, отсутствующие в первичном Фурье-спектре, то точ-

ность вычисления функции Бесселя оказывается более приоритетной, чем скорость ее расчета.

Литература

1. Математика // Большой энциклопедический словарь/ Гл. ред. Ю.В. Прохоров: 3-е изд. М., 1998.

2. Потёмкин В.Г. Система инженерных и научных расчётов МЛТЬЛБ 5.x: В 2 т. М., 1999.

3. Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов / Пер. с англ. А.Л. Зайцева, Э.Г. На-заренко, Н.Н. Тетёкина. М., 1978..

Морской государственный университет им. адмирала Г.И. Невельского, г. Владивосток 1 июля 2005 г.

УДК 007:631.2

О МОДЕЛИРОВАНИИ НЕЧЕТКИХ ЭКСПЕРТНЫХ ЗНАНИЙ ПО ТЕХНОЛОГИЧЕСКОЙ РЕГУЛИРОВКЕ КОМБАЙНА

© 2005 г. В.А. Алуханян, Л.В. Борисова, В.П. Димитров

Введение

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

Для сложных технических систем, каким является зерноуборочный комбайн, построить адекватное математическое описание процесса принятия решений при функционировании машины и техническом обслуживании очень сложно. Специфика выбора модели принятия решений, предназначенной для решения задач технического обслуживания уборочной техники, заключается в следующем. Во-первых, достаточно сложные условия эксплуатации зерноуборочной техники в меньшей мере поддаются точному количественному описанию, а следовательно, большая часть информации о стратегиях принятия решений (ПР) (которая представлена в словесной форме) исходит непосредственно от эксперта. Во-вторых, алгоритмы ПР предназначены для оперативного использования в условиях уборки, т. е. должны работать в реальном времени, поэтому применение точных методов оптимизации, как правило, исключается вследствие их трудоемкости. В-третьих, то обстоятельство, что алгоритмы должны работать в качестве «советчика» механизатора, предъявляет к ним требование учитывать качественную информацию, представленную в лингвистической форме [1].

Анализ предметной области [2 - 5] позволил выделить два типа задач принятия решений. К первому типу относятся задачи, в результате решения которых происходит выбор значений параметров рабочих органов зерноуборочного комбайна (частоты вращения молотильного барабана и т.д.). К этому типу примыкают задачи, в которых по имеющимся значениям регулируемых параметров и факторов внешней среды необходимо определить значения показателей качества технологического процесса (например, потери зерна и др.).

Ко второму типу относятся задачи определения причин, вызвавших отклонения значений показателей качества от допустимых. Другими словами по известному внешнему признаку отклонения необходимо найти причину, которая вызвала это отклонение (корректировка технологических регулировок). К данному типу относятся задачи поиска неисправностей в агрегатах и системах комбайна.

Идентификация предметной области

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

Обозначим через X, У, множество значений факторов внешней среды, существенно влияющих на выбор выходного параметра V. В нашем случае в качестве выходных параметров выступают регулировочные параметры рабочих органов зерноуборочного комбайна. В общем процесс принятия решений имеет несколько выходных параметров (конечное множест-