CC BY
20СТАТИСТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПЕРЕНОСА ОПТИЧЕСКИХ СИГНАЛОВ ЧЕРЕЗ ВЗВОЛНОВАННУЮ ПОВЕРХНОСТЬ
О, С, Ухинова1, Б, А. Каргин1,2
1 Институт, вычислительной математики и математической геофизики СО РАН, 630090, Новосибирск 2 Новосибирский государственный университет, 630090, Новосибирск
УДК 517.968
Б01: 10.24411/9999-016А-2019-10085
В работе рассмотрены две модели оценки оптической передаточной функции системы океан-атмосфера, предложенные в [5]. После определенных уточнений были проведены численные эксперименты с использованием второй из предложенных моделей (для первой модели численные эксперименты были проведены также в работе [5]). Выло проведено сравнение предложенных методов оценки, которое показало преимущество последней модели над предыдущей.
Ключевые слова: метод Монте-Карло, система океан-атмосфера, оптическая передаточная функция.
Введение
Строгое решение подавляющего числа прямых задач оптики системы океан-атмосфера возможно лишь в статистической формулировке. Это обусловлено тем, что поле оптического излучения системы океан-атмосфера в значительной мере формируется под воздействием взволнованной морской поверхности, отражающей и преломляющей приходящее электромагнитное излучение. Характерная особенность этого воздействие обусловлена случайным характером ветрового волнения морской поверхности. В работе для описания взволнованной поверхности моря принята так называемая плоская "фацетная" модель, которая впервые была применена к расчету переноса солнечной радиации в системе атмосфера-океан в [1,2]. Также в работе рассмотрены алгоритмы метода Монте-Карло для моделирования переноса излучения через взволнованную морскую поверхность в применении его для вычисления оптической передаточной функции (ОПФ) системы океан-атмосфера. Представлены некоторые количественные оценки влияния ряда факторов на ОПФ системы для одной из моделей, а также проведено сравнение двух рассмотренных моделей в применении в вычислению оценок ОПФ для различных скоростей ветра.
1 Постановка задачи
Постановка задачи полностью совпадает с задачей в работе [5]. Верхняя граница г = Н плоской атмосферы освещается параллельным потоком фотонов 106(и> — шо) в напралении ш0. Взаимодействие света с веществом атмосферы и океана определяется заданием коэффициентов ослабления а (г) и рассеяния ае(г), индикатрисы рассеяния д(г, р), нормированной условием д(г, ^)с!^ = 1, где ^ — косинус угла рассеяния. Введем обозначения: = (х, у) — проекция радиуса-вектора г = (х, у, х) па горизонтальную плоскость; ш = (а, Ь, с) — едшвчлый вектор направления с направляющими косинусами а, Ь,с; П — множество всех направлений с с € [—1,1], — единичная полусфера с с € [0,1], — единична полусфера с с € [—1, 0]. При попадании света на плоскость г = 0 он испытывает отражение, которое характеризуется коэффициентом яркости ^(г^,равным отношению яркости света, идущего из точки (г^, 0) в направлении ш к освещенности единичной площадки в точке (г±, 0), ориентированной перпендикулярно к ш'. В плоскости г = к
Работа выполнена в рамках государственного задания ИВМиМГ СО РАН (проект 0315-2019-0002).
ISBN 978-5-901548-42-4
расположена граница раздела двух сред, при взаимодействии с которой свет испытывает преломление и отражение. Эта граница представляет собой случайную поверхность, составленную из набора элементарных площадок, центры которых лежат в плоскости г = Н, а нормали к площадкам I — случайные единичные векторы с плотностью распределения вероятности р(1), I = (1Х, 1у, 1г) € (для такой модели взволнованной морской поверхности в соответствующей литературе утверждено название "фацотная модель").
Рис. 1: Схема просесса переноса излучения через границу раздела двух сред 1 детерминированная плоская граница. 2 случайная граница.
Обозначим через Ф(г^, я, ш) полную яркость светового поля в точке (г±, г) в направлении ш. Световой режим в плоскопараллелыгой среде, освещенной параллельным потоком фотонов, описывается известным интегро-дифференциальным уравнением переноса с соответствующими условиями (см. [3]. стр. 5). Для описания светового поля в системе океан-атмосфера к этим краевым условиям необходимо добавить выражение, определяющее взаимодействие света со случайной границей раздела двух сред:
Ф(г^,к,ш) = / ¿1 р(1)К(ш',1)Ф(г^,к,ш')5[ш' - шотр(ш,1)]вш' + (1)
шеп П- п
/ ¿1 р(1)[1 - К(ш',1)]Ф(г^,к,ш')6[ш' - (ш,1)]Ли' ■! П- -П1
Здесь К(ш, I) — френелевский коэффициент отражения,
шотр(ш, I) = ш - 2(ш, 1)1,
(ш, I) = и(ш + 71), 7 = *1дп[(и>, 1)]^1 ^ - 1 + (ш,1)2 - (ш, I),
п, (ш',1) < 0; 1/п, (и',1) > 0;
п — коэффициент преломления воды относительно воздуха,
В случае, когда поверхность раздела двух сред является статической, р(1) = 6(1 - к) и
Ф(г±, Н,ш)
К(штр,к)Ф(г^, р) + и2[1 - К(шпр,к)]Ф(г±, к,шпр),
шеп
где к = (0,0,1).
Яркость Ф(г, является суперпозицией световой дымки системы океан-атмосфера Фо(-г, ш) и яркости отраженного плоскостью г = 0 излучения Ф(г^,г,ш). На основе теории линейных систем [4] введем ОПФ системы океан-атмосфера Т(р,х,ш) такую, что
Ф(, к,ш)
= Т(Р, %, Ш)
шеп
Ф(р, 0, ш)
шеп
(2)
где символом Л обозначен Фурье-образ, р = (рх,ру)- вектор пространственных частот. Функция Т(р,х,ш) определяется только оптическими и геометрическими свойствами системы океан-атмосфера и не зависит от условий освещения и пространственного спектра яркости Ф(г±, 0, и)шеП_. Таким образом, функция Т(■)
позволяет для произвольного распределения яркости Ф(■)
вблизи плоскости г = 0 с помощью прямого и
обратного преобразования Фурье определить пространственно-угловое распределение яркости Ф(■) на любой высоте над поверхностью океана. И наоборот, по измерениям поля яркости системы с помощью функции Т( ■ ) = 0
2 Численный алгоритм
В данной работе рассматривается дополнительная сопряженная задача, в которой источником света является светящаяся плоскость г = 0 Яркость последней обозначим через 7(€ О-. Световое поле в такой системе описывается упоминавшимся интегро-дифференциальным уравнением переноса с краевыми условиями
Ф (г±,Н,ш)
= 0, Ф(г±, 0,ш)
шеп+
= 7 (г±,ш)
шеП-
шеП-
вместе с соотношением (1).
Т( ■ )
соб оценки методом Монте-Карло пограничной кривой, описывающей световое поле в плоском слое мутной среды от источника с прямолинейным разрывом яркости. Алгоритм построен на сочетании методов "коррелированной выборки" и "сопряженных траекторий". Не вдаваясь в общую схему построения траекторий фотонов и вычисления соответствующих оценок (см., например [3]) отметим лишь специфичные для рассматриваемой задачи особенности.
Выберем для определенности систему координат хух так, чтобы граница разрыва яркости 7(■) была ориентирована перпендикулярно осп ОХ, т.е.
к ) V ) /0 ,ж< 0;
I 7 (ш) ,х > 0.
В этом случае Ф(г±, г*,ш*) = Ф(х, г*, ш*). Зададим на оси ОХ сетку узлов {х^, г = 1,..., Ж}, в которых вычисляются значения пограничной кривой Ф^ = Ф(х¿, х*, ш*) для направления ш*. Согласно [6] траектории фотонов строятся из произвольной точки Р = (хр, ур, г *) плоскости г = г* в направле нии — ш* и закапчиваются либо вылетом в полупространство г > либо попаданием фотона на плоскость г = 0. В случае пересечения траектории фотона с плоскостью г = к в направлении из' вычисляется коэффициент К(ш', в), где
из распределения с плотностью р(/). С вероятностью Р(ш, в) происходит зеркальное отражение и фотон приобретает направление ш = из' — 2(из', в)в. С вероятностью 1 — К(из', в) происходит преломление в направлении ш = из'/у — /в. Коэффициент К(ш, I) вычисляется по известной формуле Френеля (см., например, [7]).
Использование коррелированных траекторий и локальной оценки "по направлению" [3] позволяет по одним и тем же траекториям оцепить функцию Ф(х, х*, из*) для всего набора точек {х^, г = 1,..., N} одновременно. Случайные оценки ^ вкладов в функционалы Ф¿(-) от столкновения в точке (х, у, г) с направлениями ш' ти ш до и после столкновения с учетом (1) приобретают вид:
а(г)
IV ехр(-т)[1 - Я(и,в)]п21 (р, -т) ,г е (к, Н];
а(г
IV ехр(-т(р, -ш) ,г е (0,к].
г - к к
Хг + х - хр + (а, ш)-.-г- + (а, т)
XI + х - хр + (а, си)
(-ш к) (-ш,к)
(-т, к)
,г е (к, Н];
, г е (0,к],
где т = ш/п - ^ = [(ш, а) + \/п2 - 1]/п, а = (1, 0, 0),
/ а[г +^ш,к)]Л + / а[г + ^(ш,к) +^т,к)]Л, 'о ■! о
т = ^ где Ь1 = (г - к)/(-ш, к), Ь 2 = к/(-т,к) ,г е (к,Н];
Г '1
/ а[х+г(ш,к)](И, где 11 = г/(-ш, к) ,г е (0,к],
о
при (и, к) < ^^и (и, к) > 0 значения ^ = 0.
В (3) ^ — "вес" фотона до столкновения (см. [6]). В точке (х,у, г) пересечения траектории фотона с плоскостью г = к в направлении ш имеем:
(3)
6 =
Ш[1 - Е(ш, в)] ехр(-т(р, -т)п2 , т определяется так же как в (В),ш е в)] ехр(-г (р, -т) ,т = ш - 2(ш, в) в, ш е
(4)
где р = х^ + х - хр + (ат)к/(-тк), а т определяются также, как в (3) для г е (0, к].
По результатам оценки пограничной кривой Ф(х, г, ш) в узлах сетки {ац, % = 1,..., Ж} нетрудно вычислить ОПФ системы. В пашем случае р = (рх, 0) и имеем [4,8]:
Т (рх ,г,ш)=Та(-)ехр[гТф (•)],
где Та(рх,2,ш) = [Я2(^)+С2(■)]1/2 —амплитудно-частотная и Тф(рх,г,ш) = агсвт[5(■)/Та(■)] ^ фазочастотная характеристики системы, а Я(рх,г,ш) и С(рх,г,ш) соответственно синус- и косинус- преобразования Фурье от функции
дФ (■) ~
А(х, г,ш) = „ /[Фх,ш) - Ф(-то, г,ш)]. ох
Поскольку пограничная кривая вычисляется в конечном числе узлов, то точность вычисления функций Я(■) и С(■) существенно зависит от выбора сетки {х^}. Вблизи границы раздела яркости рекомендуется брать более мелкую сетку для более точного вычисления функции А(х, г, ш), поскольку она вычисляется методом аппроксимации. Кроме того, в данной работе при нахождении ОПФ используется дискретное прямое и обратное преобразование Фурье, результат которого является периодической функцией, период которой также напрямую зависит от выбора сетки и вычисляется по формуле Ь = 2-к/Ах, где Ах = - х¿. Поэтому необходимым условием для выбора сетки является то, что период полученной после преобразования Фурье функции включал в себя интересующий нас интервал области исследования оптической передаточной функции.
3 Результаты
Здесь приведен пример тестовых расчетов функций Та(■) с использованием описанной методики. Расчеты выполнены для длины волны А=0.5 мкм. Индикатриса рассеяния в атмосфере была принята рэлеевской ( [10], стр. 18), в воде использовалась линейная интерполяция индикатрисы, заданной таблично по результатам наблюдений для чистой воды ( [10], стр. 34).
Показатель преломления воды относительно воздуха был принят равным п = 1.33. Плотность распределения нормали к взволнованной поверхности определялась в соответствии с [9] двумерной гауссовской
плотностью распределения уклонов элементарных площадок с параметрами, зависящими от скорости ветра над поверхностью океана:
р(гх ) =
1
2пахау
ехр
2 -( М
\Ъх) \ Яу)
2 -
1Х = -zxД/г2 + + 1
1У = -zy+ ^ + 1,
^ + г2 + 1
к = 1/
а2х = 0.001 + 3.16 • 10-3 • V , а2у = 0.003 + 1.85 • 10-3 • V, V - скорость ветра и м/сек.
На рис. 2 представлены результаты оценок функции Та, полученные для скоростей ветра 2, 6, 10, 14 м/с и двух направлений. Синей линией обозначены графики функций для направления ветра поперек границы раздела яркости на дне, красной линией — вдоль границы раздела яркости.
Для качественного сравнения двух алгоритмов, описанных в работе [5], были произведены численные эксперименты при различных показателях скорости ветра.
Рис. 2: Функция Та для глубины 10 м в зависимости от скорости и направления ветра.
Рис. 3: Функция Та для глубины 10 м, скорости ветра У=2 м/с (слева) и У=6 (справа), вычисленная с использованием двух различных алгоритмов.
На рис. 3 представлены результаты оценок функции Та для скоростей ветра 2 м/с и 6 м/с, полученные с использованием двух разных алгоритмов: (1) — оценка, полученная с использованием формул (3) и (4); (2) — оценка, полученная с использованием алгоритма, исследованного ранее в работе [5]. В первом случае моделировалось lQ6 траекторий, а во втором — lQ7. Расчеты показывают, что первый алгоритм существенно превосходит по эффективности второй алгоритм, т.е. при одних тех же затратах времени вычислений статистическая погрешность расчетов первым алгоритмом значительно меньше статистической погрешности расчетов вторым алгоритмом.
Заключение
В работе получены результаты исследования поведения оптической передаточной функции для различных параметров окружающей среды. Результаты исследования показали, что ОПФ уменьшается с увеличением скорости ветра, а также очень сильно зависит от его направления. Также произведено сравнение двух возможных алгоритмов для оценки ОПФ. Результаты показали, что использование в объединении с остальными алгоритмами метода "сопряженных траекторий" при решении сопряженной задачи (в которой источником света считается морское дно) показывает более точную оценку ОПФ при одних и тех же затратах времени вычисления.
Список литературы
[1] Мулламаа Ю.-А.Р. Атлас оптических характеристик взволнованной поверхности моря. Институт физики и астрономии АН ЭССР. Тарту. 1964.Plass G.N., Kattawar G.W. Radiative Transfer in an Atmosphere-Ocean System // Appl. Opt. 1969. V. 8. № 2. P. 455-466.
[2] Plass G.N., Kattawar G.W. Radiative Transfer in an Atmosphere-Ocean System // Appl. Opt. 1969. V. 8, № 2. P. 455-466.
[3] Марчук Г. П., Михайлов Г. А., Назаралиев М. А. и др. Метод Монте-Карло в атмосферной оптике. // Новосибирск: Наука, 1976, с.280.
[4] Папулис А. Теория систем и преобразований в оптике. // М.: Мир, 1971, с. 495.
[5] Каргин Б.А., Ухинова О.С. Расчет методом Монте-Карло оптической передаточной функции системы океан-атмосфера // Труды международной конференции "Вычислительная математика и математическая геофизика", посвященной 90-летию со дня рождения академика А. С. Алексеева, ИВМиМГ СО РАН, Новосибирск. 8-12 октября, 2018 г., с.189-195.
[6] Каргин Б. А., Кузнецов С. В., Михайлов Г. А. Оценка методом Монте-Карло функции передачи контраста яркости через светорассеивающую среду. // Изв. АН СССР, ФАО, 1979, т.15, №10, с.1027-1035.
[7] Бори М., Вольф Э., Основы оптики.// М.:Наука, 1973, с. 720.
[8] Мишин И. В., Сушкевич Т.А. Оптическая пространственно-частотная характеристика атмосферы и ее приложения. // Исследования Земли из космоса, 1980, №3. с. 58-68.
[9] Сох С., Münk W. H., The measurement of the roughness of the sea curface from photographs of the sun's glitter. //J. Opt. Soc. America, 1954, 44, No. 11, p. 838-850.
[10] Кондратьев К.Я., Кантер P.P., Каргин Б.А., Поздняков Д.В. Численное моделирование в задачах оптического дистанционного зондирования внутренних водоемов Л.: Наука, 1987, 62 с.
Ухинова Ольга Сергеевна — к.ф.-м.н., науч.сотр. Института вычислительной математики и математической геофизики СО РАН;
e-mail: olsuÚosmf.sscc.ru.
Каргин Борис Александрович — д.ф.-м.н., зав. лабораторией Института вычислительной математики и математической геофизики СО РАН; Новосибирский государственный университет;
e-mail: bkarginÚosmf.sscc.ru.
Дата поступления — 30 апреля 2019 г.
