Акустический журнал, 2019, T. 65, № 4, стр. 540-550

ВВЕДЕНИЕ

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

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

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

В настоящей работе в рамках наихудшего сценария построен адаптивный алгоритм NM–MUSIC (“normal mode based MUSIC”), предназначенный для локализации источника в акустическом волноводе с неточно известными параметрами, учитывающий эффекты рассогласования между принятым модовым составом и его расчетной моделью. Представлены результаты сравнительного анализа эффективности данного способа оценивания с методом максимума правдоподобия и традиционным методом MUSIC, осуществляющим обработку в пространстве элементов антенны. Приведена экспериментальная апробация предложенного алгоритма, демонстрирующая его работоспособность в мелководной акватории Ладожского озера.

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

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

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

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

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

Здесь ${{\varphi }_{m}}({{z}_{n}})$ – собственные функции волновода в месте расположения $n$-го гидрофона, ${{\kappa }_{m}}$ и ${{\alpha }_{m}}$ – постоянная распространения и коэффициент затухания в грунте m-ой моды. В свою очередь, для модели дна в виде однородного жидкого поглощающего полупространства с плотностью ${{\rho }_{b}}$ и продольной скоростью звука ${{c}_{l}} = {{c}_{b}}(1 - i\alpha )$ модовые коэффициенты ${{\alpha }_{m}}$ могут быть рассчитаны по формуле (см., например, [24]):

Здесь $H$ – глубина канала, ${{k}_{0}}$ – опорное волновое число звуковой волны, ${{\rho }_{w}}$ – плотность воды, ${{n}_{b}}$ – акустический показатель преломления в дне, а параметр α характеризует поглощающие свойства донных осадков и связан с измеряемым коэффициентом затухания $\beta $ соотношением $\beta = 40\pi \lg e\alpha $ дБ/λ.

В рамках волнового описания для вектора отклика АР ${\mathbf{g}}({\mathbf{\theta }})$ имеем

Здесь ${\mathbf{a}}({\mathbf{\theta }})$ – вектор размерности $M \times 1,$ компонентами которого являются амплитуды мод, определяемые формулой (2), а ${\mathbf{U}}$ – матрица модовой структуры размерности $N \times M$ с элементами ${{U}_{{nm}}} = {{\varphi }_{m}}({{z}_{n}}).$

При $N > {\kern 1pt} \,M$ матричное соотношение (3) представляет собой переопределенную линейную систему уравнений относительно вектора ${\mathbf{a}}({\mathbf{\theta }}).$ Решение этой системы, минимизирующее норму невязки $\left\| {{\mathbf{g}} - {\mathbf{Ua}}({\mathbf{\theta }})} \right\|,$ дается известным выражением

в котором ${{{\mathbf{U}}}^{{\dag }}}$ – матрица, псевдообратная к ${\mathbf{U}}.$ Обратим внимание, что для сравнительно коротких АР матрица взаимной ортогональности мод ${{U}^{T}}{\mathbf{U}}$ оказывается вырожденной, поэтому для нахождения ${{({{U}^{T}}{\mathbf{U}})}^{{ - 1}}}$ необходима регуляризация, например, на основе ее сингулярного спектрального разложения.

В общем случае случайный вектор ${{{\mathbf{n}}}_{l}},$ фигурирующий в (1), является суммой двух компонент: пространственно белого шума и распределенной помехи, представляющей собой динамические шумы морской среды, которые (также как и полезный сигнал) формируются преимущественно модами дискретного спектра. В рамках такого сценария результирующая пространственная ковариационная матрица аддитивного шума записывается в виде

Здесь $\sigma _{w}^{2}$ – дисперсия собственного шума АР, ${{{\mathbf{I}}}_{N}}$ – единичная матрица $N \times N,$ $\sigma _{m}^{2}$ и ${{K}_{m}}$ – уровень модового шума и его нормированная ковариационная матрица размерности $M \times M,$ а ${{( \cdot )}^{ + }}$ и $\left\langle \cdot \right\rangle $ означают операции эрмитового сопряжения и статистического усреднения, соответственно. Отметим, что для расчета матрицы ${{K}_{m}}$ можно воспользоваться моделью, предложенной в [25] (см. также [24], глава IX).

При выполнении неравенства $N > {\kern 1pt} \,M > J$ с вычислительной точки зрения пространственную обработку удобно реализовать в модовом пространстве, в котором вектор наблюдения ${{{\mathbf{y}}}_{l}}$ связан с исходным вектором (1) следующим образом:

Ковариационная матрица соответствующего M‑мерного вектора ${{{\mathbf{y}}}_{l}}$ равна

где ${\mathbf{A}}({\mathbf{\theta }}) = [{\mathbf{a}}({{{\mathbf{\theta }}}_{1}}), \ldots ,{\mathbf{a}}({{{\mathbf{\theta }}}_{J}})]$ – матрица $M \times J,$ составленная из J модовых векторов, зависящих от положений источников, а ${\mathbf{S}} = \left\langle {{{{\mathbf{s}}}_{l}}{\mathbf{s}}_{l}^{ + }} \right\rangle $ – матрица $J \times J,$ описывающая взаимные корреляции излученных сигналов. В реальных ситуациях истинная матрица ${{{\mathbf{\Gamma }}}_{{\mathbf{y}}}}$ неизвестна, поэтому при практической реализации модовых алгоритмов обработки вместо нее используется оценочная матрица ${{{\mathbf{\hat {\Gamma }}}}_{{\mathbf{y}}}}$ размерности $M \times M,$ равная где ${{{\mathbf{\hat {\Gamma }}}}_{{\mathbf{x}}}}$ – выборочная матрица размерности $N \times N$ исходного вектора наблюдения.

Одним из наиболее распространенных способов локализации является метод MUSIC [26], являющийся обобщением метода Писаренко [27] на случай произвольной размерности шумового подпространства. При приеме сигнала на фоне коррелированных помех с известной ковариационной матрицей ${{{\mathbf{\Gamma }}}_{{\mathbf{n}}}}$ эта процедура основана на использовании информации, содержащейся в системе собственных векторов обобщенной задачи на собственные значения

Здесь ${{\lambda }_{1}}, \ldots ,{{\lambda }_{N}}$ – положительные собственные числа, пронумерованные в порядке убывания (т.е. ${{\lambda }_{1}} \geqslant {{\lambda }_{2}} \geqslant \cdots \geqslant {{\lambda }_{N}}$), а ${{{\mathbf{\hat {\psi }}}}_{1}}, \ldots ,{{{\mathbf{\hat {\psi }}}}_{N}}$ – отвечающие им собственные векторы, причем ${\mathbf{\hat {\psi }}}_{i}^{ + }{{{\mathbf{\Gamma }}}_{{\mathbf{n}}}}{{{\mathbf{\hat {\psi }}}}_{j}} = {{\delta }_{{ij}}}.$

Если ранг сигнальной матрицы ${\mathbf{S}}$ равен $J,$ то число максимальных собственных чисел уравнения (7) совпадает с числом источников. Соответствующие им собственные векторы $\{ {\mathbf{\hat {\psi }}}\} _{{i = 1}}^{J}$ принадлежат сигнальному подпространству, а остальные $N - J$ младших собственных векторов $\{ {\mathbf{\hat {\psi }}}\} _{{i = J + 1}}^{N}$ формируют шумовое подпространство. В случае, когда матрица ${{{\mathbf{K}}}_{{\mathbf{m}}}},$ описывающая взаимные корреляции мод, возбуждаемых распределенными шумовыми источниками, имеет полный ранг, равный $M,$ на $N - M$ векторов $\{ {\mathbf{\hat {\psi }}}\} _{{i = M + 1}}^{N}$ накладываются дополнительные условия вида (см., например [28])

Для рассматриваемого способа обработки искомые координаты источников находятся из условия ортогональности ожидаемого вектора ${\mathbf{g}}({\mathbf{\theta }})$ к векторам шумового подпространства

и процедура оценивания сводится к поиску $J$ минимумов целевой функции ${{{\mathbf{g}}}^{ + }}({\mathbf{\theta }}){{{\mathbf{\hat {\Pi }}}}_{{\mathbf{n}}}}{\mathbf{g}}({\mathbf{\theta }}),$ где ${{{\mathbf{\hat {\Pi }}}}_{{\mathbf{n}}}} = {{{\mathbf{\hat {\Psi }}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Psi }}}_{{\mathbf{n}}}^{ + },$ а ${{{\mathbf{\hat {\Psi }}}}_{{\mathbf{n}}}} = [{{{\mathbf{\hat {\psi }}}}_{{J + 1}}}, \cdots ,{{{\mathbf{\hat {\psi }}}}_{M}}]$ – матрица размерности $N \times \,\,(M - J).$ (При написании (9) учтено представление (3) и соотношение (8).)

Очевидно, что приведенный алгоритм (9) путем замен ${\mathbf{g}}({\mathbf{\theta }}) \to {\mathbf{a}}({\mathbf{\theta }}),$ $\;{{{\mathbf{\hat {\Gamma }}}}_{{\mathbf{x}}}} \to {{{\mathbf{\hat {\Gamma }}}}_{y}}$ и ${{{\mathbf{\Gamma }}}_{{\mathbf{n}}}} \to {\mathbf{\Sigma }}$ может быть непосредственно перенесен в модовое пространство, в котором $M$-мерные собственные векторы $\{ {\mathbf{\hat {\varphi }}}\} _{{i = 1}}^{M}$, формирующие сигнальное и шумовое подпространства, являются решениями обобщенной задачи

и удовлетворяют условиям ортогональности ${\mathbf{\hat {\varphi }}}_{i}^{ + }{\mathbf{\Sigma }}{{{\mathbf{\hat {\varphi }}}}_{j}} = {{\delta }_{{ij}}}.$ Отметим, что соответствующие векторы ${{{\mathbf{\hat {\varphi }}}}_{i}}$ связаны с решениями уравнения (7) соотношениями ${{{\mathbf{\hat {\varphi }}}}_{i}} = {{{\mathbf{U}}}^{T}}{{{\mathbf{\hat {\psi }}}}_{i}},$ $(i = 1, \ldots ,M).$

Модовый критерий MUSIC непосредственно следует из условия (9) при подстановке в него (3):

где ${{{\mathbf{\hat {\Pi }}}}_{{\mathbf{n}}}} = {{{\mathbf{\hat {\Phi }}}}_{{\mathbf{n}}}}{\mathbf{\hat {\Phi }}}_{{\mathbf{n}}}^{ + },$ а ${{{\mathbf{\hat {\Phi }}}}_{{\mathbf{n}}}} = {{{\mathbf{U}}}^{T}}{{{\mathbf{\hat {\Psi }}}}_{{\mathbf{n}}}} \equiv [{{{\mathbf{\hat {\varphi }}}}_{{J + 1}}}, \cdots ,{{{\mathbf{\hat {\varphi }}}}_{M}}]$ – матрица размерности $M \times (M - J),$ составленная из собственных векторов уравнения (10), отвечающих наименьшим собственным числам.

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

Ниже мы построим адаптивный алгоритм NM–MUSIC, предназначенный для локализации источника в акустическом волноводе с неточно известными параметрами. С этой целью введем в рассмотрение $M$-мерные векторы ${{{\mathbf{\tilde {\varphi }}}}_{i}}$, связанные с ${{{\mathbf{\hat {\varphi }}}}_{i}}$ равенством ${{{\mathbf{\hat {\varphi }}}}_{i}} = {{{\mathbf{\Sigma }}}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}}{{{\mathbf{\tilde {\varphi }}}}_{i}},$ и сформируем из них матрицы ${{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{s}}}}$ и ${{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{n}}}}$ вида ${{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{s}}}} = [{{{\mathbf{\tilde {\varphi }}}}_{1}}, \ldots ,{{{\mathbf{\tilde {\varphi }}}}_{J}}]$ и ${{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{n}}}} = [{{{\mathbf{\tilde {\varphi }}}}_{{J + 1}}}, \cdots ,\,{{{\mathbf{\tilde {\varphi }}}}_{M}}].$ Отметим, что векторы ${{{\mathbf{\tilde {\varphi }}}}_{i}}$ взаимно ортогональны (${\mathbf{\tilde {\varphi }}}_{i}^{ + }{{{\mathbf{\tilde {\varphi }}}}_{j}} = {{\delta }_{{ij}}}$) и являются собственными векторами матрицы ${{{\mathbf{\tilde {\Gamma }}}}_{{\mathbf{y}}}} = {{{\mathbf{\Sigma }}}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}}{{{\mathbf{\hat {\Gamma }}}}_{{\mathbf{y}}}}{{{\mathbf{\Sigma }}}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}}.$ (Последняя представляет собой выборочную матрицу модового вектора наблюдения ${{{\mathbf{\tilde {y}}}}_{l}},$ получаемого из (4) в результате операции выбеливания ${{{\mathbf{\tilde {y}}}}_{l}} = {{{\mathbf{\Sigma }}}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}}{{{\mathbf{y}}}_{l}}$.) В свою очередь, матрицы ${{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{s}}}}$ и ${{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{n}}}}$ удовлетворяют соотношениям ${\mathbf{\tilde {\Phi }}}_{{\mathbf{s}}}^{ + }{{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{s}}}} = {{{\mathbf{I}}}_{J}},$ ${\mathbf{\tilde {\Phi }}}_{{\mathbf{n}}}^{ + }{{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{n}}}} = {{{\mathbf{I}}}_{{M - J}}}$ и ${\mathbf{\tilde {\Phi }}}_{{\mathbf{s}}}^{ + }{{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{n}}}} = {\mathbf{0}},$ вытекающим из условия ортогональности векторов $\{ {\mathbf{\tilde {\varphi }}}\} _{{i = 1}}^{M}.$ Учитывая, что ${{{\mathbf{\hat {\Phi }}}}_{{\mathbf{n}}}} = {{{\mathbf{\Sigma }}}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}}{{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{n}}}},$ для проекционной матрицы ${{{\mathbf{\hat {\Pi }}}}_{{\mathbf{n}}}},$ фигурирующей в (11), получим ${{{\mathbf{\hat {\Pi }}}}_{{\mathbf{n}}}} = {{{\mathbf{\Sigma }}}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}}{{{\mathbf{\tilde {\Pi }}}}_{{\mathbf{n}}}}{{{\mathbf{\Sigma }}}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}},$ где ${{{\mathbf{\tilde {\Pi }}}}_{{\mathbf{n}}}} = {{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{n}}}}{\mathbf{\tilde {\Phi }}}_{{\mathbf{n}}}^{ + }.$ В результате критерий ¹(11) может быть переформулирован следующим образом:

где $\tilde {a}({\mathbf{\theta }}) = {{{\mathbf{\Sigma }}}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}}a({\mathbf{\theta }}).$ В этом случае положения $J$ максимумов выходной мощности ${{P}_{{NM - MU}}}$ служат оценками искомых координат источников.

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

ПОСТРОЕНИЕ АДАПТИВНОГО МОДОВОГО АЛГОРИТМА NM–MUSIC

При построении адаптивной процедуры NM–MUSIC, основанной на наихудшем сценарии приема, будем предполагать возможность контролируемого отклонения ожидаемого модового вектора ${{{\mathbf{\tilde {a}}}}_{0}}({\mathbf{\theta }})$ от истинного $\tilde {a}({\mathbf{\theta }})$: евклидова норма вектора рассогласования не должна превышать заданную величину, т.е. ${{\left\| {\tilde {a}({\mathbf{\theta }}) - {{{\tilde {a}}}_{0}}({\mathbf{\theta }})} \right\|}^{2}} \leqslant \varepsilon ,$ где $\varepsilon $ – положительный параметр регуляризации. Адаптация к неизвестным условиям приема состоит в нахождении нормированного робастного вектора $\tilde {a}({\mathbf{\theta }},\varepsilon ),$ обеспечивающего максимум выходной мощности (12) (или минимум целевой функции ${{{\mathbf{\tilde {a}}}}^{ + }}({\mathbf{\theta }}){{{\mathbf{\tilde {\Pi }}}}_{{\mathbf{n}}}}{\mathbf{\tilde {a}}}({\mathbf{\theta }})$) и удовлетворяющего указанному выше неравенству:

Условия нормировки ${{\left\| {{\mathbf{\tilde {a}}}({\mathbf{\theta }})} \right\|}^{2}} = {{\left\| {{{{{\mathbf{\tilde {a}}}}}_{0}}({\mathbf{\theta }})} \right\|}^{2}} = 1$ позволяют исключить тривиальное решение оптимизационной задачи и переписать (13) в виде

Здесь аргумент ${\mathbf{\theta }}$ опущен для краткости записи.

Для нахождения оптимального вектора ${\mathbf{\tilde {a}}}$ составим функцию Лагранжа

где $\mu $ и $\nu $ – неопределенные множители, удовлетворяющие условиям $\mu \geqslant 0$ и $\nu > 0.$ Дифференцируя $L$ по $\tilde {a}$ (при фиксированных $\mu $ и $\nu $) и приравнивая результат к нулю, для искомого модового вектора получаем

Для вычисления множителя Лагранжа $\mu $ воспользуемся соотношением

вытекающим из ограничения на норму вектора рассогласования, откуда

В свою очередь, множитель $\nu $ находится из условия ${{\left\| {{\mathbf{\tilde {a}}}} \right\|}^{2}} = 1$ при подстановке в него (14) совместно с (15):

Для определения $\nu $ воспользуемся известной формулой обращения матриц, согласно которой ${{({\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 {{{\mathbf{I}}}_{M}},$ ${\mathbf{B}} = {{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{n}}}},$ ${\mathbf{C}} = {\mathbf{\tilde {\Phi }}}_{{\mathbf{n}}}^{ + }$ и учитывая, что ${\mathbf{\tilde {\Phi }}}_{{\mathbf{n}}}^{ + }{{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{n}}}} = {{{\mathbf{I}}}_{{M - J}}},$ получим

где ${{{\mathbf{\tilde {\Pi }}}}_{{\mathbf{s}}}} = {{{\mathbf{I}}}_{M}} - {{{\mathbf{\tilde {\Pi }}}}_{{\mathbf{n}}}}$ – проекционная матрица на сигнальное подпространство. Подстановка (17) в (16) приводит к квадратному уравнению относительно $\nu $: где ${{V}_{0}} \equiv {{V}_{0}}({\mathbf{\theta }}) = {\mathbf{\tilde {a}}}{}_{0}^{ + }({\mathbf{\theta }}){{{\mathbf{\tilde {\Pi }}}}_{{\mathbf{n}}}}{{{\mathbf{\tilde {a}}}}_{0}}({\mathbf{\theta }})$ имеет смысл целевой функции неадаптивного алгоритма NM–MUSIC. Положительный корень уравнения (18) существует при $\varepsilon < 2[1 - \sqrt {1 - {{V}_{0}}({\mathbf{\theta }})} ]$ и равен

Знание $\nu $ позволяет на основании (14), (15) и (17) найти робастный модовый вектор

рассчитать целевую функцию ${{V}_{{NM - MU}}}({\mathbf{\theta }},\varepsilon )$ = = ${{{\mathbf{\tilde {a}}}}^{ + }}({\mathbf{\theta }},\varepsilon ){{{\mathbf{\tilde {\Pi }}}}_{{\mathbf{n}}}}{\mathbf{\tilde {a}}}({\mathbf{\theta }},\varepsilon )$ адаптивного алгоритма NM–MUSIC и в итоге оценить искомые координаты источников путем поиска максимумов выходной мощности ${{P}_{{NM - MU}}}({\mathbf{\theta }},\varepsilon )$ = ${1 \mathord{\left/ {\vphantom {1 {{{V}_{{NM - MU}}}({\mathbf{\theta }},\varepsilon )}}} \right. \kern-0em} {{{V}_{{NM - MU}}}({\mathbf{\theta }},\varepsilon )}}$ или из эквивалентного критерия

где $f({\mathbf{\theta }},\varepsilon ) \equiv f{\kern 1pt} ({{V}_{0}}({\mathbf{\theta }}),\varepsilon )$ = = $1 - {\varepsilon \mathord{\left/ {\vphantom {\varepsilon 2}} \right. \kern-0em} 2} + \sqrt {\varepsilon (1 - {\varepsilon \mathord{\left/ {\vphantom {\varepsilon 4}} \right. \kern-0em} 4})[1 - {{{{V}_{0}}({\mathbf{\theta }})]} \mathord{\left/ {\vphantom {{{{V}_{0}}({\mathbf{\theta }})]} {{{V}_{0}}({\mathbf{\theta }})}}} \right. \kern-0em} {{{V}_{0}}({\mathbf{\theta }})}}} .$ (При получении (20) учтено, что ${{{\mathbf{\tilde {\Pi }}}}_{{\mathbf{n}}}}{{{\mathbf{\tilde {\Pi }}}}_{{\mathbf{s}}}} = 0$ в силу взаимной ортогональности ${{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{s}}}}$ и ${{{\mathbf{\tilde {\Phi }}}}_{{\mathbf{n}}}}$.) При $\varepsilon = 0,$ $f({\mathbf{\theta }},\varepsilon ) = 1$ и целевая функция ${{V}_{{NM - NM}}}({\mathbf{\theta }},\varepsilon )$ совпадает с функцией ${{V}_{0}}({\mathbf{\theta }})$ неадаптивного модового алгоритма NM–MUSIC, и следовательно, обычный метод NM–MUSIC является частным случаем предложенного способа оценивания. При фиксированном значении ${{V}_{0}}({\mathbf{\theta }})$ величина $f({\mathbf{\theta }},\varepsilon )$ является убывающей функцией параметра $\varepsilon .$ (Для иллюстрации на рис. 1 показано поведение $f{\kern 1pt} ({{V}_{0}}({\mathbf{\theta }}),\varepsilon )$ в зависимости от ε при ${{V}_{0}}({\mathbf{\theta }}) = 0.18.$)

Оценка ${\mathbf{\hat {\theta }}},$ при которой целевая функция в (20) достигает минимума, удовлетворяет соотношению $V_{{NM - MU}}^{'}({\mathbf{\hat {\theta }}},\varepsilon ) = 0,$ где $V_{{NM - MU}}^{'}$ означает взятие градиента ${{V}_{{NM - MU}}}$ по ${\mathbf{\hat {\theta }}}.$ Вычисление приводит к результату

Очевидно, что при выполнении неравенства

минимизация ${{V}_{{NM - NM}}}({\mathbf{\theta }},\varepsilon )$ эквивалентна минимизации ${{V}_{0}}({\mathbf{\theta }}).$ Простой анализ показывает, что неравенство (21) будет удовлетворено при любых ${\mathbf{\theta }},$ если $f{\kern 1pt} ({{V}_{0}}({\mathbf{\theta }}),\varepsilon ) > 0$ или ${{\varepsilon }^{2}} - 4\varepsilon + 4{{V}_{0}}({\mathbf{\theta }}) < 0,$ откуда

Таким образом, при соблюдении (22) адаптивный алгоритм NM–MUSIC приводит к тем же оценкам координат, что и его неадаптивная версия. Однако при конечном $\varepsilon ,$ принадлежащим интервалу (22), наличие множителя ${{f}^{2}}({\mathbf{\theta }},\varepsilon )$ в критерии (20) заметно уменьшает целевую функцию ${{V}_{{NM - NM}}}({\mathbf{\theta }},\varepsilon ),$ обеспечивая резкий максимум выходной мощности в области локализации источника. Последнее означает, что робастная процедура обладает повышенной разрешающей способностью и позволяет значительно подавить ложные пики на выходе процессора, затрудняющие получение корректного решения обратной задачи.

При $\varepsilon > 2[1 - \sqrt {1 - {{V}_{0}}({\mathbf{\theta }})} ]$ множитель Лагранжа $\nu ,$ определяемый формулой (19), становится отрицательным. В этом случае решение исходной оптимизационной задачи (13) не зависит от $\varepsilon $ и дается выражением

Это решение обращает в ноль целевую функцию в (13), удовлетворяет условию нормировки и для рассматриваемых ε автоматически гарантирует выполнение ограничения на допустимую норму рассогласования. Очевидно, что в процессе поиска координат параметр регуляризации ε должен удовлетворять неравенству (22), поскольку в противном случае $V{}_{{NM - MU}}({\mathbf{\theta }},\varepsilon ) \equiv 0.$

Для реализации предложенного адаптивного модового алгоритма используется следующая последовательность операций.

1. Для номинальных параметров акустического волновода и заданной геометрии АР строятся матрица модовой структуры ${{{\mathbf{U}}}_{0}}$ и псевдообратная к ней (с использованием при необходимости процедуры регуляризации); в области поиска рассчитываются ожидаемые модовые векторы ${{a}_{0}}({\mathbf{\theta }}).$

2. По принятой реализации входного процесса с помощью формулы (6) (в которой матрица ${\mathbf{U}}$ заменена на ${{{\mathbf{U}}}_{0}}$) вычисляется выборочная ковариационная матрица ${{{\mathbf{\hat {\Gamma }}}}_{{\mathbf{y}}}}$ модового вектора наблюдения, и в соответствии с определением ${\mathbf{\hat {\Sigma }}} = {\mathbf{U}}_{0}^{{\dag }}{{{\mathbf{\Gamma }}}_{{\mathbf{n}}}}{{({\mathbf{U}}_{0}^{{\dag }})}^{T}}$ находится оценка шумовой матрицы в модовом представлении.

3. Формируется матрица ${{{\mathbf{\tilde {\Gamma }}}}_{{\mathbf{y}}}} = {{{\mathbf{\hat {\Sigma }}}}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}}{{{\mathbf{\hat {\Gamma }}}}_{{\mathbf{y}}}}{{{\mathbf{\hat {\Sigma }}}}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}}$ и проводится ее спектральное разложение. Найденные собственные векторы шумового подпространства используются для построения проекционной матрицы ${{{\mathbf{\tilde {\Pi }}}}_{{\mathbf{n}}}}.$

4. На основании (20) определяется адаптивная целевая функция ${{V}_{{NM - MU}}}({\mathbf{\theta }},\varepsilon )$ модового алгоритма NM–MUSIC, находятся ее минимумы, служащие оценками искомых координат.

В заключение этого раздела отметим, что при использовании адаптивного метода MUSIC, осуществляющего обработку в пространстве элементов АР, соответствующие координаты источников могут быть оценены из условия (20), в котором необходимо сделать замены

где ${{{\mathbf{g}}}_{0}}({\mathbf{\theta }})$ – ожидаемый вектор отклика АР, а ${{{\mathbf{\tilde {\Pi }}}}^{ \bot }}$ – проекционная матрица, составленная из собственных векторов матрицы ${\mathbf{\Gamma }}_{{\mathbf{n}}}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}}{{{\mathbf{\hat {\Gamma }}}}_{{\mathbf{x}}}}{\mathbf{\Gamma }}_{{\mathbf{n}}}^{{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. \kern-0em} 2}}},$ принадлежащих шумовому подпространству.

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

Приведем результаты статистического моделирования, иллюстрирующие работоспособность предложенного способа оценивания (20), и сравним его эффективность с адаптивным алгоритмом MUSIC, осуществляющим обработку принятого сигнала в пространстве элементов АР, а также с модовым методом максимума правдоподобия, реализованном в [29] применительно к морским условиям. Согласно последнему методу положения источников определяются из следующего критерия:

где $Tr\,( \cdot )$ означает след матрицы, а ${{{\mathbf{P}}}_{{{\mathbf{A}}({\mathbf{\theta }})}}} = {\mathbf{A}}({\mathbf{\theta }}){{[{{{\mathbf{A}}}^{ + }}({\mathbf{\theta }}){{{\mathbf{\hat {\Sigma }}}}^{{ - 1}}}{\mathbf{A}}({\mathbf{\theta }})]}^{{ - 1}}}{{{\mathbf{A}}}^{ + }}({\mathbf{\theta }}){{{\mathbf{\hat {\Sigma }}}}^{{ - 1}}}.$

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

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

Истинный профиль скорости звука, измеренный в данном регионе, показан на рис. 2 сплошной линией (пунктиром показан профиль, используемый при расчетах). В рамках численного эксперимента глубина моря принималась равной $H = 160$ м; дно моделировалось жидким поглощающим полупространством с плотностью ${{\rho }_{b}} = 1.8$ г/см³, скоростью звука ${{c}_{b}} = 1750$ м/с и коэффициентом поглощения в грунте $\beta = 0.13$ дБ/λ, а при расчете ожидаемого модового вектора в качестве номинальных геоакустических параметров использовались значения $H = 162.5$ м, ${{\rho }_{b}} = 1.75$ г/см$^{3}$, ${{c}_{b}} = 1735$ м/с и $\beta = 0.1$ дБ/λ. Полное число мод $M$ для рассматриваемых акустических характеристик канала и несущей частоты составляло 28.

Предполагалось, что прием осуществлялся линейной вертикальной антенной (с центром на глубине 80 м), состоящей из 32 элементов, расположенных через 3 м. Дистанция между источниками и антенной составляла соответственно 15 и 10 км. При моделировании комплексные огибающие излученных сигналов и компоненты шума рассматривались как статистически независимые гауссовские случайные процессы, при этом сигнальная матрица задавалась в виде ${\mathbf{S}} = diag(\sigma _{1}^{2},\sigma _{2}^{2}),$ где $\sigma _{1}^{2}$ и $\sigma _{2}^{2}$ – уровни излучения, а для ковариационной матрицы модового шума, следуя [28, 29], использовалось представление ${{{\mathbf{K}}}_{m}} = {{{\mathbf{I}}}_{M}}.$ Для определенности считалось, что уровень динамического шума $\sigma _{m}^{2}$ в 4 раза превосходит уровень $\sigma _{w}^{2}$ пространственно белого шума. Входные отношения сигнал/шум, определяемые соотношениями

считались одинаковыми, т.е. SNR₁ = SNR₂ = = SNR. Поиск источников по дальности осуществлялся в диапазоне $(5 - 18)$ км с шагом 10 м, а по глубине – в интервале $(1 - 141)$ м с шагом 1 м.

Для рассматриваемого сценария на рис. 3 построен спектр собственных значений матрицы ортогональности мод $U_{0}^{T}{{{\mathbf{U}}}_{0}}.$ Из приведенного графика видно, что из 28 компонент этого спектра первые 15 практически одинаковы, а последние 9 близки к нулю. Следовательно, соответствующая матрица $U_{0}^{T}{{{\mathbf{U}}}_{0}}$ плохо обусловлена и при ее обращении необходима регуляризация. Ниже при нахождении матрицы, псевдообратной к ${{{\mathbf{U}}}_{0}}$, использовалось сингулярное разложение, учитывающее 15 наибольших собственных чисел. (В этом случае матрица $U_{0}^{T}{{{\mathbf{U}}}_{0}}$ близка к единичной, а выделяемые моды ортогональны на апертуре антенны.)

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

где ${{\hat {r}}_{{q,j}}}$ и ${{\hat {z}}_{{q,j}}}$ – оценки координат j-го источника для q-ой реализации вектора наблюдения, а общее число независимых реализаций $Q$ бралось равным 1000.

Для используемых методов оценивания на рис. 4a представлены результаты расчета указанной вероятности в зависимости от $SNR{\kern 1pt} ,$ при этом выборочная ковариационная матрица формировалась по $L = 100$ временным отсчетам общей длительностью 60 секунд. На рис. 4б показана зависимость ${{\hat {P}}_{{CL}}}$ от числа выборок $L$, по которым оценивается ковариационная матрица входного процесса, при $SNR = 0$ дБ. При вычислениях значения ${\delta }r$и ${\delta }z$ принимались соответственно равными 400 и 2 м. Кривая 1 на рис. 4 отвечает методу максимума правдоподобия (23), а кривые 2 и 3 – робастным алгоритмам MUSIC (при $\varepsilon = 0.3$), осуществляющим обработку в пространстве элементов АР и модовом пространстве, соответственно. Для рассматриваемого численного эксперимента, как следует из рис. 4, обычный алгоритм MUSIC незначительно уступает в эффективности модовому аналогу, однако его вычислительные затраты в разы больше, чем у MN–MUSIC. Что касается метода максимума правдоподобия, использующего неадаптивные ожидаемые модовые векторы, то он не в состоянии обеспечить гарантированной локализации источника для всех рассматриваемых значений $SNR$ и $L$.

ЭКСПЕРИМЕНТАЛЬНАЯ АПРОБАЦИЯ МЕТОДА

Для верификации предложенного метода локализации были использованы экспериментальные данные, полученные в августе 2014 года на Ладожском озере в стационарных условиях. Профиль скорости звука в месте проведения работ показан на рис. 5. На глубину 12 м был опущен источник, работавший одновременно на пяти частотах в килогерцовом диапазоне. Прием осуществлялся на вертикальную антенну (с центром на глубине 10.1 м), состоящую из 96 элементов, расположенных эквидистантно с шагом 0.2 м. Дистанция между источником и антенной была равной 150 м. Входное отношение сигнал/шум составляло порядка 15 дБ. Точные значения параметров дна в рассматриваемой акватории априори неизвестны и при моделировании ожидаемого модового вектора считалось, что осадочные породы представляли собой илы с характерными значениями плотности 1.2 г/см$^{3}$ и скорости звука 1450 м/с. Коэффициент затухания в грунте $\beta $ при вычислениях брался 0.1 дБ/λ. Выборочная ковариационная матрица оценивалась по первым 300 отсчетам, взятым из 2-х минутного фрагмента записи. При расчетах поиск источника по дальности осуществлялся в диапазоне (0–300) м с шагом 1 м, а по глубине – в интервале (0–18) м с шагом 0.5 м.

Ниже мы приведем решение задачи локализации источника, работавшего на несущей частоте 5025 Гц²², полагая что полезный сигнал регистрируется на фоне белого шума. Для данной частоты и указанных акустических характеристик водоема полное число распространяющихся мод составляло 12, профили которых в месте установки антенны представлены на рис. 6а. Поскольку в описываемом эксперименте приемная антенна перекрывала весь канал, то при выделении модового состава из исходного вектора наблюдения не требовалась процедура регуляризации. Средние интенсивности разрешаемых мод (нормированные на суммарную интенсивность), оцениваемые как диагональные элементы выборочной матрицы ${{{\mathbf{\hat {\Gamma }}}}_{{\mathbf{y}}}},$ показаны на рис. 6б. При построении ожидаемого модового вектора исключались компоненты с номерами $m > 6,$ так что в соответствующих алгоритмах локализации учитывались только наиболее энергонесущие моды, в наименьшей степени взаимодействующие со дном.

На рис. 7а изображена нормированная (на максимальное значение) выходная мощность ${{P}_{{{\text{ML}}}}}({\mathbf{\theta }}),$ построенная согласно (23) с использованием ожидаемого модового вектора ${{{\mathbf{a}}}_{0}}({\mathbf{\theta }}),$ рассчитанного для номинальных параметров волновода. Для сравнения на рис. 7б и 7в показано поведение нормированной выходной мощности адаптивных процессоров MUSIC, осуществляющих обработку в модовом пространстве и в пространстве элементов антенны. (При расчетах параметр регуляризации ε в рассматриваемых адаптивных алгоритмах выбран 0.15 и 0.9, соответственно.)

Из приведенных графиков видно, что в первых двух случаях положение абсолютного максимума выходного эффекта наблюдается при $\hat {r} = 148$ м и $\hat {z} = 11$ м, что довольно близко к истинным значениям координат источника. Однако применение метода максимума правдоподобия приводит к существенному уширению основного пика (области локализации), а следовательно, к значительному снижению разрешающей способности АР. В то же время при использовании традиционного алгоритма MUSIC, как следует из рис. 7в, наблюдается появление дополнительного ложного пика при $\hat {r} = 118$ и $\hat {z} = 11$ м, уровень которого превышает уровень правильного. Это можно объяснить тем, что традиционный метод по умолчанию учитывает все моды волновода, включая высокие. Последнее обстоятельство не позволяет обеспечить устойчивость данного способа локализации к детерминированному рассогласованию, обусловленному неточным знанием геоакустических характеристик донных осадков. Отметим также, что адаптивная процедура MN–MUSIC без привлечения модовой фильтрации приводит к точно такому же некорректному решению, что и стандартный алгоритм MUSIC.

ЗАКЛЮЧЕНИЕ

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

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

Работа выполнена в рамках государственного заказа ИПФ РАН (тема № 0035–2014–0011).