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

Перепишем еще раз краевую задачу для плоского слоямутной среды толщиной 0 при освещении его плоским мононаправленным(ПМ) источником в направлении ˆl0  1   ,0,   :200 dL( , ˆl )  L( , ˆl ) x(ˆl, ˆl) L( , ˆl)dˆl,d4 L( , ˆl )  (ˆl  ˆl 0 ), L( , ˆl ) 0.0,0  0 ,  0(2.1)Для того чтобы выделить особенность, представим решение в виде суммы40L( , ˆl )  La ( , ˆl )  Lr ( , ˆl ),(2.2)где La ( , ˆl ) – анизотропная часть (содержит особенность, и вычисляетсяаналитически), Lr ( , ˆl ) – регулярная (гладкая) часть, которая является гладкойфункцией, и может быть представлена конечной суммой.Подставив (2.2) в УПИ (2.1), получим новую краевую задачу для гладкойчасти решения Lr ( , ˆl ) Lr ( , ˆl ) x(ˆl, ˆl) Lr ( , ˆl ) dˆl  ( , ˆl );4 Lr ( , ˆl ) 0; Lr ( , ˆl )  La ( , ˆl ),0,0,00(2.3)dL ( , ˆl )( , ˆl )    a La ( , ˆl ) d4(2.4)где x(ˆl, ˆl) L ( , ˆl)dˆlaесть функция источников от невязки.

Изменения ГУ в (2.3) связано с тем, чтоприближенное решение не может удовлетворять точным ГУ. Выражение (2.4)нельзя использовать для прямых расчетов, поскольку La ( , ˆl ) содержит δфункцию, которая связана с прямым нерассеянным излучением. Выделим ее:La ( , ˆl )  e0 (ˆl  ˆl0 )  La ( , ˆl ),(2.5)что приведет выражение (2.4) к видуdL ( , ˆl )( , ˆl )    a La ( , ˆl ) d4 x(ˆl, ˆl) L ( , ˆl)dˆl  4 ea 0x(ˆl 0 , ˆl ).

(2.6)Поскольку La ( , ˆl ) является гладкой функцией, ее расчет не представляетпроблем [Budak et al., 2014].Важным и не очевидным является вопрос, какую часть тела яркостисчитать анизотропной? В угловом пространстве это сделать затруднительно,поскольку неясен критерий оценки величины производной, которую можнобыло бы считать соответствующей анизотропной части решения. Наиболееудачной здесь оказалась идея перехода к угловому спектру, обладающему41известным свойством: чем больше анизотропная угловая зависимость функции,тем более пологим является ее спектр. Следовательно, нужно понимать подLa ( , ˆl ) часть решения соответствующего наиболее медленной изменяющейсячасти углового спектра. Метод на основе данной идеи получил названиемалоугловой модификации метода сферических гармоник (МСГ), былразработан и подробнейшим образом описан в работах [Будак и Савенков,1980; Будак и Савенков, 1982; Будак и др., 1983; Будак и Сармин, 1990],поэтому мы не станем останавливаться подробно на выкладках, а приведемтолько окончательное выражение:dLa ( ,  , )  k 2k  1  d k Z k ( )Q km ( 0 )Q km (  )eim d0 k 0 m k 4где Qln (  ) 0e 02k  1xk Qmk ( 0 )Qmk (  )eim ,k 0 m k 4k(2.7)(l  n)! nPl (  ) – перенормированные полиномы Лежандра, Pln ( ) –(l  n)!присоединенные полиномы Лежандра, Z k ( )  e dk 0 e 0 , dk  1  xk .Однако важно отметить, что данный метод является наиболее общимсреди всех малоугловых методов [Будак и Федосов, 1985], а самое главное – онпозволяет наиболее точно выделить анизотропную часть решения, тем самымделая регулярную практически изотропной по углу функцией, что важно причисленной реализации [Kokhanovsky et al., 2010; Sokoletsky et al., 2014].

Исходяиз этих соображений, анизотропную часть решения в настоящей диссертациимы будем представлять именно в МСГ.Вследствие выделения анизотропной части решения, функция Lr ( , ˆl )будет являться гладкой по углу, а потому ее представление возможно поконечному базису. Например, если для дискретизации УПИ используется МДО,то такое представление будет иметь вид42ML r ( , i , )   (2   0,m )cos( m )C m ( , i ), C  ( )  C m ( , i ),(2.8)m 0C( )  C ( ), C ( ) гдеT–вектор-столбецдискретныхзначенийкоэффициентов разложения яркости в ряд Фурье по азимуту,  j  0.5( j  1) j,wj – нули и веса гауссовой квадратуры порядка N/2 для дискретизацииинтеграла рассеяния по зенитному углу.Это позволяет заменить интеграл рассеяния конечной суммой, а краеваязадача (2.3) преобразуется в краевую задачу для матричного неоднородноголинейного дифференциального уравнения первого порядка с постояннымикоэффициентами [Budak and Korkin, 2008; Budak et al., 2011; Budak et al., 2014]:d C( )  BC( )  M 1 ( ),d(2.9)где wB  M (1  AW), W   i4010 i, Mwi 0K0,x()(2k  1) xk Pk (cos  ),i k 0KA   (2k  1) Pkm ( i ) xk Pkm (  j ) .k 0Решение матричного уравнения (2.9) представимо в виде суммы общегорешения соответствующего однородного уравнения и частного решениянеоднородного [Budak et al., 2004; Budak et al., 2014; Fokina et al., 2014]:0 C(0)  P(0, 0 )C( 0 )   P(0, )M 1 ( , 0 )d ,(2.10)P(t , )  e B( t )(2.11)0гдеесть решение соответствующего однородного уравнения.

Она описывает связьраспределений яркости в двух точках t и  среды и носит название«пропагатор» [Budak et al., 2014].43Проблемы реализации решения (2.10) связаны с наличием в пропагаторе(2.11) как отрицательных, так и положительных экспонент, что физическисоответствует распространению потоков сверху-вниз и снизу-вверх. Этоприводит к быстрому ухудшению обусловленности матрицы и делаетневозможным расчеты для световых полей с оптическими толщами >1. Дляустранения этого эффекта используется масштабное преобразование [Karp etal., 1980]:0 SU C(0)  HU C( 0 )  S  et U 1M 1 (t, 0 ) dt,11(2.12)0где матричная экспонента представлена в виде разложения по собственнымвекторам и собственным значениям матрицы B [Sykes, 1951]:eB 0  Ue 0 U 1,(2.13)где U – матрица собственных векторов матрицы B ,   diag( ,  ) –диагональная матрица собственных значений, отсортированных так, что1e 00  u11 u12     , U  , H  , S    0uu 2122 0 e 010.1 Выражение (2.12) содержит только отрицательные экспоненты, чтоделает его устойчивым при любых значениях толщи.Для дальнейших расчетов необходимо вычислить интеграл от функцииисточников, стоящий в правой части выражения (2.12).

Для этого рассмотримвыражение для азимутальных коэффициентов функции источников, принимаяво внимание выражения (2.6) и (2.7) [Budak et al., 2014]: K 2k  1 m ( ,  )    1 d k Z k ( )Qkm ( 0 )Q km (  )  0 k m 4 e0  K 2k  11xk Qkm ( 0 )Qkm (  ). 0  k m 4(2.14)Дискретизируем это выражение в соответствии с выбранной схемойМДО:44 m ( )  M/ 0  1  2k4 1 (d ZKkk mk xk e 0)qk ( )Q km .(2.15)В соответствии с (2.12) и (2.15) функция источников примет вид [Budak etal., 2004]:02k  1qk (1  xk ) J k  J 0 U 1 1  0 M 1 Qmk , (2.16)k m 4KJ  S  e U M (t , 0 )dt  t110 где1J k  diag   i 0  d k dk  0 (He0 1  0 S), J 0  diag  (He  i 0  1 0 S).

(2.17)Вектор-столбцы C (0), C ( 0 ) в выражении (2.12) описывают потокиизлучения, падающего на слой, и определяются граничными условиями.Векторы C (0), C ( 0 ) описывают потоки отраженного и прошедшего слойизлучения. Разрешим уравнение (2.12) относительно потоков выходящих изслоя и получим решение в форме рассеивателей [Budak et al., 2011; Budak et al.,2014]: C (0)   F   R T   C (0)  ,FTRC()C()   0   0     F где     hJ, F  u11  u12 e 0 u12 R T h,h    0 0TReuuu 222122  e(2.18)1e 0 u11  .u 21 Выражение (2.18) представляет собой связь выходящих из слоя потоков спадающими.

Матрицы R и T содержат дискретные значения коэффициентовяркости по отражению и пропусканию соответственно, а вектор F описываетсобственное излучение слоя, обусловленное невязкой МСГ с точным решениемУПИ.В работе [Будак и др., 2011] показано, что выделение анизотропной частии дискретизация УПИ всегда приводит одному и тому же уравнению УПИ дляслоя (2.9), которое имеет единственное решение либо через пропагатор (2.10),45либо в эквивалентной форме рассеивателей (2.18). Это подтверждаютсравнения между собой известных алгоритмов [Kokhanovsky et al., 2010;Sokoletsky et al., 2014]: все результаты совпали с машинной точностью, а этоозначает, что сравнение проводилось для различных реализаций одного и тогоже решения.

Различия составляли лишь методы выделения анизотропной части,о чем подробнее будет написано в параграфе 3 настоящей главы. Сравнениепоказало, что выделение анизотропной части на основе МСГ позволяетсократить время расчета более чем на два порядка.Сложность расчетов по полученному решению (2.18) определяетсяразмерами входящих в него матриц, которые зависят от трех целых чисел: N, Mи K, где N – количество дискретных ординат, M – количество азимутальныхгармоник, K – количество членов разложения индикатрисы по полиномамЛежандра.В общем случае все эти три величины примерно равны, однако приудачном выделении анизотропной части можно добиться, что угловаязависимость регулярной части будет близка к изотропной.

В этом случаеK >> N >> M,что существенно ускоряет решение.Расчеты программы MDOM в двух случаях: для одного и того жезначения K, N1=4N2, M1=32M2 показывают практически полное совпадениераспределений яркости по углу за исключением нескольких острых пиков. Приэтом время счета разнится более чем в 200 раз. Это связано с тем, чторегулярная часть решения действительно близка к изотропной, но с некотороймелкой рябью.

Встает вопрос: сколько дискретных ординат требуется дляописания этой ряби?Поскольку регулярная часть решения является гладкой функцией, то ееразложение в ряд по полиномам Лежандра имеет конечное число членов N:2k  1 mLk Pk (  ).2k 1NLm (  )  46(2.19)Все полиномы Лежандра порядка менее чем N точно выражаются черезполином N+1 порядка:N 1Pk (  )   Pk ( i )i 1PN 1 (  ),(   i ) PN 1 ( i )(2.20)где i корни полинома PN+1().Соответственно, это приводит к выражениюNN 1PN 1 (  )2k  1 mPN 1 (  )Lm ( ,  )  Lk Pk ( i )   Lm ( , i ), (2.21) 1 ( i ) k 0 2(   i )PN 1 ( i )i 1 (   i )PNi 1N 1что соответствует интерполяционной формуле Лагранжа для функции L(,)[Budak et al., 2015].Эти соотношения являются аналогом теоремы Котельникова о выборкахуглового спектра углового распределения по полиномам Лежандра. Отсюдаследует, что для сходимости в равномерной метрике, шаг дискретизациидолжен соответствовать половине углового размера самой мелкой детали,которую необходимо воспроизвести в распределении яркости.Для решения обратных задач требуется алгоритм решения УПИ вравномерной метрике не хуже одного процента со временем счета не более 1секунды для одной длины волны.