Акустический журнал, 2022, T. 68, № 6, стр. 689-696

Локализация источника в переменном по трассе волноводе в условиях неполной информации о пространственной изменчивости среды распространения

А. Г. Сазонтов ab*

a Институт прикладной физики РАН
Нижний Новгород, Россия

b Нижегородский государственный университет им. Н.И. Лобачевского
Нижний Новгород, Россия

* E-mail: sazontov@ipfran.ru

Поступила в редакцию 20.06.2022
После доработки 25.07.2022
Принята к публикации 27.07.2022

Полный текст (PDF)

Аннотация

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

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

ВВЕДЕНИЕ

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

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

В настоящей работе построен адаптивный алгоритм пониженного ранга RARE (“rank reduction”), позволяющий оценить искомые координаты источников без знания истинных параметров неоднородного акустического волновода. Проводимое рассмотрение рассчитано на наихудший сценарий приема и предполагает ограниченность по норме вектора рассогласования между принятым звуковым полем и его расчетной моделью, при этом процедура адаптации заключается в нахождении робастного вектора отклика АР, обеспечивающего максимум выходной мощности процессора и удовлетворяющего наложенным на него ограничениям (см. [3 ] и цитируемую там литературу). Представлены результаты сравнительного анализа эффективности предложенного способа оценивания с обычным методом MUSIC, не учитывающим эффекты взаимодействия мод.

ПОСТАНОВКА ЗАДАЧИ. ИСХОДНЫЕ СООТНОШЕНИЯ

Пусть в точках с координатами ${{{\mathbf{\theta }}}_{1}} = {{({{r}_{1}},{{z}_{1}})}^{T}}, \ldots ,$ ${{{\mathbf{\theta }}}_{J}} = {{({{r}_{J}},{{z}_{J}})}^{T}}$ волноводного канала находятся $J$ источников звука, излучающих частично–когерентные узкополосные сигналы с одинаковой несущей частотой ${{\omega }_{0}}.$ (Верхний индекс $T$ обозначает операцию транспонирования.) Прием осуществляется линейной вертикальной АР, состоящей из $N$ элементов, расположенных на горизонтах $\{ {{z}_{n}}\} _{{n\, = \,1}}^{N}$. (Начало координат по дальности выбрано в месте установки АР.)

В узкополосном приближении результирующее поле на входе АР характеризуется $N$–мерным вектором наблюдения ${{{\mathbf{x}}}_{l}}$:

(1)
${{{\mathbf{x}}}_{l}} = \sum\limits_{j{\kern 1pt} = {\kern 1pt} {\kern 1pt} 1}^J {{\kern 1pt} {\mathbf{g}}{\kern 1pt} {\kern 1pt} ({{{\mathbf{\theta }}}_{j}})} {{s}_{j}}(l) + {{{\mathbf{n}}}_{l}},\,\,\,\,\,l = 1,2, \ldots ,L.$

Здесь $l$ – номер выборочного отсчета, ${\mathbf{g}}({{{\mathbf{\theta }}}_{j}}) = [\left. {G(0,{{z}_{1}}} \right|{{{\mathbf{\theta }}}_{j}}), \ldots ,G(\left. {0,{{z}_{N}}} \right|{{{\mathbf{\theta }}}_{j}}){{]}^{T}}$ – вектор отклика АР при приеме сигнала от $j$-го источника с комплексной огибающей ${{s}_{j}}(l),$ $\{ G(\left. {0,{{z}_{n}}} \right|{{{\mathbf{\theta }}}_{j}})\} _{{n = 1}}^{N}$ – функции Грина, связывающие координаты соответствующего источника с координатами приемных элементов, ${{{\mathbf{n}}}_{l}}$ – вектор аддитивного белого шума, а $L$ – объем входной выборки. Задача состоит в построении адаптивного алгоритма обработки, позволяющего по принятой выборке $\{ {{{\mathbf{x}}}_{l}}\} _{{l = 1}}^{L}$ оценить искомые координаты источников без знания пространственной изменчивости акустического волновода.

При дальнейшем анализе будем считать, что число источников известно, а случайные векторы ${{{\mathbf{s}}}_{l}}$ и ${{{\mathbf{n}}}_{l}}$ взаимно не коррелированы и характеризуются следующими ковариационными матрицами

$\left\langle {{{{\mathbf{s}}}_{l}}{\mathbf{s}}_{q}^{ + }} \right\rangle = {\mathbf{S}}{\kern 1pt} {{\delta }_{{lq}}},\,\,\,\,\;\left\langle {{{{\mathbf{n}}}_{l}}{\mathbf{n}}_{q}^{ + }} \right\rangle = \sigma _{n}^{2}{\mathbf{I}}{{\delta }_{{lq}}},$

где ${\mathbf{S}} \in {{C}^{{J \times J}}}$ – матрица, описывающая взаимные корреляции излученных сигналов, $\sigma _{n}^{2}$ – неизвестный уровень шума, ${\mathbf{I}}$ – единичная матрица размерности $N \times N,$ ${{\delta }_{{lq}}}$ – символ Кронекера, а ${{( \cdot )}^{ + }}$ и $\left\langle \cdot \right\rangle $ означают операции эрмитового сопряжения и статистического усреднения, соответственно.

Ниже при нахождении ожидаемой реплики сигнала мы воспользуемся волновым подходом, в рамках которого функция Грина $G(0,\left. {{{z}_{n}}} \right|{{{\mathbf{\theta }}}_{j}})$ может быть представлена в виде суперпозиции конечного числа $M$ распространяющихся локальных нормальных мод. В адиабатическом приближении для $G(0,\left. {{{z}_{n}}} \right|{{{\mathbf{\theta }}}_{j}})$ имеем [4 ] :

(2)
$\begin{gathered} G(0,\left. {{{z}_{n}}} \right|{{{\mathbf{\theta }}}_{j}}) = \sum\limits_{m\, = 1}^M {{{a}_{m}}({\mathbf{\theta }}{}_{j}){{\varphi }_{m}}(0,{{z}_{n}})} , \\ {{a}_{m}}({\mathbf{\theta }}{}_{j}) = \frac{{{{\varphi }_{m}}({{r}_{j}},{{z}_{j}})}}{{\sqrt {8{{\pi }}{\kern 1pt} {{k}_{m}}(0){{r}_{j}}} }}\exp \left[ {i\int\limits_0^{{{r}_{j}}} {{{k}_{m}}} (r)dr + {{i{{\pi }}} \mathord{\left/ {\vphantom {{i{{\pi }}} 4}} \right. \kern-0em} 4}} \right]. \\ \end{gathered} $

Здесь ${{\varphi }_{m}}(0,{{z}_{n}})$ и ${{\varphi }_{m}}({{r}_{j}},{{z}_{j}})$ – собственные функции $m$-ой моды в месте расположения $n$-го приемного элемента и $j$-го источника излучения, а ${{k}_{m}}(r)$ – соответствующее горизонтальное волновое число в сечении $r = const.$ В указанном приближении вектор отклика АР записывается в виде: ${\mathbf{g}}({{{\mathbf{\theta }}}_{j}}) = {\mathbf{U}}\,{\mathbf{a}}({{{\mathbf{\theta }}}_{j}}),$ где ${\mathbf{U}}$ – матрица модовой структуры размерности $N \times M$ с элементами ${{U}_{{nm}}} = {{\varphi }_{m}}(0,{{z}_{n}}),$ а ${\mathbf{a}}({\mathbf{\theta }}) \in {{С}^{{M \times 1}}}$ – модовый вектор с компонентами $a{}_{m}({\mathbf{\theta }}),$ определяемыми формулой (2).

Взаимодействие нормальных волн в переменном по трассе канале приводит к связи между их амплитудами, которую можно учесть путем введения матрицы связи ${{{\mathbf{C}}}_{j}}$ размерности $M \times M.$ В результате для ${\mathbf{g}}({\mathbf{\theta }}{}_{j})$ следует

(3)
${\mathbf{g}}({{{\mathbf{\theta }}}_{j}}) = {\mathbf{U}}{{{\mathbf{C}}}_{j}}{\mathbf{a}}({{{\mathbf{\theta }}}_{j}}).$

В общем случае матрица ${{{\mathbf{C}}}_{j}}$, фигурирующая в определении (3), содержит ${{M}^{2}}$ элементов и не имеет определенной структуры. Однако, учитывая, что в неоднородном по трассе волноводе взаимодействие в основном происходит между модами близких номеров (а следовательно, элементы матрицы связи, удаленные от главной диагонали, вносят пренебрежимо малый вклад в ее формирование), ${{{\mathbf{C}}}_{j}}$ можно аппроксимировать ленточной матрицей, содержащей лишь $P \ll {{M}^{2}}$ ненулевых элементов. Последние удобно рассматривать как компоненты вектора

(4)
${{{\mathbf{c}}}_{j}} = {{[{{c}_{1}}({{{\mathbf{\theta }}}_{j}}),{{c}_{2}}({{{\mathbf{\theta }}}_{j}}), \cdots ,{{c}_{P}}({{{\mathbf{\theta }}}_{j}})]}^{T}} \in {{С}^{{P \times 1}}}.$

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

В рамках рассматриваемой аппроксимации матрицы связи справедливо непосредственно проверяемое соотношение [5–7 ] :

(5)
${\kern 1pt} {{{\mathbf{C}}}_{j}}{\mathbf{a}}({{{\mathbf{\theta }}}_{j}}) = {\mathbf{T}}({{{\mathbf{\theta }}}_{j}}){{{\mathbf{c}}}_{j}},$
где ${\kern 1pt} {\mathbf{T}}({\mathbf{\theta }}) = [{{{\mathbf{E}}}_{1}}{\mathbf{a}}({\mathbf{\theta }})\,{{{\mathbf{E}}}_{2}}{\mathbf{a}}({\mathbf{\theta }}) \cdots {{{\mathbf{E}}}_{P}}{\mathbf{a}}({\mathbf{\theta }})] \in {{C}^{{M \times P}}},$ а

$\begin{gathered} {{[{{{\mathbf{E}}}_{p}}]}_{{\,nm}}} = \left\{ \begin{gathered} 1,\,\,\,\,если\,\,[{{{\mathbf{C}}}_{j}}]{{{\kern 1pt} }_{{nm}}}\,\, = {{c}_{p}}({{{\mathbf{\theta }}}_{j}}); \hfill \\ 0,\,\,\,\,в\,\,противном\,\,\;случае, \hfill \\ \end{gathered} \right. \\ p = 1, \ldots ,P. \\ \end{gathered} $

Простое доказательство формулы (5) основано на возможности представления матрицы ${{{\mathbf{C}}}_{j}}$ в виде

${\kern 1pt} {{{\mathbf{C}}}_{j}} = \sum\limits_{p = 1}^P {{{{\mathbf{E}}}_{p}}{{c}_{p}}({{{\mathbf{\theta }}}_{j}})} .$

В результате

$\begin{gathered} {\kern 1pt} {{{\mathbf{C}}}_{j}}{\kern 1pt} {\kern 1pt} {\mathbf{a}}({{{\mathbf{\theta }}}_{j}}) = \sum\limits_{p = 1}^P {{{{\mathbf{E}}}_{p}}{\mathbf{a}}({{{\mathbf{\theta }}}_{j}}){{c}_{p}}({{{\mathbf{\theta }}}_{j}})} = \\ = [{{{\mathbf{E}}}_{1}}{\mathbf{a}}({{{\mathbf{\theta }}}_{j}}) \cdots {{{\mathbf{E}}}_{P}}{\mathbf{a}}({{{\mathbf{\theta }}}_{j}})]\left[ {\begin{array}{*{20}{c}} {{{c}_{1}}({{{\mathbf{\theta }}}_{j}})} \\ \vdots \\ {{{c}_{P}}({{{\mathbf{\theta }}}_{j}})} \end{array}} \right] \equiv {\mathbf{T}}({{{\mathbf{\theta }}}_{j}}){{{\mathbf{c}}}_{j}}. \\ \end{gathered} $

С учетом (4) и (5) исходный вектор наблюдения (1) может быть переписан следующим образом

(6)
${{{\mathbf{x}}}_{l}} = {\mathbf{G}}({\mathbf{\theta }}){{{\mathbf{s}}}_{l}} + {{{\mathbf{n}}}_{l}},\,\,\,\,\,\,l = 1,2, \ldots ,L.$

Здесь ${\mathbf{G}}({\mathbf{\theta }}) = {\mathbf{U}}[{{{\mathbf{C}}}_{1}}{\mathbf{a}}({{{\mathbf{\theta }}}_{1}}), \ldots ,{{{\mathbf{C}}}_{J}}{\mathbf{a}}({{{\mathbf{\theta }}}_{J}})] \equiv $ ${\mathbf{U}}[{\mathbf{T}}({{{\mathbf{\theta }}}_{1}}){\kern 1pt} {{{\mathbf{c}}}_{1}},$ $ \ldots ,{\mathbf{T}}({{{\mathbf{\theta }}}_{J}}){\kern 1pt} {{{\mathbf{c}}}_{J}}] \in {{С}^{{N \times J}}}$ – передаточная матрица канала, зависящая от искомого вектора ${\mathbf{\theta }} = {{({\mathbf{\theta }}_{1}^{T}, \ldots ,{\mathbf{\theta }}_{J}^{T})}^{T}}$ размерности $2J{\kern 1pt} \times 1,$ определяющего пространственные положения источников, ${{{\mathbf{s}}}_{l}} = {{[{{s}_{1}}(l), \ldots ,{{s}_{J}}(l)]}^{T}} \in {{C}^{{J \times 1}}}$ – вектор комплексных огибающих излученных сигналов.

Для статистически независимых сигналов и помех ковариационная матрица вектора (6) дается соотношением:

(7)
${{{\mathbf{\Gamma }}}_{{\mathbf{x}}}} = {\mathbf{G}}({\mathbf{\theta }}){\mathbf{S}}{{{\mathbf{G}}}^{ + }}({\mathbf{\theta }}) + \sigma _{n}^{2}{\mathbf{I}}{\kern 1pt} .$

В реальных ситуациях истинная матрица ${{{\mathbf{\Gamma }}}_{{\mathbf{x}}}}$ неизвестна, поэтому на практике вместо нее используется оценочная матрица ${{{\mathbf{\hat {\Gamma }}}}_{{\mathbf{x}}}}$ размерности $N \times N,$ равная

${{{\mathbf{\hat {\Gamma }}}}_{{\mathbf{x}}}} = \frac{1}{L}\sum\limits_{l = 1}^L {{{{\mathbf{x}}}_{l}}{\mathbf{x}}_{l}^{ + }} .$

Одним из наиболее распространенных способов локализации является алгоритм MUSIC [8 ] . Эта процедура основана на использовании информации, содержащейся в системе собственных векторов $\{ {\mathbf{\hat {\psi }}}\} _{{i = 1}}^{N}$ выборочной матрицы ${{{\mathbf{\hat {\Gamma }}}}_{{\mathbf{x}}}}{\kern 1pt} :$

${{{\mathbf{\hat {\Gamma }}}}_{{\mathbf{x}}}}{{{\mathbf{\hat {\psi }}}}_{i}} = {{\lambda }_{i}}{{{\mathbf{\hat {\psi }}}}_{i}},\,\,\,\,i = 1, \ldots N,$
где ${{\lambda }_{1}}, \ldots ,{{\lambda }_{N}}$ – положительные собственные числа, пронумерованные в порядке убывания, т.е. ${{\lambda }_{1}} \geqslant {{\lambda }_{2}} \geqslant \cdots \geqslant {{\lambda }_{N}}$. Первые $J$ старших собственных векторов формируют сигнальное подпространство ${{{\mathbf{\hat {\Psi }}}}_{{\mathbf{s}}}} = [{{{\mathbf{\hat {\psi }}}}_{1}} \ldots {{{\mathbf{\hat {\psi }}}}_{J}}] \in {{C}^{{N \times J}}}$, а $N - J$ оставшихся векторов – шумовое ${{{\mathbf{\hat {\Psi }}}}_{{\mathbf{n}}}} = [{{{\mathbf{\hat {\psi }}}}_{{J + 1}}}, \ldots ,{{{\mathbf{\hat {\psi }}}}_{N}}] \in {{C}^{{N \times (N - J)}}}.$

Для рассматриваемого способа обработки искомые координаты источника и неизвестные элементы матрицы связи могут быть найдены из условия ортогональности ожидаемого вектора отклика (3) к векторам шумового подпространства

(8)
$\{ {{{\mathbf{\hat {\theta }}}}_{j}},{{{\mathbf{\hat {c}}}}_{j}}\} _{{j = 1}}^{J} = \arg \mathop {\min }\limits_{{\mathbf{\theta }},\,{\mathbf{c}}}^ {{\left\| {{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + }{\mathbf{UC}}({\mathbf{c}}){\mathbf{a}}({\mathbf{\theta }})} \right\|}^{2}},$
и процедура оценивания сводится к поиску $J$ минимумов целевой функции в (8) в трехмерной области параметров. Последнее требует больших вычислительных затрат.

Принимая во внимание представление (5), переформулируем исходный критерий (8) следующим образом

(9)
где ${\mathbf{Q}}{\kern 1pt} ({\mathbf{\theta }}) = {{{\mathbf{V}}}^{ + }}({\mathbf{\theta }}){{{\mathbf{\Psi }}}_{{\mathbf{n}}}}{{{\mathbf{\hat {\Psi }}}}_{{\mathbf{n}}}}{\mathbf{V}}({\mathbf{\theta }}) \in {{С}^{{P \times P}}},$ а ${\mathbf{V}}({\mathbf{\theta }}) = $ $ = {\mathbf{UT}}({\mathbf{\theta }}) \in {{C}^{{N \times {\kern 1pt} P}}}.$ Отметим, что поскольку $rank\{ {\mathbf{Q}}\} = \min {\kern 1pt} {\kern 1pt} (N - J,M,P),$ то при выполнении условия $N - J > M,$ $Q{\kern 1pt} ({\mathbf{\theta }})$ будет матрицей полного ранга, если $P < M,$ и следовательно, в рамках данного способа локализации число оцениваемых элементов матрицы связи не должно превышать $M - 1.$

Входящий в (9) неизвестный вектор ${\mathbf{c}}$ может быть найден лишь с точностью до произвольного комплексного множителя, поэтому для однозначного определения ${\mathbf{c}}$ можно наложить на него дополнительное ограничение ${{\left\| {\mathbf{c}} \right\|}^{2}} = 1$ (исключающее тривиальное решение ${\mathbf{c}} = 0$). В этом случае минимум квадратичной формы ${{{\mathbf{c}}}^{ + }}{\mathbf{Q}}{\kern 1pt} {\kern 1pt} ({\mathbf{\theta }}){\kern 1pt} {\kern 1pt} {\mathbf{c}}$ реализуется при условии совпадения ${\mathbf{c}}$ с собственным вектором, отвечающим наименьшему собственному значению ${{{{\lambda }}}_{{\min }}}\{ {\kern 1pt} {\kern 1pt} Q{\kern 1pt} {\kern 1pt} ({\mathbf{\theta }})\} $ матрицы $Q{\kern 1pt} ({\mathbf{\theta }}).$ В результате пространственные положения источников могут быть найдены из следующего критерия:

(10)
$\begin{gathered} \{ {\mathbf{\hat {\theta }}}\} _{{j = 1}}^{J} = \arg \mathop {\max }\limits_{\mathbf{\theta }}^ {\kern 1pt} {{P}_{{{\text{RARE}}}}}({\mathbf{\theta }}), \\ {{P}_{{{\text{RARE}}}}}({\mathbf{\theta }}) = {1 \mathord{\left/ {\vphantom {1 {{{{{\lambda }}}_{{\min }}}\{ Q{\kern 1pt} {\kern 1pt} ({\mathbf{\theta }})\} }}} \right. \kern-0em} {{{{{\lambda }}}_{{\min }}}\{ Q{\kern 1pt} {\kern 1pt} ({\mathbf{\theta }})\} }}. \\ \end{gathered} $

Метод оценивания (10) аналогичен алгоритму пониженного ранга RARE [9, 10 ] , используемому в технике частично калиброванных АР, содержащих $P$ подрешеток. При этом роль векторов, характеризующих амплитудно-фазовую калибровку этих подрешеток, играют векторы, определяемые соотношением (4). Важно подчеркнуть, что применение (10) позволяет решить обратную задачу без использования трудоемкой процедуры одновременного поиска как искомых координат, так и неизвестных элементов соответствующей матрицы связи.

При практической реализации данного способа локализации в качестве матрицы ${\mathbf{V}}({\mathbf{\theta }})$ (вследствие неполной информации о канале распространения) используется некоторая оценочная матрица ${{{\mathbf{V}}}_{0}}({\mathbf{\theta }}),$ рассчитываемая для номинальных акустических характеристик волновода. Однако в изменчивых и всегда не полностью известных условиях морской среды существующее рассогласование между $V({\mathbf{\theta }})$ и ${{V}_{0}}({\mathbf{\theta }})$ приводит к значительному ухудшению работоспособности приведенного метода. Ниже мы построим робастную версию алгоритма, позволяющую повысить устойчивость процедуры оценивания и частично скомпенсировать эффект подобного несоответствия.

ПОСТРОЕНИЕ АДАПТИВНОЙ ВЕРСИИ АЛГОРИТМА ПОНИЖЕННОГО РАНГА

При построении адаптивной процедуры RARE, основанной на наихудшем сценарии приема, будем предполагать, что истинная матрица ${\mathbf{V}}({\mathbf{\theta }})$ отличается от ${{V}_{0}}({\mathbf{\theta }})$ некоторой ошибкой, норма Фробениуса которой ограничена заданной величиной: $\left\| {V({\mathbf{\theta }}) - {{V}_{0}}({\mathbf{\theta }})} \right\|_{F}^{2} \leqslant \,\varepsilon $, где $\varepsilon $ – положительный параметр регуляризации. Адаптация к условиям рассогласования состоит в нахождении робастной матрицы ${\mathbf{V}}({\mathbf{\theta }},\varepsilon ),$ удовлетворяющей указанному ограничению и обеспечивающей минимум целевой функции ${{{\mathbf{c}}}^{ + }}{{{\mathbf{V}}}^{ + }}({\mathbf{\theta }}){\mathbf{\Psi }}_{{\mathbf{n}}}^{ + }{{{\mathbf{\Psi }}}_{{\mathbf{n}}}}{\mathbf{V}}({\mathbf{\theta }}){\kern 1pt} {\kern 1pt} {\mathbf{c}}$ для всех возможных значений нормированных векторов ${\mathbf{c}}{\kern 1pt} :$

(11)
$\begin{gathered} \mathop {\min }\limits_{\mathbf{V}} \left\{ {{{{\mathbf{c}}}^{ + }}{{{\mathbf{V}}}^{ + }}({\mathbf{\theta }}){{{\mathbf{\Psi }}}_{{\mathbf{n}}}}{\mathbf{\Psi }}_{{\mathbf{n}}}^{ + }{\mathbf{V}}({\mathbf{\theta }}){\kern 1pt} {\mathbf{c}}} \right\}\,\,\,\,{\text{при}} \\ \left\| {V({\mathbf{\theta }}) - {{V}_{0}}({\mathbf{\theta }})} \right\|_{F}^{2} \leqslant \varepsilon . \\ \end{gathered} $

Для исключения в (11) тривиального решения $V({\mathbf{\theta }}) = 0$ необходимо наложить дополнительное условие на допустимую норму ожидаемой матрицы ${{V}_{0}}({\mathbf{\theta }})$: $\left\| {{{V}_{0}}({\mathbf{\theta }})} \right\|_{F}^{2} > {{\varepsilon }}$.

Для нахождения искомой матрицы $V$ составим функцию Лагранжа

$L(V,{{\nu }}) = {{c}^{ + }}{{V}^{ + }}{{{\mathbf{\hat {\Psi }}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + }Vc + {{\nu }}\left( {\left\| {V - {{V}_{0}}} \right\|_{F}^{2} - {{\varepsilon }}} \right)$
с неопределенным вещественным множителем, удовлетворяющим условию ${{\nu }} > 0.$ (Здесь аргумент ${\mathbf{\theta }}$ опущен для краткости записи.) Дифференцируя $L$ по $V$ (при фиксированном $\nu $) и приравнивая результат к нулю, получим

(12)
${{{\mathbf{\hat {\Psi }}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + }VD + {{\nu }}V - {{\nu }}{\kern 1pt} {{V}_{0}} = 0,\,\,\,\,\,D = {\mathbf{c}}{{{\mathbf{c}}}^{ + }}.$

Умножая (12) справа на матрицу $D$ и учитывая, что (в силу ${{\left\| c \right\|}^{2}} = 1$) ${{D}^{2}} = D,$ имеем

$\left[ {\left( {{{{{\mathbf{\hat {\Psi }}}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + } + \nu {\mathbf{I}}} \right){\mathbf{V}} - \nu {{{\mathbf{V}}}_{0}}} \right]{\mathbf{D}} = 0.$

Поскольку это равенство должно выполняться для всех рассматриваемых векторов, то приведенное уравнение будет заведомо удовлетворено, если

(13)
$\begin{gathered} V = {{\nu }}{{\left( {{{{{\mathbf{\hat {\Psi }}}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + } + {{\nu }}I} \right)}^{{ - 1}}}{{V}_{0}} \equiv \\ \equiv {{V}_{0}} - {{(1 + \nu )}^{{ - 1}}}{{{{\mathbf{\hat {\Psi }}}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + }{{V}_{0}}. \\ \end{gathered} $

При написании (13) использовано соотношение

${{\left( {{{{{\mathbf{\hat {\Psi }}}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + } + \nu {\mathbf{I}}} \right)}^{{ - 1}}} = \frac{1}{\nu }\left[ {{\mathbf{I}} - \frac{{{{{{\mathbf{\hat {\Psi }}}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + }}}{{\nu + 1}}} \right],$
являющееся следствием известной формулы обращения матриц
${{({\mathbf{A}} + {\mathbf{BC}})}^{{ - 1}}} = {{{\mathbf{A}}}^{{ - 1}}} - {{{\mathbf{A}}}^{{ - 1}}}{\mathbf{B}}\,{{({\mathbf{I}} + {\mathbf{C}}{{{\mathbf{A}}}^{{ - 1}}}{\mathbf{B}})}^{{ - 1}}}{\mathbf{C}}{{{\mathbf{A}}}^{{ - 1}}},$
в которой ${\mathbf{A}} = {{\nu }}{\kern 1pt} {\mathbf{I}},$ ${\mathbf{B}} = {{{\mathbf{\hat {\Psi }}}}_{{\mathbf{n}}}}$, ${\mathbf{C}} = {\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + }$ (и при этом ${\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + }{{{\mathbf{\hat {\Psi }}}}_{{\mathbf{n}}}} = {\mathbf{I}}$).

Множитель Лагранжа ${{\nu }}$ находится из условия, вытекающего из ограничения на норму матрицы рассогласования:

$\left\| {V - {{V}_{0}}} \right\|_{F}^{2} = {{(1 + {{\nu }})}^{{ - 2}}}\left\| {{{{{\mathbf{\hat {\Psi }}}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + }{{V}_{0}}} \right\|_{F}^{2} = {{\varepsilon }},$
откуда

${{\nu }} = - 1 + \sqrt {Tr{\kern 1pt} {\kern 1pt} {{({\mathbf{V}}_{0}^{ + }{{{{\mathbf{\hat {\Psi }}}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + }{{{\mathbf{V}}}_{0}})} \mathord{\left/ {\vphantom {{({\mathbf{V}}_{0}^{ + }{{{{\mathbf{\hat {\Psi }}}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + }{{{\mathbf{V}}}_{0}})} {{\varepsilon }}}} \right. \kern-0em} {{\varepsilon }}}} .$

В свою очередь, из неравенства $\nu > 0$ следует, что в процессе локализации параметр регуляризации ${{\varepsilon }}$ не должен превышать определенной величины

${{\varepsilon }} < {{F}_{0}}({\mathbf{\theta }}),\,\,\,\,{{F}_{0}}({\mathbf{\theta }}) = \sqrt {Tr{\kern 1pt} [{\mathbf{V}}_{0}^{ + }({\mathbf{\theta }}){{{{\mathbf{\hat {\Psi }}}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + }{{{\mathbf{V}}}_{0}}({\mathbf{\theta }})]{\kern 1pt} } .$

Таким образом, знание множителя $\nu (\varepsilon )$ при $\varepsilon < {{F}_{0}}\left( \theta \right)$ позволяет на основании (13) найти оптимальную матрицу ${\mathbf{V}}({\mathbf{\theta }},\varepsilon )$, рассчитать адаптивную матрицу

${\mathbf{Q}}({\mathbf{\theta }},{{\varepsilon }}) = {{{\mathbf{V}}}^{ + }}({\mathbf{\theta }},{{\varepsilon }}){{{\mathbf{\hat {\Psi }}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + }{\mathbf{V}}({\mathbf{\theta }},{{\varepsilon }}) = {{f}^{2}}({\mathbf{\theta }},{{\varepsilon }}){{{\mathbf{Q}}}_{0}}({\mathbf{\theta }}),$
где $f({\mathbf{\theta }},{{\varepsilon }}) = 1 - \sqrt {{{{\varepsilon }} \mathord{\left/ {\vphantom {{{\varepsilon }} {{{F}_{0}}({\mathbf{\theta }})}}} \right. \kern-0em} {{{F}_{0}}({\mathbf{\theta }})}}} ,$ а ${{{\mathbf{Q}}}_{0}}({\mathbf{\theta }}) = $ $ = {\mathbf{V}}_{0}^{ + }({\mathbf{\theta }}){{{\mathbf{\hat {\Psi }}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + }{{{\mathbf{V}}}_{0}}({\mathbf{\theta }}),$ и в итоге оценить искомые координаты источников путем поиска $J$ максимумов выходной мощности адаптивного процессора

(14)
$\begin{gathered} \{ {\mathbf{\hat {\theta }}}\} _{{j = 1}}^{J} = \arg \mathop {\max }\limits_{\mathbf{\theta }}^ {{P}_{{RARE}}}({\mathbf{\theta }},{{\varepsilon }}), \\ {{P}_{{RARE}}}({\mathbf{\theta }},{{\varepsilon }}) \equiv {1 \mathord{\left/ {\vphantom {1 {{{f}^{2}}}}} \right. \kern-0em} {{{f}^{2}}}}({\mathbf{\theta }},{{\varepsilon }}){{{{\lambda }}}_{{\min }}}\{ {{{\mathbf{Q}}}_{0}}({\mathbf{\theta }})\} . \\ \end{gathered} $

Отметим, что при $\varepsilon = 0,$ $f({\mathbf{\theta }},\varepsilon ) = 1$ и предложенный способ оценивания переходит в неадаптивную версию алгоритма.

РЕЗУЛЬТАТЫ СТАТИСТИЧЕСКОГО МОДЕЛИРОВАНИЯ

Приведем результаты статистического моделирования, иллюстрирующие работоспособность предложенного способа оценивания (14), и сравним его эффективность с традиционным алгоритмом MUSIC, не учитывющим эффекты взаимодействия мод.

Для определенности рассмотрим мелководный канал с характерной летней гидрологией (изображенной на рис. 1), в котором звуковое поле создается двумя некоррелированными источниками, расположенными на глубинах 45 и 75 м, и излучающими узкополосные сигналы с несущей частотой 250 Гц.

Рис. 1.

Профиль скорости звука, используемый при моделировании.

В рамках численного эксперимента глубина моря в месте постановки приемной АР бралась равной $H = 120$ м, дно моделировалось жидким поглощающим полупространством с плотностью ${{\rho }_{b}} = 1.8$ г/см3, скоростью звука ${{c}_{b}} = 1750$ м/с и коэффициентом поглощения в грунте $\beta = 0.13$ дБ/${{\lambda }}$, а при расчете ожидаемого модового вектора и матрицы модовой структуры в качестве геоакустических параметров использовались значения $H = 121$ м, ${{\rho }_{b}} = 1.75$ г/см3, ${{c}_{b}} = 1725$ м/с и $\beta = 0.1$ дБ/${{\lambda }}$. Полное число мод для рассматриваемых акустических характеристик канала и несущей частоты составляло 20. Подчеркнем, что для рассматриваемой постановки задачи заданными считались лишь номинальные параметры волновода в месте установки АР, при этом предполагалось, что акустические характеристики канала изменялись вдоль трассы распространения неизвестным образом.

При вычислениях в качестве ${{{\mathbf{С}}}_{j}}$ использовались тридиагональные матрицы вида

${{{\mathbf{C}}}_{j}} = \left( {\begin{array}{*{20}{c}} 1&{{{{{\alpha }}}_{j}}}&0& \cdots &0 \\ {{{{{\gamma }}}_{j}}}&1&{{{{{\alpha }}}_{j}}}& \ddots & \vdots \\ 0&{{{{{\gamma }}}_{j}}}&1& \ddots &0 \\ \vdots & \ddots & \ddots & \ddots &{{{{{\alpha }}}_{j}}} \\ 0& \cdots &0&{{{{{\gamma }}}_{j}}}&1 \end{array}} \right),\,\,\,\,j = 1,2,$
при этом для каждого из элементов ${{{{\alpha }}}_{j}}$ и ${{{{\gamma }}}_{j}}$ бралась одна из реализаций, генерируемых случайным образом в соответствии с формулами: ${{c}_{j}} = {{x}_{j}} + i{{y}_{j}},$ где ${{x}_{j}}$ и ${{y}_{j}}$ – равномерно распределенные случайные числа из интервала $[ - 0.5,0.5]$.

Предполагалось, что прием осуществлялся линейной вертикальной антенной (с центром на глубине 60 м), состоящей из 32 элементов, расположенных через 3 м. Дистанция между источниками и антенной составляла соответственно 15 и 10 км. В процессе моделирования вектора наблюдения методом Монте-Карло комплексные огибающие излученных сигналов и компоненты шума рассматривались как статистически независимые гауссовские случайные процессы, при этом сигнальная матрица задавалась в виде ${\mathbf{S}} = diag(\sigma _{1}^{2},\sigma _{2}^{2}),$ где $\sigma _{1}^{2}$ и $\sigma _{2}^{2}$ – уровни излучения. Входные отношения сигнал/шум, определяемые соотношениями

$SN{{R}_{j}} = \frac{{{{\sigma }}_{j}^{2}{{{\left\| {{\mathbf{g}}({{{\mathbf{\theta }}}_{j}})} \right\|}}^{2}}}}{{{{\sigma }}_{n}^{2}N}},\,\,\,\,j = 1,2,$
считались одинаковыми, т.е. $SN{{R}_{1}} = SN{{R}_{2}} = SNR.$ Поиск источника по дальности осуществлялся в диапазоне 1–20 км с шагом 50 м, а по глубине – в интервале 1–120 м с шагом 0.5 м.

На рис. 2 показана зависимость среднеквадратических ошибок (СКО) оценивания положения источников по дальности и глубине от входного отношения сигнал/шум $SNR.$Соответствующие ошибки рассчитывались по формулам

$\begin{gathered} СКО(\hat {r}) = \sqrt {{{{(KJ)}}^{{ - 1}}}\sum\limits_{k = 1}^K {\sum\limits_{j = 1}^J {{{{(\hat {r}_{j}^{{(k)}} - {{r}_{j}})}}^{2}}} } } , \\ СКО(\hat {z}) = \sqrt {{{{(KJ)}}^{{ - 1}}}\sum\limits_{k = 1}^K {\sum\limits_{j = 1}^J {{{{(\hat {z}_{j}^{{(k)}} - {{z}_{j}})}}^{2}}} } } , \\ \end{gathered} $
в которых $\hat {r}_{j}^{{(k)}}$ и $\hat {z}_{j}^{{(k)}}$ – оценки координат $j$-го источника для $k$-ой реализации вектора наблюдения, $J = 2,$ а общее число независимых реализаций $K$ бралось равным $1000.$ Выборочная ковариационная матрица формировалась по $L = 120$ временным отсчетам.

Рис. 2.

Среднеквадратичные ошибки оценивания координат источника (а) – по дальности и (б) – по глубине в зависимости от входного $SNR$ для рассматриваемых методов обработки.

Кривые 1 на рис. 2 отвечают стандартному методу MUSIC, не учитывающему связи между модами, а кривые 2 соответствуют предложенному методу (14) (при ${{\varepsilon }} = 0.2$). Пунктирной линией изображена нижняя граница Крамера–Рао (выражение для которой приведено в Приложении), показывающая потенциально достижимую точность оценки соответствующих координат. Для рассматриваемого численного эксперимента, как следует из рис. 2, применение адаптивного алгоритма (14) позволяет частично компенсировать эффекты взаимодействия мод и тем самым значительно повысить точность оценивания координат по сравнению с обычным алгоритмом MUSIC. Последний не в состоянии приблизиться к нижней границе Крамера–Рао даже при больших значениях $SNR.$

Одной из важных характеристик алгоритма является достигаемая с его помощью вероятность правильной локализации, определяемая как доля реализаций, для которых ошибки в совместном определении положений источников по дистанции и глубине не превосходят заданных значений $\delta {\kern 1pt} r$ и ${{\delta }}{\kern 1pt} z.$ В качестве оценки соответствующей вероятности бралась величина

$\begin{gathered} {{{\hat {P}}}_{{CL}}} = \frac{1}{K}\sum\limits_{k = 1}^K {{{p}_{k}}} , \\ {{p}_{k}} = \left\{ \begin{gathered} 1,\,\,\,\,если\,\,\,\left| {\hat {r}_{j}^{{(k)}} - {{r}_{j}}} \right| < {{\delta }}{\kern 1pt} r \hfill \\ и\,\,\,\,\left| {\hat {z}_{j}^{{(k)}} - {{z}_{j}}} \right| < {{\delta }}z,j = 1,2; \hfill \\ 0,\,\,\,\,в\,\,\,противном\,\,\,случае. \hfill \\ \end{gathered} \right. \\ \end{gathered} $

Для используемых методов оценивания на рис. 3a представлены результаты расчета указанной вероятности в зависимости от $SNR{\kern 1pt} ,$ при этом число временных отсчетов $L$ бралось равным 120. На рис. 3б показана зависимость ${{\hat {P}}_{{CL}}}$ от числа выборок $L$, по которым оценивается ковариационная матрица входного процесса, при $SNR = - 3$ дБ. При вычислениях значения ${{\delta }}{\kern 1pt} r$ и ${{\delta }}{\kern 1pt} z$ принимались соответственно равными 400 и 2 м. Видно, что соответствующая вероятность весьма чувствительна к выбору алгоритма. В частности, наихудшую эффективность демонстрирует традиционный метод MUSIC, рассчитанный на прием звукового поля, сформированного невзаимодействующими модами, который не в состоянии обеспечить гарантированной локализации источников для всех рассматриваемых значений $SNR$ и $L.$

Рис. 3.

Вероятность правильной локализации источника в зависимости от (а) – входного $SNR$ и (б) – объема выборки $L$ для рассматриваемых методов обработки.

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

ЗАКЛЮЧЕНИЕ

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

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

Работа выполнена при финансовой поддержке РНФ (грант № 20-19-00383).

Список литературы

  1. Baggeroer A.B., Kuperman W.A., Mikhalevsky P.N. An overview of matched field methods in ocean acoustics // IEEE J. Oceanic Eng. 1993. V. 18. P. 401–423.

  2. Сазонтов А.Г., Малеханов А.И. Cогласованная пространственная обработка сигналов в подводных звуковых каналах (Обзор) // Акуст. журн. 2015. Т. 61. № 2. 233–253.

  3. Robust Adaptive Beamforming / Eds. by Li J. and Stoica P. John Wiley & Sons, Inc., Hoboken, New Jersey, 2006. 422 p.

  4. Бреховских Л.М., Лысанов Ю.П. Теоретические основы акустики океана. Л.: Гидрометеоиздат, 1982. 264 с.

  5. Friedlander B., Weiss A. Direction finding in the presence of mutual coupling // IEEE Trans. Antennas Propag.1991. V. 39. № 3. P. 273–284.

  6. Elbir A.M. Direction finding in the presence of direction-dependent mutual coupling // IEEE Antennas Wireless Propag. Lett. 2017. V. 16. P. 1541–1544.

  7. Ge Q., Zhang Y., Wang Y. A low complexity algorithm for direction of arrival estimation with direction-dependent mutual coupling // IEEE Commun. Lett. 2020 V. 24. № 1. P. 90–94.

  8. Schmidt R.O. Multiple emitter location and signal parameter estimation // IEEE Trans. on Antennas and Propagation. 1986. V. 34. № 3. P. 276–280.

  9. Pesavento M., Gershman A.B., Wong K.M. Direction finding in partly calibrated sensor arrays composed of multiple subarrays // IEEE Trans. on Signal Process. 2002. V. 50. № 9. P. 2103–2115.

  10. See C.M.S., Gershman A.B. Direction-of-arrival estimation in partly calibrated subarray-based sensor arrays // IEEE Trans. on Signal Processing. 2004. V. 52. № 2. P. 329–338.

  11. Stoica P., Nehorai A. MUSIC, maximum likelihood, and Cramer‑Rao bound // IEEE Trans. Acoust. Speech, Signal Processing. 1989. V. 37. № 5. P. 720–741.

  12. Nehorai A., Paldi E. Vector-sensor array processing for electromagnetic source localization. // IEEE Trans. Signal Process. 1994. V. 42. № 2. P. 376–398.

Дополнительные материалы отсутствуют.