Диссертация (Быстрые и точные алгоритмы расчета сигналов оптико-электронной системы дистанционного зондирования из космоса), страница 5

При сильно анизотропном теле яркости путь,пройденный рассеянными фотонами, мало отличается от пути, пройденногонерассеянными, что позволяет пренебречь дисперсией путей рассеянныхфотонов. Главной проблемой МУП является расчет многомерной свертки насфере индикатрисы рассеяния и яркости падающего излучения, для которойсуществуют три основных пути решения [Будак, 2007]. Первый подход былпредложен Боте, в нем процесс рассеяния частиц представляется марковскойцепью блужданий, и помимо пренебрежения дисперсией фотонов такжепредполагает большую остроту индикатрисы, по сравнению с телом яркости.Поэтому такая форма МУП не решает главную задачу – строгий учетособенностей решения.В работе [Goudsmit and Saunderson, 1940] был предложен другой подход,основанный на теореме сложения для полиномов Лежандра, которыйпренебрегает только дисперсией путей рассеянных фотонов, но при этом неможет быть обобщен ни на какие другие источники, кроме плоскогомононаправленного.В работах [Snyder and Scott, 1949] и других авторов была предложенатретья форма МУП, которая наряду с пренебрежением дисперсии путейрассеянных фотонов предполагала возможность замены свертки на сфере вобласти малых углов сверткой на плоскости.В работах [Lewis, 1950; Wang and Guth, 1951] рассматривалась связь всехтрех подходов, и было проведено их сравнение между собой по точности.Показано, что наиболее точно яркость под малыми углами учитывает второйподход, а точность первого и третьего существенно падает с ростом угла.В 1982 году был предложен новый малоугловой метод на примереточечного изотропного источника, основанный на методе сферическихгармоник (СГ) – малоугловая модификация метода сферических гармоник23(МСГ) [Будак,Савенков, 1982].

В этом случае яркость и индикатрисарассеяния раскладываются в ряд по полиномам Лежандра2k  12k  1Ck (r ) Pk (), x(ˆl, ˆl)  xk Pk (ˆl  ˆl) ,44k 0k 0L ( r , )  (1.8)где Ck(r) – коэффициенты разложения, Pk(μ) – полиномы Лежандра, xk –коэффициенты разложения индикатрисы.После чего делается предположение о том, что коэффициентыразложения яркости и индикатрисы являются непрерывными функцияминомера гармоники в предположениилинейной аппроксимациимеждупредыдущим и последующим коэффициентами.

В результате этого бесконечнаясистема связанных уравнений заменяется одним уравнением математическойфизики.Показано, что решение при использовании такого подхода описывает всеособенности решения УПИ [Будак и др., 1983; Будак и др., 1984]. В статье 1985года [Будак и др.] проводится сравнение МСГ с известными формами МУП, ипоказано, что подход через МСГ является наиболее общим, пренебрегающимтолько дисперсией путей рассеянных фотонов и обратным рассеянием, и другиеформы МУП следуют из него. В виду этого для выделения анизотропной частирешения в настоящей работе мы будем использовать именно подход МСГ,который применительно к предлагаемому методу будет подробно описан вовторой главе.Одним из исторически первых методов, наряду с двухпотоковымприближением, был метод «обычных», не модифицированных сферическихгармоник (СГ).

В этом случае искомую яркость и индикатрису рассеяния такжераскладывают в ряд по сферическим функциям, что приводит к бесконечнойсистеме дифференциальных уравнений с постоянными коэффициентами [Karpet al., 1980]. Для решения систему делают конечной, полагая все коэффициентыразложения, начиная с номеров больше некоторого N, равными нулю. Прирешении анизотропной части в МСГ, а гладкой – с помощью СГ происходит24выигрыш в скорости по сравнению со сферическими гармониками более чем в360 раз. Однако представление гладкой части в виде ряда по сферическимфункциям имеет свои недостатки. Плохо обусловленная матрица приводит косцилляциям в решении, а также в этой форме является неудобнымпредставление граничных условий [Budak et al., 2004].Среди методов решения УПИ нельзя не остановиться на методедискретных ординат (МДО), который разработал С.

Чандрасекар, развиваяидею двухпотокового приближения Шварцшильда-Шустера. Идея заключаласьв том, чтобы уточнить слишком грубый метод большим количествомусредненных потоков для возможности получения детальной информации обанизотропии и поляризации. Это переводит УПИ в конечную системуобыкновенных дифференциальных уравнений [Чандрасекар, 1953]idL(, i )1 L(, i )   a j L(,  j )  q(, i ), (i  1,...,  n) ,d2 j(1.9)где aj – веса в квадратурной формуле, основанной на делении интервала (1,+1)точками μj.

Обычно эти точки представляют собой нули полинома Лежандра mной степени, при этом предпочтение обычно отдают формуле Гаусса,поскольку в этом случае при одном и том же значении m можно получитьточность почти вдвое большую, по сравнению с другими квадратурнымиформулами [Чандрасекар, 1953].Этот метод стал представлять все больший интерес с ростомвычислительной мощности компьютеров, поскольку, будучи ограниченнытолько ею, увеличивая количество потоков, в пределе мы фактически получимточное решение УПИ.Компьютерная реализация МСГ для анизотропной части решения и МДОдля регулярной породила программу MVDOM (Modified Vectorial DiscreteOrdinate Method), созданную научной группой кафедры светотехники МЭИ подруководством профессора Владимира Павловича Будака.

В работе [Будак и др.,2011]получено, что после выделения анизотропной части решения и25дискретизации, УПИ имеет единственное решение в матричной форме. Этозначит, что ключ к ускорению алгоритма стоит искать в как можно болееточном выделении анизотропной части решения, чтобы получить регулярнуюкак можно более гладкой функцией, и иметь возможность использовать болеекрупную сетку. В этой работе также было проведено сравнение даннойпрограммы с аналогами с использованием различных математическихбиблиотекисредпрограммирования.Результатыпоказали,чтодляотносительно гладких индикатрис рассеяния другие методы не уступают МСГ,однако, в случае острых индикатрис МСГ является наиболее удачным методомвыделения анизотропной части решения, который позволяет получитьрегулярную часть практически изотропной по углу функцией.

В вопросахскорости наилучший результат был достигнут с использованием средыпрограммирования Matlab 2010b, в котором используется библиотека MKL.Также сравнение MVDOM было проведено с известными аналогами в работах[Boudak and Korkin, 2008; Kokhanovsky et al., 2010; Sokoletsky et al., 2014],которое показало, что эта программа является лучшей, однако и она неудовлетворяет требованиям гиперспектральных систем по скорости системпочти на два порядка.Среди методов решения УПИ нельзя не коснуться группы статистическихметодов, таких как методы Монте-Карло.

Общий принцип методов МонтеКарло основан на моделировании такой случайной величины, статистическиехарактеристики которой представляли бы собой искомые параметры исходнойзадачи [Соболь, 1978]. В случае решения УПИ методом Монте-Карло, процессраспространения света представляется как случайная марковская цепьстолкновений фотонов с веществом, которые приводят либо к рассеянию, либок поглощению фотонов.

Метод Монте-Карло заключается в моделированиитраекторий этой цепи и вычислении статистической оценки для искомыхфункционалов. Построение случайных траекторий для физической моделипроцесса называется прямым моделированием. Математическая суть данного26методасостоитвнахожденииоптимальныхспособовмоделированияслучайных величин [Марчук и др., 1976].Сложные задачи переноса излучения, как правило, невозможно решитьметодом прямого моделирования с необходимой точностью.

Для этогоразрабатываются методы уменьшения дисперсии оценок и алгоритмы длярасчетов специальных функционалов методом Монте-Карло. Эффективностьметодов уменьшения дисперсии сильно зависит от специфики конкретнойзадачи [Марчук и др., 1976].Несмотря на широкое применение метода Монте-Карло в задачахпереноса излучения, мы не будем останавливаться на нем подробно, посколькудлядостиженияточности,требуемойгиперспектральнымисистемами,потребуются большие временные затраты, и требования по скорости не будутудовлетворены.Важное место среди методов решения УПИ как вообще, так и в рамкахнастоящей диссертации, занимает подход, при котором искомое поле яркостипредставляется в виде суммы по кратностям рассеяния – метод итераций.

Дляэтого необходимо перейти к интегральному уравнению Пайерлса, интегрируяУПИ вдоль луча. Сделаем переход к переменным вдоль луча из точки R понаправлению l̂ . При этом произвольная точка r на луче на расстоянии  от Rвдоль луча должна удовлетворять уравнениюr  R   ˆl .(1.10)В этом случае УПИ в переменных вдоль луча (1.10) будет иметь видdL(R  ˆl, ˆl )  L(R  ˆl, ˆl ) x(ˆl, ˆl) L(R  ˆl, ˆl)dˆl .d4(1.11)Если считать правую часть уравнения (1.11) известной функцией, то онобудет являться линейным неоднородным дифференциальным уравнениемотносительно переменной , решение которого имеет вид270  (0  )L(r, ˆl )  L0 e4 0 x(ˆl, ˆl) L(r   ˆl, ˆl)d ˆld ,(1.12)где 0 – расстояние по лучу до границы объема, а L0  Le (r  0ˆl , ˆl )e 0 – прямоеизлучение, не испытавшее ни одного рассеяния и ослабленное по законуБугера.

То же самое можно представить в операторной формеL  L0TL ,(1.13)где T  D1S – оператор переноса вдоль луча.На основе принципа сжимающихся отображений решение (1.13) можнопредставить в виде ряда Неймана [Краснов, 1975; Будак, 2007]L   T n L0 ,(1.14)n 0что физически и означает разложение по кратностям рассеяния света.Соответственно яркость произвольной кратности выглядит следующимобразом:Ln (z, rˆ , ˆl )  TLn1  T n L0    4 n 11n........exp( k )  x(ˆl, ˆl1 )  ... 0 0   k 1n(1.15)n x(ˆl n1 , ˆl n ) L0 (r, ˆl n )dˆl1...dˆl n d1...d n .Следует отметить, что все программы, с которыми проводилосьсравнение MDOM, вне зависимости от подхода к решению на завершающемэтапе используют метод итераций.

Это носит название пост-процессорнойобработки решения для аппроксимации решения на произвольной сетке угловвизирования [Thomas and Stamnes, 2002]. При этом даже одна итерациясущественно уточняет решение, сглаживая его. Этот факт явился источникомидеи применить метод итераций к программе MDOM.