Акустический журнал, 2024, № 2

Российская академия наук
АКУСТИЧЕСКИЙ 
ЖУРНАЛ
Том 70     № 2     2024     Март–Апрель
Журнал основан в январе 1955 г.
Выходит 6 раз в год 
ISSN: 0320-7919
Журнал издается под руководством 
Отделения физических наук РАН
Главный редактор
И.Б. Есипов
Редакционная коллегия:
Ю.И. Бобровницкий (зам. главного редактора), 
М.Л. Лямшев (отв. секретарь),
С.В. Егерев, В.Ю. Зайцев, А.А. Карабутов, 
Т.К. Козубская, В.Ф. Копьев, А.И. Коробов, 
А.И. Малеханов, М.А. Миронов,
В.Г. Петников, Е.В. Чарная
Редакционный совет:
Ю.В. Гуляев, С.Н. Гурбатов, В.А. Зверев , С.А. Никитов, 
Л.А. Островский, О.В. Руденко, А.П. Сарвазян,
Б.Н. Четверушкин
Зав. редакцией В.А. Гусев
Научн. редакторы В.А. Гусев, А.М. Романовская
Адрес редакции: 119991 Москва, Ленинские горы, физический факультет МГУ
 
Тел.: (495) 939-29-18; E-mail: acoust-journal@physics.msu.ru
Москва
ФГБУ «Издательство «Наука»
© Российская академия наук, 2024
© Редколлегия “Акустического журнала” 
     (составитель), 2024

СОДЕРЖАНИЕ
Том 70, номер 2, 2024
КЛАССИЧЕСКИЕ ПРОБЛЕМЫ  
ЛИНЕЙНОЙ АКУСТИКИ И ТЕОРИИ ВОЛН
Корреляционный итерационный метод акустической томографии  
с некогерентными источниками поля
К. В. Дмитриев	
143
Внутренние симметричные волны Лэмба для больших фазовых скоростей
В. В. Мокряков	
156
НЕЛИНЕЙНАЯ АКУСТИКА
Влияние 3D-печати на упругие свойства нитевидных образцов полимера ABS
А. Б. Володарский, А. И. Кокшайский, Н. И. Одина, А. И. Коробов, Е. С. Михалев	
167
ФИЗИЧЕСКАЯ АКУСТИКА
Акустическая диагностика подводных выбросов, распространяющихся в виде многофазной струи
И. К. Гималтдинов, М. В. Столповский, Е. Ю. Кочанова	
174
Теплоемкость и особенности фононного спектра  
монокристаллов твердых растворов иттрий-лютециевых алюмогранатов
С. А. Никитов, А. В. Таранов, Е. Н. Хазанов, Е. В. Чарная, М. В. Лихолетова, Е. В. Шевченко	
180
Исследование частотного диапазона работы пьезоэлектрического преобразователя  
акустооптического фильтра электрическим и оптическим методами
Н. В. Поликарпова, В. Э. Пожар	
186
Компенсация аберраций при фокусировке ультразвука через череп на основе данных КТ и МРТ
Д. Д. Чупова, П. Б. Росницкий, О. В. Солонцов, Л. Р. Гаврилов, В. Е. Синицын, Е. А. Мершина,  
О. А. Сапожников, В. А. Хохлова	
193
АКУСТИКА ОКЕАНА. ГИДРОАКУСТИКА
Оценка координат подвижных морских роботов с использованием векторно-скалярных антенн,  
стационарно установленных в глубоком море
Г. М. Глебова, Г. А. Жбанков, Г. Н. Кузнецов	
206
Определение дистанции до подводного источника  
в условиях дальних зон акустической освещенности
И. Е. Лободин, А. И. Машошин	
217
Влияние слаборасходящегося акустического пучка  
на формирование пространственно-временной структуры импульсных сигналов  
в подводном звуковом канале
Ю. В. Петухов, Е. Л. Бородина	
225
Моделирование коэффициента усиления вертикальной антенны  
в мелководном волноводе со взволнованной поверхностью 
М. А. Раевский, В. Г. Бурдуковская	
232

АТМОСФЕРНАЯ И АЭРОАКУСТИКА
Характеристики скалярного частотно-волнового спектра пристеночных пульсаций давления  
в безградиентном турбулентном пограничном слое
Е. Б. Кудашев, Л. Р. Яблоник	
244
АКУСТИКА СТРУКТУРНО НЕОДНОРОДНЫХ ТВЕРДЫХ СРЕД.  
ГЕОЛОГИЧЕСКАЯ АКУСТИКА
Экспериментальное исследование медленной релаксации скорости звука в карбонатной породе
А. В. Лебедев, С. А. Манаков	
253
ФИЗИЧЕСКИЕ ОСНОВЫ ТЕХНИЧЕСКОЙ АКУСТИКИ
Исследование качества термодиффузионной сварки кристаллов  
в дисковом оптическом элементе оптоакустическим методом
В. В. Казаков, И. Б. Мухин, А. А. Курников, П. В. Субочев	
273
ПИСЬМА В РЕДАКЦИЮ
К вопросу о звучании речи и фортепиано 
В. А. Зверев , А. И. Малеханов	
283
ИНФОРМАЦИЯ
Система цикло-пролонгированного информационного обеспечения в области акустики
А. Б. Горшков, В. Г. Шамаев	
289
Памяти Виталия Анатольевича Зверева (03.11.1924—06.03.2024)	
297

АКУСТИЧЕСКИЙ ЖУРНАЛ,  2024, том 70, № 2,  с.  143–155
КЛАССИЧЕСКИЕ ПРОБЛЕМЫ ЛИНЕЙНОЙ АКУСТИКИ 
И ТЕОРИИ ВОЛН
УДК 534.2:517.4
КОРРЕЛЯЦИОННЫЙ ИТЕРАЦИОННЫЙ МЕТОД  
АКУСТИЧЕСКОЙ ТОМОГРАФИИ  
С  НЕКОГЕРЕНТНЫМИ ИСТОЧНИКАМИ ПОЛЯ
© 2024 г. К. В. Дмитриев*
Московский государственный университет имени М.В. Ломоносова, физический факультет, ГСП-1,  
Ленинские горы, Москва, 119991 Россия
*е-mail: presentatio@mail.ru
Поступила в редакцию 01.09.2023 г.
После доработки 14.11.2023 г.
Принята к публикации 19.12.2023 г.
Предложен метод восстановления акустических параметров среды с помощью итерационной 
обработки матриц когерентности акустического поля случайных источников, для части из которых известна их плотность мощности. Обсуждаются возможности повышения устойчивости 
и ускорения сходимости метода. Проводится сравнение результатов восстановления с функционально-аналитическим подходом, основанным на обработке амплитуды рассеяния.
Ключевые слова: обратная задача рассеяния, метод шумовой интерферометрии, уравнение Липпмана–
Швингера, борновское приближение, итерационные алгоритмы, функционально-аналитические 
алгоритмы
DOI: 10.31857/S0320791924020013, EDN: YNXAGE
1. ВВЕДЕНИЕ
Пассивная акустическая томография – это метод 
определения свойств исследуемой среды, не использующий активные источники звука. Он основан 
на регистрации и последующей обработке акустического излучения от присутствующих в среде источников. Задачи такого рода возникают, например, 
в медицинской томографии [1–3], где акустическое 
поле порождается тепловым движением частиц; 
в геоакустике [4, 5], где его источниками служат 
микросейсмы; в гидроакустике [6, 7], где оно включает естественные и антропогенные шумы различной природы, в гелиосейсмологии [8, 9]. Отсутствие 
активных источников звука, во-первых, существенно удешевляет такие исследования, а во-вторых, 
является более предпочтительным с точки зрения 
вопросов экологии (в гео- и гидроакустике) или 
 
уменьшения воздействия на организм (в медицине). Вместе с тем, методы пассивной томографии 
оказываются сложнее, чем методы активной томографии. Это связано с необходимостью оценки 
не только искомых свойств среды, но и описания 
создающих акустическое поле источников. Такая 
комбинированная обратная задача излучения-рассеяния представляется сильно недоопределенной, 
и для ее успешного решения целесообразно наложить те или иные дополнительные ограничения 
на класс возможных источников поля и допустимые диапазоны свойств среды.
Одним из методов пассивной акустической томографии является метод шумовой интерферометрии. В его основе лежит предположение о том, 
что акустическое поле создается некогерентными 
источниками, пространственная плотность мощности которых изотропна по отношению к системе наблюдения. Тогда корреляционная функция 
сигналов, зарегистрированных в двух точках среды, 
оказывается пропорциональной разности запаздывающей и опережающей функций Грина этой 
среды с аргументами в виде координат выбранных двух точек. Корреляционную функцию можно 
определить на основе экспериментальных данных. 
Далее она используется для постановки и решения 
обратной задачи рассеяния. В качестве недостатка 
метода шумовой интерферометрии можно указать 
ограничения, накладываемые на область его применимости. Первым ограничивающим фактором 
является требование инвариантности описывающих среду уравнений по отношению к преобразованию обращения времени, которое нарушается 
в присутствии течений или поглощения звука. Эти 
вопросы обсуждались в работах [10–12]. Второй 
фактор связан с тем, что метод перестает работать, 
если мощность шума, приходящего с некоторого 
143

ДМИТРИЕВ
причем их координаты и мощность создаваемого 
шума известны. Изменение положения судов с течением времени обеспечивает перебор ракурсов 
облучения; при этом предполагается, что свойства 
среды за время наблюдений меняются слабо. В других случаях контролируемые источники могут быть 
созданы специально, что отвечает задаче активной 
или активно-пассивной томографии. Например, 
в медицине может идти речь о использовании собственного теплового акустического излучения биологических сред совместно с некогерентным излучением дополнительной “подсветки” [2, 3]. Сделав 
такое замечание, дальнейшее рассмотрение можно вести независимо от природы контролируемых 
источников.
2. КОРРЕЛЯЦИОННЫЙ 
ИТЕРАЦИОННЫЙ МЕТОД
Рассматривается рассеяние акустического поля 
на  неоднородности, занимающей внутри всего 
координатного пространства IRd  размерности d  
область ℜ конечных размеров. Вне ℜ находится однородная фоновая среда. Предполагается, 
что сигналы сосредоточены в узкой полосе частот 
с центральной частотой ω0 , и временнáя зависимость выбирается в виде ~ exp(
)
−i
t
ω0 . Координатные зависимости скорости звука c( )
r  и амплитудного коэффициента поглощения α( )
r  описываются с помощью комплексного волнового числа 
k
c
i
( )
( )
( )
r
r
r
≡
+


. Снаружи неоднородности, 
т.е. внутри фоновой среды, оно известно: k
k
( )
r ≡
0  
при r ∈
ℜ
IRd \
. Кроме этого, считается известной 
функция ka( )
r , которая равна k0  вне неоднородности, а внутри нее выражает априорную оценку 
свойств этой неоднородности. Такое предположение не снижает общность рассмотрения, т.к. в отсутствие дополнительной информации можно положить k
k
a( )
r ≡
0  во всем пространстве.
Потенциал ϕ( )
r  акустического поля удовлетворяет уравнению Гельмгольца:
	
∆ϕ
ϕ
( )
( ) ( )
( ),
r
r
r
r
+
=
k
F
2
 r ∈IRd ,	
(1)
где F( )
r  – источники поля. Его решение, с учетом 
условия излучения на бесконечности, имеет вид 
 r r
,
,
′ ∈IRd 	
(2)
	
ϕ( )
( ,
) ( )
,
r
r r
r
r
=
′
′
′
∫G
F
d
d
IR
где G( ,
)
r r′  – запаздывающая функция Грина рассматриваемой неоднородной среды. Аналогичное 
уравнение Гельмгольца можно записать для поля 
ϕa( )
r , которое создают те  же самые источники 
в среде с волновым числом ka( )
r :
направления, выше, чем с других. Поскольку такая 
анизотропия шума встречается часто, это требует 
осторожности при практическом применении.
В случае, когда анизотропия шумового поля вызвана наличием небольшого числа (по сравнению 
с числом элементов приемной антенной решетки) 
дополнительных источников, удается “восстановить” работоспособность метода шумовой интерферометрии [13]. Для этого корреляционной обработке подвергаются все доступные в эксперименте 
пары сигналов, в результате чего в каждой относительно узкой полосе частот вычисляются их матрицы когерентности. Для каждой из них осуществляется корректировка максимальных собственных 
чисел так, чтобы минимизировать дисперсию диагональных элементов итоговых матриц после коррекции. Схожая процедура описывалась в [14, 15], 
где она использовалась для улучшения видимости 
слабых источников сигнала на фоне более мощных. Физически ее применение означает такую 
фазировку приемной антенной решетки, при которой сигналы из области пространственной локализации дополнительных источников, создающих 
анизотропную часть шумового поля, подавляются 
до среднего уровня шумов. В результате получается 
набор корреляционных функций, близких по значениям к искомым разностям функций Грина.
В случае, когда число дополнительных источников соизмеримо или превосходит количество элементов приемной антенной решетки, описанная 
процедура коррекции собственных чисел не приводит к  успеху. Тогда предлагается отказаться 
от использования метода шумовой интерферометрии и рассмотреть задачу акустической томографии в следующей постановке, которая является 
предметом настоящей статьи. Предполагается, что 
шумовое поле в среде создается стационарными 
фоновыми источниками с неизвестной пространственной плотностью мощности I( )( ).
0 r  Дополнительно к ним в среде могут появляться независимые контролируемые некогерентные источники 
с известной пространственной плотностью мощности I n
( )( ).
r  Термин “контролируемые” понимается 
в том смысле, что в эксперименте доступны как 
сигналы, записанные при наличии лишь фоновых 
источников I( )( ),
0 r  так и сигналы, записанные при 
наличии источников обоих типов, т.е. с пространственной плотностью мощности I
I n
( )
( )
( )
( ).
0 r
r
+
 
Верхний индекс n,
, который в дальнейшем будет 
называться ракурсом облучения, при этом означает 
номер пространственной конфигурации контролируемых некогерентных источников, и при n = 0  
такие источники отсутствуют.
В некоторых случаях описанная постановка задачи может быть формально отнесена к пассивной томографии. Например, в приложениях гидроакустики может рассматриваться акустическое 
поле в акватории, по которой перемещаются суда, 
	
∆ϕ
ϕ
a
a
a
( )
( )
( )
( ),
r
r
r
r
+
=
k
F
2
 r ∈IRd .	
(3)
АКУСТИЧЕСКИЙ ЖУРНАЛ
том 70
№ 2
2024

	
Корреляционный итерационный метод акустической томографии   
145
H0
1
( )()
⋅  – функция Ханкеля первого рода нулеВ этом случае вводится запаздывающая функция Грина Ga( ,
)
r r′ , и решение уравнения (3) име. Если добавить 
ет вид ϕa
a
вого порядка. В случае, когда априорная оценка ka( )
r  близка к  истинному значению k( )
r , 
d
IR
( )
( ,
) ( )
r
r r
r
r
=
′
′
′
∫G
F
d
G
G
( ,
)
( ,
)
r r
r r
′ ≈
′
a
, и для уравнения (5) справедливо 
борновское приближение [17, 18, 20, 21]:
в обе части уравнения (1) слагаемое ε
ϕ
( ) ( )
r
r , где 
ε( )
( )
( )
r
r
r
≡
−
k
k
a
	
G
G
G
G
d
( ,
)
( ,
)
( ,
) (
)
(
,
)
,
r r
r r
r r
r
r
r
r
′ ≈
′ +
′′
′′
′′
′
′′
2
2
 – функция рассеивателя относительно среды с волновым числом ka( ),
r
 получится уравнение (3) с источниками поля, равными 
ℜ
∫
a
a
a
ε
	
r r
,
,
′ ∈IRd  ′′ ∈ℜ
r
.	
(7)
F ( )
( ) ( )
r
r
r
+ ε
ϕ
. Его формальное решение с помощью функции Грина приводит к хорошо известному в квантовой теории поля уравнению Липпмана–
Швингера [16]. Оно часто служит основой для 
решения обратных задач и в акустике [17, 18]:
	
ϕ
ϕ
ε
ϕ
( )
( )
( ,
) ( ) ( )
,
r
r
r r
r
r
r
=
+
′
′
′
′
ℜ
∫
a
a
G
d
(4)
Далее рассматриваются акустические поля, порожденные шумовыми источниками F n
( )( )
r  при 
каждом n -м ракурсе облучения. Регистрация сигналов осуществляется с помощью NR  приемников. 
Каждый j -й из них можно задать, определив множество точек S j  его поверхности. Тогда принятые 
сигналы U j
n
( )  представляются в виде 
	
r ∈IRd , ′ ∈ℜ
r
. 	
n
n
	
U
d
j
( )
( )( )
,
= ∫ϕ
r
r 	
(8)
S j
Это уравнение отражает тот факт, что полное поле 
источника ϕ( )
r  складывается из падающего поля 
где ϕ( )( )
n r , r ∈IRd  – поле акустического потенциала при соответствующем ракурсе облучения. 
На  их основе формируется набор матриц когерентности, элементы которых определяются как 
ϕa( )
r , которое он создавал бы в среде с волновым 
числом ka( )
r , и рассеянного поля, выраженного интегралом по области неоднородности ℜ, где 
функция ε( )
r  отлична от нуля. Учитывая, что (4) 
выполнено при любой конфигурации источников 
поля, и полагая F( )
(
)
r
r
r
=
−′′
δ
 при произвольном 
Γij
n
i
n
j
n
U
U
( )
( )
( )
*
=
{
}
, где “*” означает комплекс′′
r , можно записать такое же по форме уравнение 
относительно функции Грина G( ,
)
r r′ :
	
G
G
G
G
d
( ,
)
( ,
)
( ,
) (
) (
,
)
,
r r
r r
r r
r
r
r
r
′ =
′ +
′′
′′
′′
′
′′
ℜ
∫
a
a
ε
(5)
ное сопряжение, а угловые скобки – усреднение 
по всем реализациям. Последовательно учитывая 
(8) и (2), а также принимая во внимание, что функции Грина являются детерминированными, и их 
можно вынести из-под операции усреднения, для 
матриц когерентности можно получить выражение
	
r r
,
,
′ ∈IRd  ′′ ∈ℜ
r
.	
( )
Γ
=
ij
n
Его можно использовать для вычисления функции 
(9)
*
( )
( , )
(
,
)
( ,
)
Γ
. 		
F
n
d
d
d
d G
G
=
′
′′
′
′′
∫
r
r
r
r
r r
r
r
r r
1
2
1
2
1
2
d
d
	
	
IR
IR
∫
∫
∫
S
S
j
i
Здесь ΓF
n
n
n
F
F
( )
( )
( )
*
( ,
)
( )
(
)
r r
r
r
1
2
1
2
≡
{
}
 – функция 
G( ,
)
r r′ , если известна функция Ga( ,
)
r r′ . В  свою 
очередь, функцию Ga( ,
)
r r′  можно определить 
из того же уравнения (5), выполнив замену переменных и рассматривая теперь в качестве априорной 
оценки однородную среду с волновым числом k0  
и функцией Грина G0( ,
)
r r′ :
	 G
G
G
G
d
a
0
0
a
a
( ,
)
( ,
)
( ,
)
(
)
(
,
)
,
r r
r r
r r
r
r
r
r
′ =
′ +
′′
′′
′′
′
′′
ℜ
∫
ε
(6)
	
r r
,
,
′ ∈IRd  ′′ ∈ℜ
r
.	
когерентности источников поля, расположенных 
в двух произвольных точках r r
1
2
,
∈IRd. Поскольку 
фоновые источники и контролируемые источники 
предполагаются независимыми, и, кроме того, контролируемые источники некогерентны, выполняется равенство
	


F
n
F
n
I
( )
( )
( )
( ,
)
( ,
)
( ) (
),
r r
r r
r
r
r
1
2
0
1
2
1
1
2
=
+
−

  
	
r r
1
2
,
,
∈IRd
Здесь εa
0
a
( )
( )
r
r
≡
−
k
k
2
2
 – функция рассеивателя для априорной оценки ka( )
r  относительно 
фоновой среды. Функция Грина G0( ,
)
r r′  известна аналитически [2, 17–19]. Например, в двумерном случае ( d = 2 ), который будет рассматриваться в дальнейшем на этапе численного модегде ΓF
( )( ,
)
0
1
2
r r
 и I n
( )( ) (
)
r
r
r
1
1
2
δ
−
 – функции когерентности фоновых и контролируемых источников 
(
)
( )
( ,
)
,
=
′ = −
−′
(
)
r r
r
r
 где 
лирования, G
i H
k
d
0
2
0
1
0
4
АКУСТИЧЕСКИЙ ЖУРНАЛ
том 70
№ 2
2024

ДМИТРИЕВ
функция Ga( ,
)
r r′  пересчитывается согласно (6). 
В  борновском приближении G
G
( ,
)
( ,
),
r r
r r
′ ≈
′
a
 
и поэтому
поля, соответственно. Тогда разность матриц когерентности при n -м и нулевом ракурсах облучения 
зависит лишь от мощности контролируемых источников I n
( )( )
r :
( )
( )
( )( ,
)
−
−
′
′′
′
′′ ≈
a
Γ
Γ
Γ
ij
n
ij
G
n
	
Γ
Γ
Γ
ij
n
ij
G
n
∫
∫
0
r
r
r r
d d
( )
( )
( )( ,
)
,
−
=
′
′′
′
′′
∫
∫
0
r r
r
r
	
(10)
S
S
d
d
j
i
S
S
j
i
	
	
где
( )
	
(13)


( )
r r
r r
′
′′ +
Γ
a
ε
( , )
( ,
)
G
n
a
.
d
d
d
G
r
r
r
r
≈
′
′′
	

G
( )
*
n
n
G
G
I
d
( )
*
( )
( ,
)
( , )
(
, )
( )
,
′
′′ ≡
′
′′
∫
r r
r r
r
r
r
r
( )
(
, )
( ,
)
r
r
r
r r
Γ
a
a
ε
G
n
G





ℜ
∫
∫
∫
S
S
j
i




+
′′
′
{
}


	
′
′′ ∈
r r
,
.
IRd 	
(11)
Величина ΓG
n
( )( ,
)
′
′′
r r
 представляет собой функцию 
когерентности полей контролируемых источников, 
которые измеряются в точках ′
r  и ′′
r . Ее можно 
представить с помощью уравнения (5) в виде суммы 
четырех слагаемых:
( )
( )
( ,
)
( )
( , )
(
, )

′
′′ =
′
′′
+
a
a
*
Сходство уравнений (7) и (13) легко прослеживается: роль падающего поля выполняет функция 
когерентности ΓG
n
a
( )( ,
)
′
′′
r r
; роль полного поля, регистрируемого на антенной решетке, выполняет 
разность Γ
Γ
ij
n
ij
( )
( )
−
0 . Для решения (13) надо предварительно вычислить ΓG
n
a
( )( ,
)
′
′′
r r
 с помощью (12), 
что на практике не очень удобно. Вместо этого 
вводятся обозначения Q
G
G
a
a
a
∫
r r
r
r
r r
r
r

G
n
n
d I
G
G
*
( ,
, )
( , )
(
, )
′
′′
≡
′
′′
r r
r
r r
r
r  
и  P
G
G
G
a
a
a
a
n
d I
G
G
r
r
r
r r
r r
d G
r
r
r
(
,
*
( )
*

+
′′
1
1
a
*
1
1)
( )
( )
( , )
( , )
*
( ,
, , )
( , )
( , )
(
, )
′
′′
≡
′
′′
r r
r r
r r
r r
r
r
1
1
1
 для произведений двух и трех функций Грина соответственно, и производится перегруппировка множителей:

∫
∫
′
+
ℜ
a
	
n
( )( ) ( , )
(
, )
r
r r
r
r
d G
d
r
r r
r
r
( , ) ( )

′′
+
+
′
a
I
G
G
1
1
1
1
a
*

∫
∫
ℜ
n
d
d
d I
Q
( )
( )
( )( )
( ,
, )
−
−
′
′′
′
′′
≈
Γ
Γ
a
∫
∫
∫
0
r
r
r
r
r r
r
*
Ω
ij
n
ij
S
S
j
i
	
(14)
r
r
r
,
)
(
)
d G
d G
r
r r
r
r
( , ) ( )
(

+
′
′
2
2
1
1
1
2
a
a
*

′
×
∫
ℜ
∫
ℜ


ε
1
1
a
n
r
r
r
r
r
r r
r r
d I
d
P
d
d
r
≈
′
′′
( )
*
( )
*
( )
( )
( ,
, , )
d I
G
G
n
r
r
r r
r r
( ) ( , )
( , ),
×
1
2
P
∫
∫
∫
′
′′
+
+
′′
ε
∫



1
a
*
′
r
r
S
S



ℜ
∫
r r r
, , )
.


( )
(
,
Ω
1
j
i
	
1
	

где ′
′′ ∈
r r
,
IRd . Первое слагаемое
	

G
n
n
d I
G
G
a
a
a
*
( )
( )
( ,
)
( )
( , )
(
, )
′
′′ =
′
′′
∫
r r
r
r
r r
r
r 	 (12)

Можно видеть, что уравнение (14) линейно относительно действительной ′
≡
ε
ε
( )
Re ( )
r
r  и  мнимой ′′
≡
ε
ε
( )
Im ( )
r
r  частей функции ε( )
r , которые 
рассматриваются как независимые переменные. 
После дискретизации исследуемой среды интегрирование заменяется матричным произведением, 
и уравнение (14) сводится к простой системе линейных уравнений. Функции P
a  в правой части 
отличаются порядком аргументов и комплексным 
сопряжением, т.е. соответствующие матрицы – 
гильбертово сопряженные.
Подводя итоги, можно сформулировать корреляционный итерационный метод определения 
волнового числа внутри неоднородности, что эквивалентно раздельному восстановлению скорости 
звука и амплитудного коэффициента поглощения, 
в виде следующего алгоритма.
1. При нулевом и при каждом n -м ракурсе облучения в эксперименте определяются значения 
матричных элементов Γij
n
( ) , и формируются разности Γ
Γ
ij
n
ij
( )
( )
−
0 . Первая оценка полагается равной 
 
известно, поскольку оно зависит от  известной 
мощности I n
( )( )
r  и функций Грина Ga( ,
),
r r′
 которые могут быть рассчитаны с помощью (6). Остальные слагаемые выражают нелинейную зависимость 
от  искомой функции ε( )
r . Такая нелинейность 
присутствует явно в последнем слагаемом, содержащем произведение двух таких функций, и неявно во втором и третьем слагаемых, поскольку в них 
входит неизвестная функция Грина исследуемой 
среды. В этом смысле уравнение (11) схоже с уравнением Липпмана–Швингера (5), которое также 
нелинейно относительно искомой функции ε( ).
r  
Для решения уравнения (11) предлагается применить итерационный подход, который нередко используется [17, 18, 22, 23] при решении уравнения 
типа (5). Он состоит в последовательном улучшении априорной оценки функции k( ),
r
 для чего 
в уравнениях (5) или (11) используется линеаризующее их борновское приближение (7); после этого 
ε
ε

1
2
2
( )
( )
( )
r
r
r
≡
=
−
a
0
a
k
k
.
АКУСТИЧЕСКИЙ ЖУРНАЛ
том 70
№ 2
2024

	
Корреляционный итерационный метод акустической томографии   
147
в) Уравнение (14) решается относительно функции ε( )
( )
( )
r
r
r
=
−
k
k
m
2
2
 . Поскольку выполнено 
тождество k
k
k
k
k
k
m
m
0
0
2
2
2
2
2
2
−
=
−









+
−









=
( )
( )
( )
( )
r
r
r
r



εm ( )
r
r
+ ε
2. Организуется бесконечный цикл, в рамках которого на каждой m -й итерации, начиная с m = 1, 
делаются следующие действия.
а) Имеющаяся в наличии оценка функции рассеивателя (r) 
ε


m
m
k
k
( )
( )
r
r
=
−
0
k
k
k
k
k
k
m
m
0
0
2
2
2
2
2
2
−
=
−









+
−









=
( )
( )
( )
( )
r
r
r
r



εm ( )
( )
r
r
+ ε
, следующая оценка функции ε
m +1( )
r  
полагается равной ε
ε
ε


m
m
+
=
+
1( )
( )
( )
r
r
r .
3. Чтобы определить, насколько оценка ε
m ( )
r  
соответствует входным данным, вводится безразмерная величина δ
ε
Γ( )

( ), равная 
2
2
, где km
 ( )
r  – это соответствующая ей оценка волнового числа, подставляется 
в уравнение (6) вместо величины εa( )
r , и рассчитывается функция Грина Ga( ,
)
r r′ .
б) С ее помощью определяются входящие в уравнение (14) произведения Qa( ,
, )
′
′′
r r
r  и P
a( ,
, , )
′
′′
r r
r r
1
.
	
2
0
( )
( )
n
d
d
d I
Q
r
r
r
r
r r
r
a
i j n
ij
n
ij
i j n
, ,
, ,
.
∑
∑
−
Γ
Γ 0 2
	
δ
ε
Γ
Ω
Γ
Γ
( )
( )
( ,
, )
( )
( )
( )
 ≡
−
−
′
′′
′
′′
∫
∫
∫
ij
n
ij
S
S
j
i
Цикл заканчивается, когда она становится ниже 
заданного порогового значения. При этом относительная невязка
метода конечных элементов. В этом случае удается 
уменьшить объем вычислений и сократить время 
расчетов, но требуется аккуратное задание условий 
на границе области наблюдения поля.
	
δ ε
ε
ε
ε
( )
( )
( )
( )


≡
−
2
2
d
d 	
(15)
ℜ
ℜ
∫
∫
r
r
r
r
r
3. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ 
РАБОТЫ КОРРЕЛЯЦИОННОГО 
ИТЕРАЦИОННОГО МЕТОДА
С целью проверки работоспособности метода 
моделировалась задача томографирования двумерной неоднородности ℜ квадратной формы, расположенной в начале координат и имеющей размеры 
характеризует близость полученной оценки ε
( )
r  
к истинной функции рассеивателя ε( )
r .
Уравнение (14), лежащее в  основе корреляционного итерационного метода, можно в  некоторых случаях упростить. Так, если контролируемые источники точечные и при каждом n -м 
ракурсе облучения расположены в точке r( )
n , то 
5
5
0
0
λ
λ
×
, где 

0
0
2
≡
k  – длина волны в окружающей фоновой среде. Скорость звука c x y
( , ) и амплитудный коэффициент поглощения α( , )
x y  внутри неоднородности задавались в виде 
I
I
n
n
n
( )
( )
( )
( )
r
r
r
=
−
(
)
δ
, и приемники тоже точеч

ные и имеют координаты rj , то с учетом этого 
вместо (14) можно записать
	
c x y
c
s
x
y
( , )
,
=
+

1



 










0
0





0
0





( )
( )
0
r r r
,
,
и
a


ij
n
ij
n
i
j
n
( )
( )
−
−
(
) ≈
(16)
	
		
( )
d
P
( )
,
, ,
	
	
a
i
j
n
I
Q
≈
(
) +
r
r
r r r r

*
( )
( )
, , ,
.
r
r r r r
1 P
j
i
n
a
*(
)






ℜ
∫
	




( , )
,
,
x y
k s
y
x
=











0 0
0
0
соответственно, где c
k
0
0
= ω
 – скорость звука в фоновой среде. Безразмерные параметры α0  
и β0  варьируются, обеспечивая изменение контраста неоднородности. Функция s x y
( , )  задает 
пространственное распределение параметров внутри неоднородности. Для удобства данные об исследованных в настоящей статье неоднородностях 
и значения описанных характеристик объединены 
в таблицу. NR = 20  точечных приемников располагались в  точках с  полярными координатаСледует обратить внимание, что итерационные 
алгоритмы решения обратной задачи рассеяния, 
как правило, предполагают численное решение 
прямой задачи рассеяния на каждой итерации [17, 
18]. В рамках описанного алгоритма она возникает 
на втором шаге, когда требуется определить G a( ,
).
r r′  
Погрешности, которые привносятся при этом, оказывают существенное влияние на точность итоговой оценки. В частности, значительный вклад оказывают процессы многократного рассеяния внутри 
каждого элемента разрешения. Эти вопросы были 
затронуты в [24]. Альтернативной возможностью является использование для решения прямой задачи 
ми rj
R
j
N
=
−
(
)
100
2
1
0
λ
π
;
(
)
;  NF = 20  точечных 
АКУСТИЧЕСКИЙ ЖУРНАЛ
том 70
№ 2
2024

ДМИТРИЕВ

Таблица. Параметры неоднородностей
M
№
s x y
( , )
α0
β0
f ( ,
)
θ θ′
∆φ
π
2




1
exp
.
2
(
)
1
2
0 5
1
2
0.1
0.3
13
3π
0.9
1.3
2
0 25
2
2
0 25
.
.
exp
.
−
+
−
−
(
)
×
×







−
−
+
+
−
x
iy
x
iy
0 5
1

























x
y
2
2
0
0.4
21
3π
1.2
1.9





2
2
0 75
2
exp
.
−
+
×













x
y
x
2
2
2
2
0 75
2π
0
0.4
12
3π
1.2
1.5
3
exp
.
cos(
)
−
+
×






контролируемых источников – в точках с полярмаксимума ∆φ
ζ
θ
= max ( , )
a
 преобразования Радона 
ными координатами r( )
;
(
. )
n
F
n
N
=
−
(
)
100
2
0 5
0


. 
Шаг дискретизации был выбран равным λ0 8 . При 



( , )
( )
a
k
c
=
−










Rad
0
r
 от подынтегрального вы- 
ражения, которое можно рассчитать быстро.
Третья характеристика равна максимальному 
G
G
 
модулю отношения M ≡
′ −
′
′
′
max
( ,
)
( ,
)
G
a
( ,
)
,
r r
r r
r r
r r
a
моделировании здесь и далее априорная информация не использовалась, т.е. полагалось k
k
a( )
r =
0.
Чтобы численно определить сложность восстановления каждой конкретной моделируемой неоднородности, вводится несколько ее характеристик, 
по-разному связанных со свойствами среды. Первая из них – норма амплитуды рассеяния f ( ,
)
θ θ′  , 
2
2
. Амплитуда 
равная f
f
d d
( ,
)
( ,
)
 
 
 


′
≡
′
′
∫
∫
2
0
0
рассеянного поля G
G
( ,
)
( ,
)
r r
r r
′ −
′
a
 к падающему 
полю Ga( ,
)
r r′ . В зависимости от того, является 
ли она много меньшей единицы, меньшей единицы или превосходит единицу, неоднородность 
относят к  классу слабых, средних или сильных 
рассеивателей соответственно [17]. Сильные рассеиватели наиболее сложны для восстановления, 
и итерационные процедуры, в том числе предлагаемый в настоящей статье метод, не гарантируют успеха. Более того, в двумерном случае, когда 
объем исходных данных конечен, в классе сильных 
рассеивателей решение задачи не единственно [18].
На рис.  1а представлено пространственное 
распределение относительной скорости звука c x y c
( , )
0  для неоднородности № 1 с α0
0 1
=
. ;  
рассеяния определяется в пространствах разной 
размерности d  аналогично тому, как это сделано 
в ряде работ [18, 27]; θ  и ′
θ  – полярные углы, под 
которыми распространяются падающая и рассеянная волны. При численном расчете f ( ,
)
θ θ′
 используется уравнение (6).
Вторая характеристика определяется как максимальный по  различным траекториям ℓ дополнительный набег фазы, который получается 
при распространении сигнала через неоднородность по сравнению с набегом фазы в фоновой 
β0
0 3
= . . Она описывается функцией 
2


ℓ
k
c
d
0
r
. Черта над max  
среде: ∆φ
ω
=
−











∫
max
( )
ℓ
ℓ
( , )
exp
.
s x y
x
iy
≡
−
+
−
−
(
)
1
2
0 5
1
2
0 25
.
×







−








	
2


x
iy
0 5
1
+
+
−
(
)
exp
−
−
2
.
1
2
.
,
2
0 25
×















означает, что эта операция в данном случае нестандартная. Она заключается в том, что вначале определяется наибольшее по абсолютной величине значение аргумента, а затем результату 
приписывается тот же знак, что и у аргумента. 
Такое сохранение знака позволяет учесть, что 
набег фазы в среде может быть как положительным, так и отрицательным. Если для оценки полагать каждую траекторию ℓ прямой линией, 
то  такая процедура эквивалентна вычислению 
представляющей собой сумму двух “колец” и имеющей достаточно сложный вид. Форма пространственного распределения коэффициента 
поглощения такая же: она получается из формы 
АКУСТИЧЕСКИЙ ЖУРНАЛ
том 70
№ 2
2024

	
Корреляционный итерационный метод акустической томографии   
149
(а)
(б)
(в)
100
1.2
1.0
1.4
1
1
10−1
2
0.8
1.2
2
0.6
c/c0
B
ε/k0
2
1.0
3
10−2
0.4
δ(ε
m), δΓ(ε
m)
0.2
A
3
4
10−3
0
y/λ0
0.8
3
3
1
1
−1
−1
−3
−3
−1.5
−2.5
1.5
2.5
1
0.5
−0.5
−0.2
4
10
100
500
x/λ0
x/λ0
m
Рис. 1. (а) – Пространственное распределение относительной скорости звука c x y c
( , )
0  для неоднородности № 1. 
(б) – Действительные (линии 1 и 2) и мнимые (линии 3 и 4) части нормированных на k0
2  оценки ε
500( , )
x y  (линии 1 и 3) и искомой функции рассеивателя ε( , )
x y  (линии 2 и 4) вдоль отрезка AB.  (в) – Зависимости величин  
δ( m) (линии 1 и 2) и δГ( m) (линии 3 и 4) от номера итерации m  при точных входных данных (линии 2 и 4) и при 
наличии помех (линии 1 и 3).
сравнительно небольшие помехи, с одной стороны, замедляют скорость сходимости метода, 
а с другой стороны, приводит к тому, что результатом восстановления является неоднородность, 
существенно отличающаяся от  искомой. Данное обстоятельство является типичным при решении обратных задач в присутствии сильного 
рассеивателя.
распределения скорости звука симметрией относительно прямой y
x
=
. Несмотря на небольшой 
волновой размер, она относится к классу сильных рассеивателей. Результат ее восстановления 
с помощью корреляционного итерационного метода приведен на рис. 1б. Здесь представлены выполненные вдоль отрезка AB  и нормированные 
на k0
4. ВОЗМОЖНОСТИ УСКОРЕНИЯ 
ИТЕРАЦИОННОЙ ПРОЦЕДУРЫ
2  сечения полученной после 500 итераций 
оценки ε
500( , )
x y  и искомой функции рассеивателя ε( , ).
x y
 Хорошее визуальное согласование кривых в сочетании с низким достигнутым значением 
При численном моделировании область ℜ, 
занимаемая рассматриваемой неоднородностью, 
разбивается на отдельные элементы с размером, 
равным выбранному шагу дискретизации и характеризующиеся функцией рассеивателя εr , где r  
в данном случае – индекс такого элемента. Тогда (14) сводится к системе линейных уравнений 
стандартного вида Anij;r  r = Bnij , где n
NF
= 1;
 и 
i j
NR
,
;
= 1
. Эта система решается с помощью ре−
A
A
E
A
B
H
H
1
, 
гуляризации Тихонова: ε
κ
 =
+
(
)
δ ε
(
)
.
500
0 03
≈
 говорит о работоспособности метода 
и его применимости для восстановления подобных 
неоднородностей. На рис. 1в линиями 2 и 4 изображены зависимости от номера итерации m  величин δ ε
(
)
m  и δ
ε
Γ(
),
m
 соответственно. Видно, что 
обе они достаточно быстро стремятся к нулю.
Для проверки устойчивости метода в его входные данные, т.е. матрицы Γij
n
( )  вносились помехи. 
С этой целью ко всем принятым сигналам U j
n
( )  добавлялись независимые комплексные случайные 
величины, действительные и мнимые части которых были распределены нормально с дисперсией, 
( )
равной 
χ2
2
2N N
U
F
R
j
n
,
∑
; χ =
−
10 2. Зависимости 
j n
δ ε
(
)
m  и δ
ε
Γ(
)
m  для данных с помехами изображены на рис. 1в линиями 1 и 3. Видно, что обе они 
лежат выше исходных зависимостей; при этом 
зависимость δ ε
(
)
m  не стремится к нулю с увеличением номера итерации, а зависимость δ
ε
Γ(
)
m  – 
стремится. Все это говорит о том, что введенные 
где верхний индекс “H” означает гильбертово сопряжение; E  – единичная матрица соответствующего размера, и κ  – коэффициент регуляризации. На практике выбор значения κ  может быть 
затруднен. С одной стороны, если этот коэффициент мал, то решение теряет устойчивость, и относительная невязка δ( m) при восстановлении неоднородности может быстро возрастать. С другой 
стороны, большое значение κ  приводит к уменьшению производимой на каждой итерации добавки к восстанавливаемой функции рассеивателя, 
АКУСТИЧЕСКИЙ ЖУРНАЛ
том 70
№ 2
2024