Известия РАН. Физика атмосферы и океана, 2020, T. 56, № 1, стр. 66-75

ВВЕДЕНИЕ

Метод фазовых экранов давно и широко применяется для численного моделирования распространения волн различной природы в неоднородных средах. В качестве примеров укажем задачи моделирования распространения оптического (лазерного) излучения в турбулентной атмосфере [1–5], задачи моделирования распространения дециметровых волн при радиозатменном зондировании атмосферы [6–8]. Этот метод в англоязычной литературе называется split-step, метод пошагового расщепления. Соответственно в русскоязычной литературе используется название “метод расщепления по физическим параметрам”. Это название отражает тот факт, что неоднородная среда в этом методе представляется как последовательность тонких слоев, а пропагатор, описывающий распространение волны через каждый слой, приближенно записывается как композиция бесконечно тонкого слоя, формирующего фазовые искажения волны, и вакуумного пропагатора, описывающего дифракцию.

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

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

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

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

ВЫВОД ОСНОВНЫХ СООТНОШЕНИЙ. ВАКУУМНЫЙ ПРОПАГАТОР

В качестве примера мы будем использовать параметры пучка, заданные для проекта DELICAT (DEmonstration of LIdar based Clear Air Turbulence detection, Демонстрация обнаружения турбулентности ясного неба при помощи лидара) [11, 12]: длина волны λ = 354.84 нм, расходимость $2\alpha = 0.3$ мрад. Сам излагаемый метод может применяться и при других параметрах пучка, в рамках условий применимости используемого приближения.

Будем рассматривать декартовы координаты $\left( {x,y,z} \right)$ и сферические координаты $\left( {r,\theta ,\varphi } \right)$, связанные между собой следующим образом:

Рассмотрим волновое поле в пространстве $u\left( {x,y,z} \right)$ с наложенным условием излучения, выделяющим волны, распространяющиеся в направлении оси $x$. Пусть для волнового поля задано граничное условие ${{u}_{0}}\left( {y,z} \right)$ в плоскости $x = 0.$ Будем пользоваться преобразованием Фурье в следующей нормировке:

где $k = \frac{{2\pi }}{\lambda }$ – волновое число. Тогда волновое поле в вакууме записывается следующим образом: где $S\left( {r,\theta ,\varphi ;\eta ,\zeta } \right)$ – фазовая функция, имеющая следующий вид:

Максимальные абсолютные значения $\eta $, $\zeta $, $\theta $ и $\varphi $ оцениваются величиной половинной расходимости пучка α = 0.15 × 10⁻³. Расстояние распространения $r$ оценивается величиной 30 км. В дальнейшем мы будем пользоваться малоугловым разложением фазовой функции, поэтому оценим, члены каких степеней следует удерживать:

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

Пользуясь этим приближением, запишем поле в пространстве в виде преобразования Фурье:

Для численного нахождения поля в неоднородной среде будем применять метод фазовых экранов. В рамках этого метода каждый экран будем считать участком сферы с некоторым (текущим) радиусом $r$. Будем вычислять поле на следующем экране с радиусом $r + \Delta r$, пользуясь вакуумным пропагатором и умножая полученное поле на множитель, приближенно учитывающий среду. Для вакуумного пропагатора несложно получить явное выражение. Для этого, пользуясь полем на исходном экране, будем вычислять эффективный пространственный спектр $\widetilde {{{u}_{0}}}\left( {r,\eta ,\zeta } \right)$:

Эффективный спектр равен тому пространственному спектру поля в плоскости $x = 0,$ при котором на данном фазовом экране в отсутствие среды формировалось бы поле $u\left( {r,\theta ,\varphi } \right).$ Далее, пользуясь формулой (7) с подстановкой $r + \Delta r$ и $\widetilde {{{u}_{0}}}\left( {r,\eta ,\zeta } \right),$ и учитывая влияние среды, мы получаем поле на следующем экране:

где $\Psi \left( {r,\Delta r,\theta ,\varphi } \right)$ – добавка оптического пути вдоль прямых лучей в слое между $r$ и $r + \Delta r$ за счет среды.

Перепишем формулу (8) в следующей форме:

Тогда поле на следующем фазовом экране можно записать аналогично:

Подставляя сюда формулу (10), получаем:

Выделим здесь интегрирование по промежуточным переменным $\eta ,\zeta $:

Выполняя интегрирование по $\eta $ и $\zeta $, получаем:

Вычисляя этот интеграл для сферической волны $u\left( {r,\theta {\kern 1pt} ',\varphi {\kern 1pt} '} \right) = \operatorname{const} ,$ получим:

Пользуясь представлением:

формулу (14) можно переписать в следующем виде: где $P\left( {r,\Delta r;\zeta ,\eta } \right)$ – вакуумный пропагатор в малоугловом приближении:

Это представление соответствует стандартному методу фазовых экранов (split-step). Пересчет поля с экрана на экран выполняется в три шага. На первом шаге вычисляется пространственный спектр поля в исходном фазовом экране. На втором шаге спектр умножается на вакуумный пропагатор и при помощи обратного преобразования Фурье вычисляется поле на следующем фазовом экране без учета среды. На третьем шаге это поле умножается на фазовый множитель, учитывающий фазовый набег в среде между экранами.

Множитель $\exp \left( {ik\Delta r} \right)$ можно опустить, поскольку он дает постоянную добавку к фазе в каждом экране. В этом случае вакуумный пропагатор можно записать в следующем виде:

Перепишем формулу (17) в операторном виде:

где $\hat {P}\left( {r,\Delta r} \right)$ – оператор описывающий распространение поля с экрана на экран. Рассмотрим 3 последовательных фазовых экрана ${{r}_{1}}$, ${{r}_{2}}$ и ${{r}_{3}}$. Рассмотрим композицию операторов $\hat {P}\left( {{{r}_{2}},\Delta {{r}_{{23}}}} \right) \circ \hat {P}\left( {{{r}_{1}},\Delta {{r}_{{12}}}} \right),$ где $\Delta {{r}_{{23}}} = {{r}_{3}} - {{r}_{2}}$ и $\Delta {{r}_{{12}}} = {{r}_{2}} - {{r}_{1}}.$ При этом будем считать, что среда является вакуумом, т.е. $\Psi \left( {r,\Delta r,\theta ,\varphi } \right) = 0.$

Выделим здесь интегрирование по $d\varphi {\kern 1pt} 'd\theta {\kern 1pt} '$:

Вычисляя внутренний интеграл по $d\varphi {\kern 1pt} 'd\theta {\kern 1pt} '$, получаем:

Это позволяет переписать интеграл (22) следующим образом:

Выполняя интегрирование по $d\zeta d\eta ,$ и учитывая стандартное свойство $\delta {\text{ - }}$функции, получаем:

Рассмотрим выражение $P\left( {{{r}_{1}},\Delta {{r}_{{12}}};\zeta {\kern 1pt} ',\eta {\kern 1pt} '} \right)$ × $ \times \,\,P\left( {{{r}_{2}},\Delta {{r}_{{23}}};\frac{{{{r}_{1}}}}{{{{r}_{2}}}}\zeta {\kern 1pt} ',\frac{{{{r}_{1}}}}{{{{r}_{2}}}}\eta {\kern 1pt} '} \right)$:

где $\Delta {{r}_{{13}}} = {{r}_{3}} - {{r}_{1}}$ – расстояние между первым и третьим экранами. Учитывая это тождество, из (25) получаем:

Этим доказано естественное групповое свойство вакуумных операторов $\hat {P}\left( {r,\Delta r} \right)$:

Из точности группового свойства следует, что данный пропагатор является точным решением приближенной формы параболического уравнения, в котором в качестве поперечных координат используются ${y \mathord{\left/ {\vphantom {y r}} \right. \kern-0em} r}$ и ${z \mathord{\left/ {\vphantom {z r}} \right. \kern-0em} r},$ а также учтено регулярное убывание амплитуды сферической волны как ${1 \mathord{\left/ {\vphantom {1 r}} \right. \kern-0em} r}.$ Альтернативный подход к учету расходимости волны возможен с использованием линзовых координат [13–15]. Тем не менее в указанной литературе линзовые координаты применяются для описания пучков в однородной нелинейной диспергирующей среде. Возможность применения этого преобразования к произвольной неоднородной среде требует дополнительного исследования.

Метод сферических фазовых экранов может также применяться для исследования сфокусированных пучков, в том числе в условиях самофокусировки. Пропагатор (18) при $\Delta r < 0$ описывает не расходящийся, а сходящийся пучок.

ФАЗОВЫЙ НАБЕГ В ТОНКОМ СФЕРИЧЕСКОМ СЛОЕ

Каждый слой случайно-неоднородной среды толщиной $\Delta r$ между фазовыми экранами $r$ и $r + \Delta r$ описывается функцией $\Psi \left( {r,\Delta r,\theta ,\varphi } \right),$ которая понимается как реализация 2-мерного случайного поля как функция от $\theta ,\varphi $ при заданных $r$ и $\Delta r$. Рассмотрим соответствующую реализацию 3-мерного поля индекса рефракции $N\left( {x,y,z} \right)$. Тогда можно записать:

Будем пользоваться для вычисления этой функции малоугловым приближением. Поскольку $r{{\alpha }^{2}}$ = 0.675 мм, то величинами 2-го порядка можно пренебрегать. Таким образом, можно записать следующее приближенное выражение:

Мы будем пренебрегать регулярной частью индекса рефракции, считая его константой, дающей в каждом экране некоторый известный несущественный постоянный фазовый сдвиг. Поле $N\left( {x,y,z} \right)$ будет считаться трехмерным статистически однородным и изотропным случайным полем со спектральной плотностью ${{\Phi }_{N}}\left( {\mathbf{\kappa }} \right) = {{\Phi }_{N}}\left( \kappa \right),$ где ${\mathbf{\kappa }} = \left( {{{\kappa }_{x}},{{\kappa }_{y}},{{\kappa }_{z}}} \right)$ – трехмерный вектор пространственных частот, $\kappa = \left| {\mathbf{\kappa }} \right|$. Для корреляционной функции имеет место стандартное соотношение:

Рассмотрим корреляции фазового пути:

Выделим здесь интегрирование по :

Этот интеграл равен:

где мы ввели обозначения и Эта величина не является функцией от $\Delta \varphi $ и $\Delta \theta $. Однако в малоугловом приближении поправки, включающие ${{\kappa }_{y}}\varphi {\kern 1pt} '\,\, + {{\kappa }_{z}}\theta {\kern 1pt} '$ и либо не в виде разностей, либо с множителем $\Delta r$, как в числителе, так и в знаменателе можно приближенно исключить. При этом $\Delta r$ считается малым по сравнению с $r$. Тогда выражение (34) можно переписать в следующей приближенной форме:

Пользуясь функцией $\operatorname{sinc} x = {{\sin x} \mathord{\left/ {\vphantom {{\sin x} x}} \right. \kern-0em} x}$, это выражение можно переписать в следующей форме:

Это позволяет преобразовать интеграл (32) следующим образом:

Получим выражение для спектральной плотности фазового пути ${{\Phi }_{\Psi }}\left( \mu \right)$, где ${\mathbf{\mu }} = \left( {{{\mu }_{\varphi }},{{\mu }_{\theta }}} \right)$ – вектор пространственных частот. Имеет место следующее соотношение:

Если положить ${{\mu }_{\varphi }} = r{{\kappa }_{y}}$ и ${{\mu }_{\theta }} = r{{\kappa }_{z}},$ то мы придем к следующему выражению для спектральной плотности фазового пути:

Эта формула записана в приближении тонкого слоя среды, поэтому она не учитывает то, что в силу уширения пучка в угловые частоты ${\mathbf{\mu }}$ дают вклад не только поперечные частоты ${{\mathbf{\mu }} \mathord{\left/ {\vphantom {{\mathbf{\mu }} r}} \right. \kern-0em} r}.$ Для слоя среды толщиной порядка внешнего масштаба неоднородностей и более можно записать по аналогии с соответствующей формулой для плоского слоя [1, 2]:

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ

Для турбулентных флуктуаций индекса рефракции мы примем в качестве модели колмогоровско-кармановский спектр:

где $C_{n}^{2}$ – структурная константа, ${{\kappa }_{0}} = {{2\pi } \mathord{\left/ {\vphantom {{2\pi } {{{L}_{0}}}}} \right. \kern-0em} {{{L}_{0}}}},$ ${{L}_{0}}$ – внешний масштаб, ${{\kappa }_{m}} = {{5.92} \mathord{\left/ {\vphantom {{5.92} {{{\lambda }_{0}}}}} \right. \kern-0em} {{{\lambda }_{0}}}},$ и ${{\lambda }_{0}}$ – внутренний масштаб турбулентности.

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

Для генерации реализаций случайных некоррелированных фазовых экранов мы пользуемся дискретной формой соотношения (38) [16]:

где ${{\Psi }_{{j'l'}}} = \Psi \left( {{{\varphi }_{{j'}}},{{\theta }_{{l'}}}} \right)$ и ${{\tilde {\Psi }}_{{jl}}}$ является дискретным преобразованием Фурье от ${{\Psi }_{{j'l'}}}$:

Таким образом, ${{\tilde {\Psi }}_{{jl}}}$ являются некоррелированными случайными величинами со случайными фазами и среднеквадратичными значениями:

Для генерации случайной реализации ${{\tilde {\Psi }}_{{jl}}}$ рассмотрим независимые случайные величины $\xi _{{jl}}^{{(1)}}$ и $\xi _{{jl}}^{{(2)}}$, равномерно распределенные в интервале [0, 1]. Тогда реализация ${{\tilde {\Psi }}_{{jl}}}$ дается формулой:

поскольку комплексная случайная величина $\sqrt { - \ln \xi _{{jl}}^{{(1)}}} \exp \left( {2\pi i\xi _{{jl}}^{{(2)}}} \right)$ имеет гауссово распределение со среднеквадратичным значением, равным 1. С учетом (40) можно записать: где

Реализация оптической толщины экрана ${{\Psi }_{{jl}}}$ равна обратному дискретному преобразованию Фурье от ${{\tilde {\Psi }}_{{jl}}}$ [5]:

Для учета того, что один и тот же шаг дискретизации по угловым координатам при разных $r$ соответствует разным пространственным масштабам, мы используем адаптивную дискретизацию. На каждом фазовом экране выполняется оценка необходимого шага дискретизации в угловых координатах, и если текущий шаг дискретизации превышает эту оценку, то выполняется удвоение разрешения. Для этой цели используется интерполяция Котельникова: в пространстве Фурье-образов сетка частот дополняется так, чтобы при том же частотном разрешении максимальная частота стала вдвое выше. В добавленных узлах сетки Фурье-образ поля задается равным 0. После обратного преобразования Фурье мы получаем поле, интерполированное на пространственную сетку с вдвое меньшим шагом дискретизации.

Рассмотрим вопрос численной эффективности метода. Как для плоских, так и для сферических экранов шаг дискретизации должен быть не более характерного масштаба среды (порядка или больше внутреннего масштаба для турбулентности). Однако для плоских фазовых экранов необходимо также учесть, что на краях имеется регулярный набег фазы порядка $\varphi = {{kr{{\alpha }^{2}}} \mathord{\left/ {\vphantom {{kr{{\alpha }^{2}}} 2}} \right. \kern-0em} 2},$ составляющий в нашем примере около 5000 радиан. Необходимый шаг дискретизации равен $r\Delta \alpha = \frac{{r\Delta \varphi }}{{{{d\varphi } \mathord{\left/ {\vphantom {{d\varphi } {d\alpha }}} \right. \kern-0em} {d\alpha }}}}$ = $\frac{{\lambda \Delta \varphi }}{{2\pi \alpha }},$ где $\Delta \varphi $ – допустимая вариации фазы между дискретными выборками поля. Если принять $\Delta \varphi = \frac{\pi }{4},$ то $r\Delta \alpha = \frac{\lambda }{{8\alpha }} \approx 0.2$ мм. При указанной в статье расходимости характерный размер луча достигает 5 м, что потребует размерности сетки порядка $N = 6 \times {{10}^{8}}$ точек, что является очень большой величиной, требующей не только большого времени, но и оперативной па-мяти, не всегда доступной на современных процессорах. Фактически же для счета достаточно шага около 3 мм, что понижает размерность сетки более, чем в 100 раз, а скорость вычисления быстрого преобразования Фурье, пропорциональная $N\lg N,$ увеличивается в еще большей мере.

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

РЕЗУЛЬТАТЫ

На рис. 1–6 показаны примеры численного расчета распределений амплитуды гауссова лазерного пучка с длиной волны $\lambda = 354.84$ нм с расходимостью $2\alpha = 0.3 \times {{10}^{{ - 3}}}$ рад на различных точках трассы длиной 24 км для $2.5 \times {{10}^{{ - 15}}}\,\,{{{\text{м}}}^{{ - 2/3}}}$. Рисунки демонстрируют формирование режима спеклов. Шаг между фазовыми экранами составлял 100 м.

ЗАКЛЮЧЕНИЕ

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