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

MKL используютспециализированные математические пакеты, такие как Matlab, Maple,Mathematica [Будак и др., 2011].Существенноеускорениеалгоритмовиболеерациональноеиспользование памяти возможно с использованием алгоритмов разреженныхматриц (sparse matrices). Это эффективно при работе с матрицами, у которыхзначение большого числа элементов являются нулями. При этом хранятсятолько ненулевые элементы и дополнительная информация об индексах.Существует несколько форматов хранения разреженных матриц: compressedsparse row (CSR) format, compressed sparsed column (CSC) и Coordinate Format.Выбор формата определяется типом операции с матрицами.Программа MVDOM реализована на двух языках – Matlab и Fortran.Влияние различных факторов на время счета программы для двух тестовприведено в таб.

1. Тест №1: N = 101, K = 500, M = 32; тест №2: N = 101, K =1000, M = 32. Тесты выполнялись на компьютере IntelCore 2 Duo 3 ГГц, 2 ГбRAM, Intel Fortran Compiler 11.1 с библиотекой MKL 10.2. Использовались двакомпилятора – gfortran и ifort.62Условия реализации кодаВремя счета Времятеста №1, ссчетатеста №2, сgfortran + LAPACK240530gfortran + LAPACK+optimization230505ifort + LAPACK210490ifort + MKL115250ifort + MKL + optimization105230ifort + MKL + optimization + sparse matrices3344MATLAB 2010b2745MATLAB 2010b + CUDA2233Таблица 1 - сравнение времени счета для двух тестовТакое влияние объясняется тем, что все тесты относятся к аэрозольномурассеянию,имеющемублочно-диагональнуюматрицурассеяния,чтоопределяет не менее половины элементов нулевыми.Такжебылопроведеноотдельноеисследованиеэффективностииспользования алгоритма разреженных матриц: зависимость времени счета отK.

Благодаря этому алгоритму операции с двумерными массивами становятсяэквивалентно одномерным, что определяет линейную зависимость временисчета от K в отличие от K2 в обычном случае [Будак и др., 2011].Перспективным средством ускорения кода является использованиеграфическипроцессоровобщегоназначения.Программно-аппаратнаяархитектура (Compute Unified Device Architecture, CUDA) позволяет проводитьвычисления на графических процессорах nVIDIA GPUs, а это в свою очередьоткрываетвозможность проведениястандартныхпроцедурбиблиотекиLAPACK в рамках архитектуры CUDA. MATLAB 2010b поддерживает nVIDIACUDA GPUs, однако достоинства вычислений на GPU существенны при63использовании больших массивов, в обратной ситуации вычисления на CPUявляются более эффективными.В коде MVDOM использовался алгоритм выделения анизотропии наоснове МСГ, что существенно снизило размеры всех массивов.

Притестировании не приходилось брать число N более 300. Оценка времени счета спомощью профайлера показывает, что половина времени уходит на нахождениесобственных значений. Однако эта процедура на реализована в пакетеMATLAB GPU. Нам удалось снизить время выполнения программы на 20% засчет перемножения матриц на nVIDIA GeForce 480 GTX GPU [Будак и др.,2011].Из проведенного анализа следует, что прежде всего скорость вычисленияалгоритма определяется методом выделения анизотропной части решения.МСГ практически точно определяет излучение вперед, при этом оставляетрегулярную часть практически изотропной по углу, с мелкими пиками,несущественными с точки зрения энергетики.

Поэтому регулярную частьпредставляется возможным с высокой точностью вычислить с помощьюпростейшего метода решения УПИ – двухпотокового приближения. При этомматрицы из сотен элементов, используемые в классических методах решения,превратятся в матрицы размера 2х2.Разумеется, использование самого грубого приближения приведет к«размазыванию»этихнезначительныхпиков.Однакопоследующееиспользование итерации позволит с небольшим увеличением времени счетаучесть их, и получить хорошую сходимость в равномерной метрике.64Выводы по второй главе1.

Для эффективного решения обратных задач ОДЗ требуется реализацияматематической модели сигнала спутниковой оптико-электронной системы ДЗЗв равномерной метрике не хуже одного процента с временем счета не болееодной секунды для одной длины волны. Лучший на сегодняшний деньалгоритм MVDOM не удовлетворяет требованиям по скорости более чем на трипорядка. Поскольку MVDOM фактически является точным решением УПИ, этопозволяет говорить о несостоятельности классических методов решения с точкезрения современных гиперспектральных измерительных систем, а значит,требуется искать принципиально новый подход.2. Возможным способом ускорения сходимости алгоритма модели сигналаоптико-электронной системы ДЗЗ является метод синтетических итераций,который позволяет снизить количество ординат в МДО, и тем самымсущественно повысить скорость вычислений без потери точности.

В работевпервые данный метод был применен для описания переноса излучения.3. Поскольку МСГ позволяет получить регулярную часть решения УПИ ввиде практически изотропной по углу функции, это дает возможность вкачестве первого шага синтетических итераций использовать простейшийметод решения УПИ – двухпотоковое приближение, тем самым получитьсущественный (несколько порядков) выигрыш в скорости без потери точности.4. Наиболее существенным фактором, влияющим на скорость алгоритма,является метод выделения анизотропной части решения. При использованииМСГ регулярная часть получается практически изотропной по углу функцией.Это позволяет для ее описания применять с хорошей точностью простейшиеметоды решения УПИ, избегая использования матриц из тысяч элементов, темсамым сократить время вычисления на несколько порядков.

Последующаяитерация при небольшом увеличении времени счета позволяет добитьсяхорошей сходимости в равномерной метрике.653. МЕТОД СИНТЕТИЧЕСКИХ ИТЕРАЦИЙ3.1. Квазидвухпотоковое приближениеПод квазидвухпотоковым приближением мы будем понимать частныйслучай метода синтетических итераций, когда в качестве первого шага(приближенного метода) мы будем использовать простейший метод решенияУПИ – двухпотоковое приближение (метод Шварцшильда-Шустера).

Его идеязаключается в усреднении УПИ по двум полусферам направлений впространстве – верхней и нижней, что приводит УПИ к системе из двухобыкновенных линейных неоднородных дифференциальных уравнений. Дляэтого снова подставим полное решение в виде суммы прямой составляющей ирегулярной частиL( ,  )  L0 ( ,  )  L( ,  )(3.1)в исходное УПИ для плоского слоя, что приведет его к виду:L( ,  ) 4 x(ˆl, ˆl )L( ,  )dˆl  Q( ,  , ),(3.2)гдеQ( ,  ,  ) 4 x(ˆl, ˆl)L ( ,  )dˆl  L ( ,  )  00L0 ( ,  )(3.3)есть функция источников.Прямое излучение по определению удовлетворяет УПИ без рассеянияL0 ( ,  ) L0 ( ,  )  0,(3.4)решение которого имеет видL0 ( ,  )   (ˆl  ˆl0 )e /  .Подставим выражение (3.5) в функцию источников (3.3)Q( ,    ) 4 (ˆl  ˆl 0 )e /  /  ˆ / ˆˆˆˆˆˆ x(l, l ) (l  l0 )e dl   (l  l0 )e  66(3.5) ˆ ˆ  / 01 ˆ ˆ  / 0x ( l , l 0 )e  (ˆl  ˆl0 )e /    e /  (ˆl  ˆl 0 ) x ( l , l 0 )e.44(3.6)С учетом этого УПИ для гладкой части решения примет видДляL( ,  ) L( ,  ) 4переходак x(ˆl, ˆl)L( ,  )dˆl  4 x(ˆl, ˆl )e / 00двухпотоковомуприближению(3.7).проинтегрируемуравнение (3.7) по нижней и верхней полусферам направлений в пространствеL( ,  ) ˆdl   L( ,  )dˆl  4  L( ,  ) dˆl  L( ,  )dˆl  4  x(ˆl, ˆl )L( ,  )dˆldˆl  4 e / 0 x(ˆl, ˆl )dˆl,0ˆl, ˆl )L( ,  )dˆldˆl   e  / 0 x(ˆl , ˆl )dˆl ,x( 04(3.8)где Ω+ и Ω- означают нижнюю и верхнюю полусферы соответственно.Рассмотрим первое уравнение системы и проинтегрируем каждый членпо отдельностиL( ,  ) ˆ dl   L( ,  )dˆl  E ( ), L( ,  )dˆl  E ( ),(3.9)(3.10)44 x(ˆl, ˆl) L( , )dˆldˆl  4ˆl, ˆl ) L( ,  )dˆl  x(ˆl, ˆl ) L( ,  )dˆl dˆl x(  ˆˆˆˆˆˆˆˆx(l,l)dlL(,)dlx(l,l)dlL(,)dl    c E ( )   o E ,     / 0ex(ˆl, ˆl0 )dˆl c 0e / 0 .44(3.11)(3.12)Здесь μ↓ – средний косинус угла рассеяния в нижнюю полусферу, E ( ) –полусферическая облученность в нижнюю полусферу, βc и βo – коэффициентырассеяния в совпадающую и противоположною направлению падения67излучения полусферу соответственно, βс0 – доля излучения, падающего впереднюю полусферу.Проделаем то же самое со вторым уравнением системы (3.8):L( ,  ) ˆ dl   L( ,  )dˆl  E ( ), L( ,  )dˆl  E ( ),(3.13)(3.14)44  x(ˆl, ˆl)L( , )dˆldˆl  4ˆl, ˆl ) L( ,  )dˆl  x(ˆl, ˆl ) L( ,  )dˆldˆl x(  ˆˆˆˆˆˆˆˆx(l,l)dlL(,)dlx(l,l)dlL(,)dl      o E ( )  c E ( ),  / 0ex(ˆl, ˆl0 )dˆl o 0e / 0 ,44(3.15)(3.16)где μ↑ – средний косинус угла рассеяния в обратную (верхнюю) полусферу,E ( ) – полусферическая облученность в обратную полусферу, βo0 – доляизлучения, падающего в обратную полусферу.Дляопределениякоэффициентов рассеянияидолейпадающегоизлучения представим индикатрису рассеяния в виде разложения в ряд пополиномам Лежандра [Budak et al., 2010]x(cos  )   (2k  1) xk Pk (cos  ),k 0где γ – угол рассеяния.Имеем выражения1c 41x(ˆl, ˆl ) dˆlˆl681  x(  ) d  20(3.17)11 K  (2k  1) xk  Pk (  )d 2 k 00(3.18)иc 0 14 x(ˆl , ˆl)dˆl.(3.19)0Для полиномов Лежандра известна теорема сложенияm4 mk m ˆˆˆPk (l 0  l) Y(l)Y k 0 k (ˆl),2k  1 m k(3.20)где черта сверху означает комплексное сопряжение, а2k  1 (k  m)! m ˆYkm (ˆl ) Pk (l  zˆ )eim .4(k  m)!Тогда (3.19) примет видc 0   xkk 0mkYm  kmkm(ˆl0 )  Y k (ˆl)dˆl 22k  1 (k  m)! im  m  xk  Y (ˆl 0 )ePk (  )d  d .2(km)!k 0m  k00mk1mk(3.21)Поскольку2e im d   2 m 0 ,0то1c 0  0.5 (2k  1) xk Pk ( 0 )  Pk (  )d  .k 0(3.22)0Интеграл, входящий в выражения (3.18) и (3.22) вычисляется и имеет значения:1) от всех четных k≠0 он равен нулю;11002) при k=0  P0 (  )d    d   1 ;692(1) j 1 2 j  1!   , что3) при нечетных k=2j-1  P2 j 1 (  )d   2 jj2 (2 j  1) j  j  1! j  1!01следует из общей формулы1P2 p 22l 2( ) P2 p 22 j 102(1) j l 1(  )d   2(l  j 1)2(2l  2 j  1)(l  j  1) 2( j  p)  3! 2(l  p)  4 ! l  p ! j  p ! l  p  2 ! j  p  2 ! p1, l 1(1) j 1 2 j  1! . 2 j 12 (2 j  1) j  j  1! j  1!Однако с точки зрения численной реализации лучше избавиться отфакториалов и перейти к рекуррентному соотношению: j 1(1) j 2 2 j  1! : (1) j 1 2 j  1!  2 j 12 j 1 j 2 (2 j  1)( j  1) j ! j ! 2 (2 j  1) j  j  1! j  1!(2 j  1) j 2 j  2 j  12 j 12 j 1  j 1   j,24(2 j  1)( j  1)j2( j  1)2( j  1)(3.23)2(1) 2 1!1причем 1  2 .2  1  1 0!0! 2Следовательно1k 00 c 0  0.5 (2k  1) xk Pk ( 0 )  Pk (  )d    0.5 1   (4 j  1) x2 j 1 P2 j 1 ( 0 ) j  ,j 1(3.24) ( K 1)/2c  0.5 1   (4 j  1) x2 j 1 0j  ,j 1(3.25)иo0  1  c0 ,o  1  c .Для однопараметрической индикатрисы Хеньи-Гринстейна70x2 j 1  g 2 j 1,(3.26)и11 g2c  d 2 0 (1  g 2  2 g  )1.51 g11 g2 1 g1 g2.(3.27)В результате интегрирования по полусферам получаем систему из двухобыкновенныхдифференциальныхнеоднородныхуравненийдляполупространственных облученностей E ( ) ( c  1) E ( )   o E ( )  c 0e / 0 , 4  E ( )   E ( )  (  1) E ( )    e  / 0 ,o co0  4(3.28)или в матричной формеMгде M   0d E( )  A E( )  e / 0 F,d(3.29)0 E ( ) 1   c  o  c 0 AF,,,E()   . 1E()oc o0   Решение будем искать методом вариации произвольной постоянной.Решение соответствующего однородного уравнения имеет видE( )  exp( B )  C,(3.30)где B  M1 A .Представим константу в виде функции от τ E ( ) E( )  exp( B )  0   exp( B ) E 0 ( ) E0 ( ) и подставим ее в уравнение (3.29)Md exp( B ) E 0 ( )  Aexp( B ) E 0 ( )  e / 0 F d71(3.31)d E 0 ( )  / 0M   Bexp( B ) E 0 ( )  exp( B )F   A exp( B ) E 0 ( )  edd E 0 ( ) Aexp( B ) E 0 ( )  Mexp( B )  Aexp( B ) E 0 ( )  e / 0 F dd E 0 ( ) e / 0 exp(B ) M 1 F.d(3.32)Полученная система уравнений является системой с разделяющимисяпеременными, следовательноd E 0 ( )  e / 0 exp(B ) M 1 F d  e 1 / 0 e B M 1 F d  e(B 0 1) / 0 M 1 F dиE 0 ( )  E 0 (0)    e(B 0 1) t/ 0 dt M 1 F  0 (B 0  1) 1 e(B 0 1) / 0  1 M 1 F 0 0 (1  B 0 ) 1 1  eB e / 0  M 1 F.(3.33)Умножив полученное решение на решение однородного уравнения,получим полное решение системы (3.29)E( )  e B E 0 (0)  0e  B (1  B 0 ) 1 1  e B e  / 0  M 1 F.(3.34)Неизвестные константы найдем, подставив в уравнение (3.34) известныеграничные условия E (0)  0, E ( 0 )  0,(3.35) E ( 0 )   B  0  B1B  / 0 M 1 F.  e  E (0)   0e (1  B 0 ) 1  e e   0 (3.36)Снова представим матричную экспоненту через собственные векторы исобственные значения матрицы Be B  Ue U1 .(3.37)где   матрица собственных значений, а U  собственных векторов матрицыB:72 E ( )  0 e 0 U 1   0   U 1  0 U 1 (1  B 0 ) 1 1  eB e / 0  M 1 F.