Estimation of gas saturated sediments parameters in shallow water using vector receiver

详细

A scheme for estimating parameters of gas-saturated sediments in shallow waters based on signals from moving noise source is proposed and tested on experimental data. Advantages of joint use of information about pressure and vertical component of particle velocity in solving the considered problem are discussed.

全文:

ВВЕДЕНИЕ

Восстановление скорости звука в дне водоема является важной задачей как для акустического мониторинга акваторий [1], так и для построения моделей распространения акустических волн в волноводе, например, для экологического контроля уровней звука при изучении воздействия антропогенного шума на морские живые организмы [2]. Известно [3], что в мелководных пресноводных акваториях скорость звука в дне существенно меньше, чем скорость звука в воде. Это связано с тем, что донный грунт пресноводных водоемов, а также ряда шельфовых районов морей часто содержит газонасыщенные осадочные породы с некоторым количеством нерастворимых газов [4]. При наличии газовых пузырьков средняя плотность уменьшается незначительно, а сжимаемость увеличивается многократно, что и приводит к уменьшению скорости [3]. От изменения значений скорости звука в дне в существенной степени зависит величина пространственного затухания звукового поля [5]. Этот факт может быть использован для восстановления параметров волновода по данным в виде пространственного затухания различных составляющих акустического поля. Ранее такая возможность была продемонстрирована на основе рассмотрения пространственного затухания давления как в рамках численного моделирования, так и при обработке экспериментальных данных с цепочки гидрофонов [6]. В настоящей работе для решения этой задачи используются данные с одиночного комбинированного приемного модуля (КПМ), позволяющего проводить одновременные измерения давления и ортогональных компонент колебательной скорости [7, 8]. Следует отметить, что в настоящее время наблюдается повышенный интерес к использованию подобных модулей для решения задач геоинверсии [9].

Работа имеет следующую структуру. Вначале рассматривается помехоустойчивость оценок скорости звука в дне при использовании информации о пространственном затухании различных составляющих акустического поля. Далее описывается схема решения обратной задачи и приводятся результаты ее применения для восстановления параметров волновода при обработке экспериментальных данных, измеренных на гидроакустическом полигоне МГУ [5, 8]. Этот полигон достаточно хорошо изучен [6, 8, 10], характеризуется присутствием газонасыщенного дна, возможности восстановления параметров которого рассматриваются в настоящей работе. В заключении приводятся основные результаты и выводы работы.

ОЦЕНКА СКОРОСТИ ЗВУКА В ГАЗОНАСЫЩЕННОМ ДНЕ ПО ДАННЫМ О ПРОСТРАНСТВЕННОМ ЗАТУХАНИИ

При численном моделировании использовалась упрощенная модель волновода в виде изоскоростного водного слоя, лежащего на жидком однородном полупространстве со скоростью звука c_b, меньшей, чем в водном слое c. В случае с_b≤(0.1–0.4)c₀ выражение для потенциала скорости φ на горизонтальном удалении r между точками излучения и приема при временной зависимости $~\exp \left(i2\pi ft\right)$, f – частота в Гц, может быть записано в виде [10]:

$\begin{array}{c}\mathrm{\phi }\simeq -\frac{2\pi i}{H}\sum H_0^{\left(2\right)}\left(x_{l}r\right)\times \\ \times \sin \left(y_l{z}_s\right)\sin \left(y_l{z}_r\right),\end{array}$(1)

где $y_l=\frac{\pi l}{H}\left(1+i\frac{{\rho }_b{c}_b}{2\pi f{\rho }_{0}H}\right)$, $x_l=\sqrt{k_0^2-{\left(\frac{\pi l}{H}\right)}^2+{\left(\frac{{\rho }_b{c}_{b}l}{2{\rho }_{0}f{H}^2}\right)}^2-i\frac{\pi {\rho }_b{c}_b{l}^2}{\rho f{H}^3}}$ – комплексные вертикальная и горизонтальная компоненты волнового вектора для заданного номера моды l соответственно; суммирование проводится по номерам l распространяющихся мод; H – глубина водного слоя; $H_0^{\left(2\right)}$ – функция Ханкеля нулевого порядка второго рода; z_S, z_r – глубины погружения излучателя и приемника; ρ_b, c_b – плотность и скорость звука в подстилающей среде; ρ₀, c₀– плотность и скорость звука в воде; $k_0=2\pi f/c_0$; i – мнимая единица.

В качестве исходных данных для решения обратной задачи используется пространственное затухание давления P, вертикальной компоненты колебательной скорости V_z (P и V_z могут быть вычислены на основе j стандартными методами [3]) и вертикальной составляющей потока акустической мощности $W_z=\frac{1}{4}\left(P{V}_z^{*}+P^{*}{V}_z\right)$, “*” означает комплексное сопряжение. В качестве величины, характеризующей пространственное затухание, рассчитывается относительный уровень L(r): $L\left(r\right)=10\lg \left[\frac{{〈g^2(r,f)〉}_{\Delta f}}{{〈g^2(r_0,f)〉}_{\Delta f}}\right]$, где ${〈\mathrm{...}〉}_{\Delta f}$ обозначает усреднение по рассматриваемой полосе частот Δf, а в качестве функции g(r, f) используются значения $\left|P(r,f)\right|$, $\left|V_z(r,f)\right|$ или $\sqrt{\left|W_z(r,f)\right|}$, вычисленные для заданных r и f на основе (1); r₀ – минимальное из рассматриваемых расстояний между источником и приемником. Выбор исходных данных для решения обратной задачи в виде пространственного затухания связан с тем, что зачастую в ходе экспериментальных работ регистрируются не абсолютные значения составляющих акустического поля, а их относительные уровни, как правило, в децибелах (дБ). Тем самым информация о фазе теряется, однако сохраняется зависимость относительных уровней от расстояния до источника, то есть информация о пространственном затухании. Такая ситуация характерна, например, для проведения экологических измерений уровней антропогенных шумов с использованием сертифицированного оборудования [2].

Рассматриваемая обратная задача сводится к минимизации многомерной невязки $\epsilon \left(\overrightarrow{m}\right)=||L_{\text{obs}}-L\left(\overrightarrow{m}\right)||/||L_{\text{obs}}||$ между экспериментальными $L_{\text{obs}}$ и модельными $L\left(\overrightarrow{m}\right)$ данными, численно рассчитанными для волновода с параметрами $\overrightarrow{m}$; при вычислении $\epsilon \left(\overrightarrow{m}\right)$ используется L₂ норма. Вектор параметров $\overrightarrow{m}$ в общем случае может содержать произвольный набор характеристик как среды распространения (глубина волновода, скорость звука и плотность дна и др.), так и характеристик источника (глубина погружения, дальность, скорость перемещения и др.). За оценку ${\overrightarrow{m}}^{\text{est}}$ параметров $\overrightarrow{m}$ принимается тот вектор, который обеспечивает минимум невязки $\epsilon \left(\overrightarrow{m}\right)$: $\epsilon \left({\overrightarrow{m}}^{\text{est}}\right)=\underset{\overrightarrow{m}}{\min }\epsilon \left(\overrightarrow{m}\right)$. В текущем разделе рассматривается простейший случай, когда неизвестной является только скорость звука в дне c_b.

При численном моделировании использовались следующие параметры задачи, близкие к условиям эксперимента [5, 8]: H = 8 м, z_s = 0.3 м (приповерхностный источник), z_r = 7.5 м (приемный модуль располагается вблизи дна), ρ_b= 1700 кг/м³, ρ₀ = 1000 кг/м³, c₀ = 1500 м/с. Расстояние между источником и приемником к изменялось от 26 м до 52 м с шагом 0.5 м (r₀ = 26 м). Вычислялись уровни L(r) для давления P и вертикальных компонент V_z, W_z – L_P(r), $L_{V_Z}\left(r\right)$, L_W(r) соответственно. При восстановлении использовались исходные данные отдельно в виде L_P(r), $L_{V_Z}\left(r\right)$ и L_W(r). Дополнительно рассматривался одновременный учет пространственных затуханий L_P(r) и $L_{V_Z}\left(r\right)$ (в этом случае суммарная невязка вычислялась как среднеарифметическое невязок ε_р и ${\epsilon }_{V_z}$, полученных отдельно для L_P(r) и $L_{V_Z}\left(r\right)$ соответственно: ${\epsilon }_{P+V_z}=\left({\epsilon }_P+{\epsilon }_{V_z}\right)/2$), а также объединенная оценка для всех трех функций L_P(r), $L_{V_Z}\left(r\right)$, L_W(r) (в этом случае итоговая невязка ${\epsilon }_{P+V_z+W_z}=\left({\epsilon }_P+{\epsilon }_{V_z}+{\epsilon }_{W_z}\right)/3$). Для анализа помехоустойчивости получаемых оценок c_b к модельным данным $P(r,f)$ и $V_z(r,f)$ независимо на разных частотах f и дальностях r добавлялась нормально распределенная случайная шумовая помеха с нулевым средним и среднеквадратичным амплитудным отклонением, составляющим 30% от значений $P(r,f)$ и $V_z(r,f)$. В реальной ситуации отношение сигнал/помеха зависит от дальности и частотного диапазона, шумы могут быть коррелированными, то есть рассматриваемая модель шумов весьма упрощенная. Детальное рассмотрение возможных моделей шумов выходит за рамки настоящей работы. Примеры расчетов зависмостей $\mathrm{\epsilon }\left({с}_b\right)$ при использовании различных комбинаций исходных данных представлены на рис. 1а. Как видно на рис. 1а, невязка $\mathrm{\epsilon }\left({с}_b\right)$ имеет минимум вблизи истинного значения c_b = 170 м/с. Однако оценка скорости ${\mathrm{с}}_b^{\text{est}}$, соответствующая этому минимуму, не совпадает с истинным значением c_b, причем величина этого отличия разная при использовании различных исходных данных. Наиболее точные результаты удается получить при совместном использовании информации о пространственном затухании P, V_z, W_z. Были последовательно произведены вычисления невязок для 500 различных реализаций шума при постоянном его уровне в 30%. Усредненные по множеству реализаций шума значения невязок $〈{\mathrm{\delta c}}_b〉\equiv 〈\frac{\left|{\mathrm{c}}_{\mathit{b}}-{\mathrm{c}}_b^{\text{est}}\right|}{{\mathrm{c}}_b}〉$ представлены на рис. 1б. Как видно на рис. 1б, минимальное значение$〈{\mathrm{\delta c}}_b〉$ достигается при совместном использовании информации о пространственном затухании P, V_z, W_z. При этом пространственное затухание только давления P дает ошибку оценки скорости, более чем в два раза превышающую аналогичное значение, полученное при совместном рассмотрении P, V_z, W_z. Второе по возрастанию значение ошибки $〈{\mathrm{\delta c}}_b〉$ дает способ расчета невязки с использованием вертикальной компоненты потока акустической мощности W_z. Следует отметить, что зависимость $\mathrm{\epsilon }\left({\mathrm{с}}_b\right)$, полученная для W_z, имеет большие градиенты близи минимума (рис. 1а), по сравнению с одновременным использованием P, V_z, W_z, когда минимум невязки $\mathrm{\epsilon }\left({\mathrm{с}}_b\right)$ оказывается более пологим. Это является преимуществом при поиске глобального минимума многомерной невязки $\epsilon \left(\overrightarrow{m}\right)$ методами типа градиентного спуска по зашумленным данным. В следующем разделе при обращении экспериментальных данных будет использоваться информация о пространственном затухании W_z.

Рис. 1. Пример зависимостей ε(c_b), полученных в присутствии шумов (задавался уровень шума 30%), при использовании различных комбинаций исходных данных (а); усредненные по 500 реализациям шума значения относительных отклонений $\langle \delta c_b\rangle$ восстановленных и истинных значений скорости звука в дне (б).

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

Экспериментальные работы проводились на гидроакустическом полигоне МГУ, в акватории Клязьминского водохранилища [5, 8]. Комбинированный приемный модуль (рис. 2а) располагался в металлической рамке вблизи дна, а регистрируемые им сигналы передавались по кабелю на борт лодки, которая была заякорена вблизи места постановки (рис. 2б). Рассматривался отрезок записи, в котором присутствует сигнал от одиночного катера, проходящего мимо места постановки КПМ. Горизонтальное расстояние r между катером и КПМ определялось с помощью лазерного дальномера. По экспериментально зарегистрированным данным $P(r,t)$, $V_z(r,t)$ расчитывались значения $W_z(r,t)$ и вычислялся относительный уровень $L_{\text{obs}}\left(r\right)$: $L_{\text{obs}}\left(r\right)=10\lg \left[\frac{{〈W_z(r,t)〉}_{\mathrm{\Delta }t}}{{〈W_z(r_0,t)〉}_{\mathrm{\Delta }{t}_0}}\right]$, где ${〈\mathrm{...}〉}_{\mathrm{\Delta }t}$ обозначает усреднение по временному интервалу Dt, соответствующему дальности r; r₀ – траверзное расстояние, проходимое источником за время Δt₀.

Рис. 2. Использованный в эксперименте комбинированный приемный модуль (КПМ) (а); условная схема проведенных измерений на гидроакустическом полигоне МГУ с движущимся катером (б).

При решении обратной задачи рассматривался следующий набор неизвестных параметров, входящих в (1): $\overrightarrow{m}=\left\{c_b,{\rho }_b,H,z_s,z_r\right\}$. Характеристики водного слоя полагались известными ρ₀ = 1000 кг/м³, с₀ =1500 м/с; использовался частотный диапазон от 900 до 1000 Гц, в котором преобладал сигнал источника. Важным этапом решения обратной задачи является задание области поиска минимума невязки $\epsilon \left(\overrightarrow{m}\right)$, т.е. априорное задание областей определения искомых параметров $c_b,{\rho }_b,H,z_s,z_r$. Исходя из известных данных о полигоне, а также с учетом условий эксперимента были выбраны следующие интервалы значений: $c_b\in \left[50,500\right]$ м/с, ${\mathrm{\rho }}_b\in \left[\mathrm{1.1,}1.8\right]\cdot {10}^3$ кг/м³, $H\in \left[\mathrm{7,}9\right]$м, $z_s\in \left[\mathrm{0.1,}1\right]$ м, $z_r\in \left[\mathrm{7,}7.9\right]$ м. Для поиска минимума невязки использовался алгоритм “имитации отжига” (simulated annealing – SA) [11], зарекомендовавший себя при решении задачи геоинверсии [12]. Отличительной особенностью этого алгоритма является вероятностный характер поиска решения [11], что позволяет получать различные реализации значений параметров ${\overrightarrow{m}}^{\text{est}}$. В рассматриваемой работе алгоритм “имитации отжига” запускался 1000 раз при фиксированных экспериментальных данных $L_{\text{obs}}$. При каждом запуске достигалась сходимость итерационной процедуры алгоритма. В результате был получен ансамбль из 1000 различных реализаций параметров ${\overrightarrow{m}}^{\text{est}}$, каждый из которых описывает экспериментальные данные $L_{\text{obs}}$ с точностью не хуже 9%: $\epsilon \left({\overrightarrow{m}}^{\text{est}}\right)<0.09$. При этом минимальное значение невязки в полученном ансамбле оказалось равным ${\mathrm{\epsilon }}_{\min }\approx$0.05, а среднее арифметическое по ансамблю значение $〈\mathrm{\epsilon }〉\approx 0.06$. Эта точность приближенно характеризует ошибку, с которой были оценены $L_{\text{obs}}$, а также говорит о неидеальности используемой модели волновода, что не позволяет описать экспериментальные данные точно. При численном моделировании использовалась программная реализация алгоритма “имитации отжига” [13], параметры которой выбирались по умолчанию. Изменения параметров алгоритма практически не влияли на значения итоговых оценок. Используя полученный набор из 1000 реализаций параметров ${\overrightarrow{m}}^{\text{est}}$, можно рассчитать средние значения $\overline{\mathrm{x}}=\frac{1}{K}\sum _{k=1}^K{\mathrm{x}}_k$ и средние квадратические отклонения ${\overline{\mathrm{\sigma }}}_x=\sqrt{\frac{1}{K-1}\sum _{k=1}^K{\left(x_k-\overline{x}\right)}^2}$, К = 1000 по полученному ансамблю (здесь x обозначает c_b, ρ_b, H, z_s, z_r). Результаты этих вычислений представлены в таблице 1. Дальнейшее увеличение количества рассматриваемых реализаций ${\overrightarrow{m}}^{\text{est}}$ практически не влияло на значения получаемых оценок. Следует отметить, что детальное исследование статистических характеристик параметров $\overrightarrow{m}$ требует более аккуратного выбора элементов ансамбля их реализаций, что может быть выполнено алгоритмами типа МСМС (Monte Carlo Markov Сhain) [12], но выходит за рамки настоящей работы. Если в качестве решения обратной задачи рассматривать те параметры ${\overrightarrow{m}}_{\min }^{\text{est}}$, которые обеспечивают наименьшее значение невязки ${\epsilon }_{\min }=\epsilon \left({\overrightarrow{m}}_{\min }^{\text{est}}\right)$ в полученном ансамбле, то получаются значения, которые также представлены в таблице 1. Как видно из таблицы 1, значения ${\overrightarrow{m}}_{\min }^{\text{est}}$ близки к значениям $\overline{x}$ с учетом оцененных отклонений ${\overline{\mathrm{\sigma }}}_x$. Восстановленные параметры волновода, представленные в таблице 1, соответствуют физически реализуемым для рассматриваемого полигона [6, 8, 10], что указывает на достоверность полученных оценок.

Таблица 1. Результаты восстановления по экспериментальным данным: средние значения $\overline{\mathbf{x}}$ и средние квадратические отклонения ${\overline{\mathbf{\sigma }}}_{\mathbf{x}}$, здесь х обозначает ${\mathbf{c}}_{\mathbf{b}}\mathbf{,}\mathbf{ }{\mathbf{\rho }}_{\mathbf{b}}\mathbf{,}\mathbf{ }\mathit{H}\mathbf{,}\mathbf{ }{\mathit{z}}_{\mathbf{s}}\mathbf{,}\mathbf{ }{\mathit{z}}_{\mathbf{r}}$, полученные для ансамбля из 2500 реализаций искомых параметров $\overrightarrow{\mathit{m}}\mathbf{=}\mathbf{\{}{\mathbf{c}}_{\mathbf{b}}\mathbf{,}{\mathbf{\rho }}_{\mathbf{b}}\mathbf{,}\mathbf{H}\mathbf{,}{\mathbf{z}}_{\mathbf{s}}\mathbf{,}{\mathbf{z}}_{\mathbf{r}}\mathbf{\}}$, а также значения ${\overrightarrow{\mathit{m}}}_{\mathbf{min}}^{\mathbf{\text{est}}}$, соответствующие минимальной в рассматриваемом ансамбле невязке $\mathit{\epsilon }\mathbf{\left(}{\overrightarrow{\mathbf{m}}}_{\mathbf{min}}^{\mathbf{\text{est}}}\mathbf{\right)}\mathbf{\approx }\mathbf{0.05}$.

	$\overline{x}$	${\overline{\mathrm{\sigma }}}_{\mathrm{x}}$	${\overrightarrow{m}}_{\min }^{\text{est}}$
Скорость звука в дне, cb, м/с	1.5 ⋅ 102	0.5 ⋅ 102	1.7 ⋅ 102
Плотность дна, ρb, кг/м3	1.4 ⋅ 103	0.2 ⋅ 103	1.2 ⋅ 103
Толщина водного слоя, H, м	8.8	0.1	8.8
Глубина КПМ, zr, м	7.6	0.2	7.7
Глубина источника, zs, м	0.6	0.2	0.7

ЗАКЛЮЧЕНИЕ

В работе предложена новая схема экспериментальной оценки параметров мелководной акватории, отличительной особенностью которой является использование в качестве исходных данных для решения обратной задачи информации о пространственном затухании различных составляющих акустического поля. При обработке экспериментальных данных использовалось затухание вертикальной компоненты потока акустической мощности. Для практической реализации удобно, что в качестве источника может рассматриваться проходящее мимо судно, а не специальный контролируемый излучатель. Однако в предложенной схеме требуется дополнительно измерять дальность r до судна-источника, что не всегда удобно в эксперименте. Локализация шумового источника может быть осуществлена независимыми методами [14], что позволит еще более упростить обсуждаемый подход. Также представляет интерес рассмотрение различных частотных диапазонов при решении обратной задачи. Увеличение объема исходных данных в этом случае может благоприятно влиять на точность и помехоустойчивость получаемых оценок. Отмеченные аспекты относятся к перспективам дальнейших исследований. Полученные на текущем этапе исследований результаты обработки экспериментальных данных указывают на практическую реализуемость рассматриваемого подхода.

Исследование выполнено за счет гранта Российского научного фонда № 23-27-00271, https://rscf.ru/project/23-27-00271/.