Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.astro.spbu.ru/staff/ilin2/INTAS2/PAP/P4_4.pdf Дата изменения: Fri Nov 19 16:18:48 2010 Дата индексирования: Tue Oct 2 03:02:14 2012 Кодировка: Windows-1251 Поисковые слова: п п п п п п п п п

Вычисление параметрических производных поляризованного излучения и решение обратных задач атмосферной оптики
С. А. Ухинов



и Д. И. Юрков

14 апреля 2003 г.

Институт Вычислительной Матаматики и Математической Геофизики СО РАН, Новосибирск, 630090, Россия

Аннотация
Авторами рассмотрены вопросы вычисления производных от параметров поляризованного солнечного излучения и вопросы решения обратных задач атмосферной оптики с использованием модифицированной двойной локальной оценки метода Монте-Карло при наблюдениях с Земли. Исследована применимость итерационного процесса Ньютона-Канторовича и предложен метод, основанный на минимизации квадрата нормы разности вычисленных и наблюдаемых характеристик поляризованного излучения.

Целью настоящей статьи является продолжение разработки методов Монте-Карло для решения задач теории переноса излучения с учетом поляризации. Интерес вызывает исследование возможностей решения обратных задач атмосферной оптики по результатам наземных измерений поляризационных характеристик рассеянного солнечного излучения. В частности разрабатываются алгоритмы восстановления высотного хода коэффициентов аэрозольного рассеяния по измерениям степени поляризации излучения с вычислением соответствующих параметрических производных методом Монте-Карло. Сложность поставленной задачи состоит в том, что наземные измерения информативны, по-видимому, только при зенитых углах Солнца больших 90 (в сумерках). Погрешности же метода Монте-Карло в сумеречных расчетах на порядок больше погрешностей расчетов в дневной области Земли, а погрешности расчетов параметрических производных от функционалов на порядок больше погрешностей расчета самих функционалов. Поэтому важными этапами решения обратных задач являются как разработка и исследование алгоритмов расчета соответствующих параметрических производных, так и анализ оценок их погрешностей. В работе строятся алгоритмы определения коэффициентов аэрозольного рассеяния с использованием итерационного метода метода Ньютона-Канторовича и метода Монте-Карло для вычисления коэффициентов возникающих при этом систем линейных алгебраических уравнений. Предложен метод, аналогичный методу покоординатного спуска, полученный из принципа минимизации квадрата нормы разности вычисленных и наблюдаемых характеристик поляризационного излучения. Численно исследуется применимость и эффективность этих методов. Проводится их сравнительный анализ.

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты E-mail: sau@osmf.sscc.ru

IR 00-01-00797, 00-15-96173), Интеграционного гранта СО РАН 2000 No 43, ФЦП Интеграция.

1


1

Общие положения
L1 (X ): i = 1, ..., 4, L (X ): (1.2) (1.1)
4

Рассмотрим уравнение переноса с учетом поляризации [1,2] в

i (x) =
или в операторной форме

j =1 X

Kij (x , x)j (x )dx + Si (x),

= K + S
4

, и сопряженное к нему уравнение в

(x) =
j =1 X
или в операторной форме Здесь и

i

Kj i (x, x ) (x )dx + hi (x), j
.

i = 1, ..., 4,

= K + H

X = RЧ

, где

R

пространство координат,

пространствах вектор-функций

L1 (X ) и L (X ) = v raisupi,x |hi (x)| соответственно. (X ) Оператор K описывает изменение вектор-функции композицией операторов l и C [2]: K = C l, ||H ||L
где скалярный оператор плотности столкновений оператора

пространство направлений. Нормы в 4 задаются равенствами ||F ||L1 (X ) = i=1 X |fi (x)|dx плотности рассеяных фотонов и является



l

описывает изменение вектор-функции



при перемещении фотона от

точки рассеяния до точки столкновения, а оператор

C

описывает изменение вектор-функции

l

при изменении направления движения фотона. Следовательно ядро

K

имеет вид:

K (x , x) =
Здесь

exp{- (r, r )}s (r)P ( , , r) ( - s) |r - r |2

(1.3)

s (r) коэффициент рассеяния, (r) = c (r)+s (r) коэффициент полного взаимодействия |r - r | со средой, c (r ) коэффициент поглощения; (r, r ) = 0 (r + st)dt оптическая толщина отрезка (r , r ); s = (r - r )/|r - r |, P ( , , r ) угловая матрица рассеяния. Коэффициент рассеяния s может быть представлен в виде суммы аэрозольного a и молекулярного m коэффициентов рассеяния: s = a + m . Соответственно, физическая матрица рассеяния
имеет вид:

R(a , m , ч, r) =
Здесь

Ra (ч, r)a + Rm (ч, r)m . a + m

R

aи

R

m матрицы аэрозольного и молекулярного рассеяния соответственно.

Вектор-функция

S (x)

является физической плотностью распределения источника; е? конкретный

вид зависит от реальных источников излучения. Вектор-функция плотности рассеяний

(x)

связана с классической вектор-функцией Стокса

I (x),

описывающей интенсивность, степень поляризации, плоскость поляризации и степень

эллиптичности излучения [1] следующим образом:

I (x) =
при этом компоненты стоксовской. Будем предполагать, что оператор

[l](x) . (r)

(1.4) 1 0, 2 + 2 + 2 2 3 4

2 1. В дальнейшем произвольную функцию, удовлетворяющую этим условиям, будем называть



i функции

(x)

обладают свойствами

S

t вектор-функций Стокса, т. е. если

K является положительным St , то K St .

на воспроизводящем конусе

2


Для решения задачи вычисления линейного функционала IH = , H = , S , где , H = T (x)H (x)dx, , S = S T (x) (x)dx, строится оценка по рассеяниям метода МонтеX X Карло [3] на однородной цепи Маркова x0 , ..., xN с начальной плотностью p(x , x), имеющая
вид:

N

=
n=0
где векторные веса

T qn H (xn ),

(1.5) Qn
по рекуррентным формулам:

q

n определяются через матричные веса

q0 =

S (x0 ) , (x0 )

Q0 = {ij }, qn = QT q0 , n (x)

Qn = Qn-

1

K T (xn-1 , xn ) , p(xn-1 , xn )

n = 1, ..., N .
и переходную

(1.6) p(x , x)
плотности

Для обеспечения несмещенности оценки

на начальную

цепи Маркова следует наложить следующие условия несмещенности:

(x) = 0 p(x , x) = 0
Нетрудно заметить, что

S T (x) (x) = 0, K T (x , x) (x) = 0.

T = q0 x

0

, где

N

x0 =
n=0

Qn H (xn ). K

(1.7)
выполняются равенства:

Тогда, при выполнении соответствующих условий [1] на оператор

T IH = E = Ex0 {q0 E (x0 |x0 )},

(1.8) (1.9)

сопряженного

E x = (x),
то есть величина

x

может рассматриваться как локальная оценка решения

уравнения переноса (1.2). Дисперсия исходной оценки



связана с ковариационной матрицей

(x)

векторной оценки

x

следующим соотношением:

T V = Ex0 (q0 (x0 )q0 ) - (IH )2 .
Таким образом, из ограниченности ковариационной матрицы оценки следует конечность дисперсии T 1/ 2 пропорциональна (S (x)(x)S (x)) ,

(x)



. Если, кроме того, плотность распределения

(x)

то дисперсия оценки (1.5) очевидно минимальна.

2

Вычисление локальных характеристик поляризованного излучения
I (x )
в точке

Для вычисления вектор-функции Стокса

x = (r , ),

где

r



и

- = - ( a , b , c )

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

3


Рис. 1. Геометрия среды

OX , солнечное излучение падает на атмосферу по направлению ( -1, 0, 0 ) (см. рис. 1). Пусть поверхность, образованная лучом r - t, t 0 при вращении вокруг оси симметрии O X ; r1 = ( x1 , y1 , z1 ) точка пересечения направления с поверхностью . Точку r1 поворачиваем вокруг оси симметрии в азимутальную плоскость наблюдения. При этом точка r1 переходит в точку r2 на луче r - t, а направление в 1 . Модифицированную двойную локальную оценку метода Монте-Карло для вычисления I (x ) запишем в виде оценки по рассеяниям
Будем считать, что система Солнце-атмосфера-Земля симметрична относительно оси следующим образом:

N

N

i (x ) =
n=0

q

T n

H (xn , x )

i

=q

T 0 n=0

Qn H (xn , x ) ,
i

e E i (x ) = Ii (x ) = IHi = e , Hi .
Здесь

H (xn , x ) = s (r



n,1

)exp{- (rn , rn,1 )- (r

n,2

, r )}P (n,1 , , rn,1 )



T



1-n2 1 2 |a |1 | n,n |



(xn )

матрица

4 Ч 4; Hi (x) = (H (x, x ))i столбец матрицы H (x, x ), (xn ) индикатор того, что после рассеяния в точке rn направление n пересекает поверхность ; n = ( n1 , n2 , n3 ) 2 2 нормаль к поверхности в точке r1 ; 1 = y1 + z1 . T Введя вектор (x ) = ( 1 , 2 , 3 , 4 ) (x ) можно записать:
размера

S T (x ) (x ) = (x )
T

N

Qn H (xn , x ) =
n=0

S T (x ) x , (x ) 0

E (x ) = I (x ).

Как известно, дисперсия этой оценки расходится логарифмически. Однако при расчетах величину H (x, x ) заменяют на величину H (x, x ) = H (x, x ) , где индикатор того, что знаменатель в выражении . Для реальной атмосферы с сильно вытянутой индикатрисой это мало меняет результат, а дисперсия оценок ( )i (x ) становится конечной. двойной локальной оценки. В дальнейшем будут в основном рассматриваться производные по сечению аэрозольного рассеяния для вычисления которых использовались следующие алгоритмы [2].

H

больше



Рассмотрим вопросы, связанные с вычислением параметрических производных от модифицированно

4


Пусть в каждом из

n

непересекающихся областей

Di

, интересующие нас коэффициенты

задаются следующим образом:

ћ,i =

ћ0,i , ћ0,i + t,

при при

i = i0 , i = i0 )= i (r)

(2.1)
m i=1 c,i i (r ), а коэффициент будем обозначать индикатор

Тогда коэффициент поглощения в среде будет равен c (r m рассеяния s (r ) = i=1 s,i i (r ). Здесь и далее символом того, что зависит

r Di от t.

. Предполагается также, что моделируемая плотность перехода

p(x , x)

не

Тогда для производной модифицированной двойной локальной оценки по сечению аэрозольного рассеяния в слое

i

0 получаем следующие вычислительные формулы: N

x (x ) =
n=0

(Qn H (xn , x ) + Qn H (xn , x )), - rn |i0 K T (xn-1 , xn ) + K T (xn-1 , xn ) , p(xn-1 , xn ) - rn |i0 + |rn,2 - r |i0 )H (xn , x ) + H (xn , x ).
n-1

(2.2)

Qn = Qn-

1

-|r K T (xn-1 , xn ) + Qn-1 p(xn-1 , xn )
n,1

H (xn , x ) = -(|r
Здесь

|r - r |i0 длина части отрезка (r , r), лежащая в области Di0 , Ln,i0 = n=1 |rk - rk-1 |i0 . k 0 Величина K (x, x ) соответствует величине K (x, x ) с матрицей рассеяния R(ч) = Ra (ч)/s (r ), 0 а величина H (x, x ) величине H (x, x ) с матрицей рассеяния R(ч) = Ra (ч)/s (r1 ). Условия несмещенности и конечности дисперсии оценок ( )i (x ) и ( )i (x ), даются следующим

утверждением [2]:

Утверждение 1.

Пусть выполнены следующие условия: 1. оператор K является положительным и (K) < 1, 2. в некотором интервале - + выполены условия равномерной ограниченности по : K(i)( ) Ck < и H (i)( ) Ch < , i = 0, ..., m. Тогда оценки ()i(x) и и ()i(x) является несмещенными. Если, кроме того, выполнены условия (Kp) < 1, K(pi)(j) < , i, j = 0, 1 где
K(i)(j ) p (x) =
X

(

i ti

то дисперсия данных оценок конечна.

K (x, x ))T (x ) p(x, x )

j tj

K (x, x )

dx

производных двойной локальной оценки. Нетрудно убедиться, что производные оператора

Рассмотрим подробнее вопрос обоснования несмещенности и конечности дисперсии полученных

K

по сечениям поглощения и аэрозольного рассеяния удовлетворяют вышеприведенным условиям. Иначе обстоит дело с функцией H (x, x ). Вследствие того, что H (x, x ) L (X ) указанные / требования теряют смысл. Но для величины H (x, x ) L (X ) все рассуждения справедливы. Таким образом можно говорить лишь о том, что мы имеем формулы дифференцирования для смещенной оценки ,x (x ). Однако, вследствие ограниченности физической среды и свойств оператора

K,

можно заметить, что в некотором интервале параметров дифференцирования

сходимость

I (x ) - I (x )
0
является равномерной по параметру дифференцирования. Это позволяет предположить, что выполнено соотношение:

lim I (x ) = lim I (x ) = I (x ), 0 0
5


на основании которого можно ожидать, что производные смещенной оценки будут сходиться к производным несмещенной оценки при

0.

Обладая возможностями вычисления локальных характеристик излучения, например при помощи обсуждавшейся выше двойной локальной оценки, можно решать задачи оценивания производных различных функций от компонент вектор-функции I (x ). Из подобных функций наиболее используемой является так называемая степень поляризации P (x ), описываемая следующим образом:

P (x ) =

I2 (x ) Ч 100%. I1 (x )

Рассмотрим задачу вычисления степени поляризации и ее производных. Справедлива следующая общая теорема [4]. (i) ( Теорема 1. 1 , ..., si) i = 1, ..., n (i) E k = E k k = 1, ..., s

Пусть , выборка из s-мерного распределения с конечными средними значениями , , и положительно определенной ковариационной s Ч s матрицей {kl }. Пусть также g (x1 , ..., xs ) функция, имеющая первые производные xg = 0 gk , k = 1, ..., s во всех точках некоторой окрестности точки (E 1 , ..., E s ) и gk = gk (E 1 , ..., E s ). 1 1 0 Тогда, если по крайней мере одна из величин gk = 0, то g n n=1 1(i), ..., n n=1 s(i) имеет для i i больших n асимптотически нормальное распределение со средним g(E 1, ..., E s) и дисперсией s 1 00 k,l=1 kl gk gl . n
k

E

(1) Пусть k (x ), k (x ) = Ik (x ),
В этом

лишь потребовать

k , l = 1, 2.

..., k (x ) независимые реализации случайных величин k (x ) таких, что k = 1, 2. Здесь неважно, каким способом получена эта статистика, следует невырожденность ковариационной матрицы {E ((k (x )-Ik (x ))(l (x )-Il (x )))}, случае в качестве приближенной оценки функции g (I1 (x ), I2 (x )) = P (x ) в ? 2 (x ) gn = ? , 1 (x )

(n)

соответствии с теоремой 1 следует взять выражение

(2.3)
n i=1

где

1 ? 1 (x ) = n
следующим образом:

n i=1

1 (x ),

(i)

1 ? 2 (x ) = n

2 (x ).

(i)

При этом дисперсия асимптотического распределения такой случайной величины будет выражаться

Vgn =

1 V2 (x ) - 2P (x )cov (1 (x ), 2 (x )) + P 2 (x )V1 (x ) . 2 n I1 (x )

Для того, чтобы при помощи теоремы 1 получить выражения для оценивания производных степени поляризации, рассмотрим следующую функцию:

g (x1 , x2 , x3 , x4 ) =

x4 x2 x3 - 2. x1 x1

Нетрудно видеть, что, если компоненты Ii (x ) вектор-функции Стокса зависят от некоторого параметра , то g I1 (x ), I2 (x ), I (x ), I2 (x ) является аналитическим выражением для 1 P (x ). В этом случае в качестве оценки g можно использовать величину:

? ? (x ) (x ) gn = ?2 - gn ?1 , 1 (x ) 1 (x )
где

(2.4)
n i=1

1 ? 1 (x ) = n

n i=1

1 (x ),
6

(i)

1 ? 2 (x ) = n

2 (x ).

(i)


При этом дисперсия асимптотического распределения такой оценки будет выражаться следующим образом:

Vgn = V1 + I2 - 2P I1 2 I1 +2P

I2 - 2P I1 2 I1 2

2

+ V2

I1 2 I1

2

+ V1

P I1

2

+ V2

1 + 2 I1

I1 P 1 cov (1 , 2 ) + 2 cov (1 , 1 ) - 2 cov (1 , 2 ) + 2 I1 I1 I1

I1 I I2 cov (2 , 1 ) - 2 1 cov (2 , 2 ) - 2 3 cov (1 , 2 ) (x ) 3 3 I1 I1 I1

Таким образом, мы получили, что условия несмещенности и конечности дисперсии асимптотических распределений оценок степени поляризации и ее параметрических производных (2.3) и (2.4) практически совпадают с условиями несмещенности и конечности дисперсии соответствующих оценок вектор-функции Стокса I (x ) и ее производных, приведенными выше.

3

Оценки производных однократно рассеянного излучения
m (1)

Нетрудно видеть, что для вклада однократно рассеянного излучения справедлива следующая формула:

I

(x ) =
i=1

Ii ,


ti

Ii =
Здесь

ti-1

s (Ai )exp{-0 (r - t ) - (r , r - t )}R(Ai , -a )



S 0 0 0

dt

номера сло?в, через которые проходит луч r - t , t 0; ti-1 и ti точки пересечения границ слоя Ai с вышеуказанным лучом; 0 (r - t ) оптическая длина отрезка, образованного точкой r - t и точкой r (t) пересечения внешней границы атмосферы T лучом r - t + l ( 1, 0, 0 ) , l 0; R(Ai , -a ) матрица рассеяния в Ai -м слое.

Ai , i = 1 , . . . , m

Тогда справедлива следующая формула для вычисления производной вклада однократно рассеянного излучения по коэффициенту поглощения в

i

0 -м слое:

I

(1)

(x ) = c

m i=1

Ii (x ) , c

Ii (x ) =- c

ti ti-1

s (Ai )exp{-0 (r - t ) - (r , r - t )} Ч


Ч [|r0 (t) - (r - t )|i0 + |r - (r - t )|i0 ] R(Ai , -a )

S 0 0 0

dt

Для производной вклада однократно рассеянного излучения по коэффициенту аэрозольного рассеяния в

i

0 -м слое справедливы аналогичные выражения:

I

(1)

(x ) = a

i

Ii (x ) , a

7


Ii (x ) = a

ti ti-1

s (Ai )exp{-0 (r - t ) - (r , r - t )} Ч


Ч - (|r0 (t) - (r - t )|i0 + |r - (r - t )|i0 ) R(Ai , -a ) + i0 (Ai )

Ra (Ai , -a ) s (Ai )

S 0 0 0

dt

Здесь

i0 (Ai )

индикатор того, что

i0 = A

i.

С помощью данных формул, используя какие-либо процедуры численного интегрирования, можно проводить вычисления производных вкладов однократного рассеяния по коэффициентам

c

и

a

.

4

Сравнительный анализ метода зависимых испытаний и метода прямого дифференцирования для вычисления производ по параметрам
В данном случае следует изучить вопрос трудоемкости данного алгоритма в сравнении с

Рассмотрим метод зависимых испытаний в его применении к вычислению производных рассеянного излучения по сечениям взаимодействия со средой. трудоемкостью метода прямого дифференцирования уравнения переноса. Пусть

{x (h)}

множество оценок решения сопряженного уравнения переноса с учетом

поляризации, смоделированных на одних и тех же цепях Маркова для сред, в которых соответствующий (h)- (0) Отметим, что коэффициент взаимодействия равен ћ + h. Введем обозначение: (h) = h сходимость

E (h) - E
h0
прямо следует из свойства

(0) ћ K
уравнения переноса и мажорантного

(K) < 1

интегрального оператора

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

V (h) - V
h0

(0) ћ
T V = Ex0 (q0 (x0 )q0 ) - (IH )2 ковариационных матриц
, задачу следующим

Нетрудно видеть, что, используя известное соотношение сходимости дисперсий можно свести к задаче сходимости образом:

T T V = V(q0 x0 ) = Ex0 (q0



(x0 )q0 ) -

IH (h) - IH (0) h

2

Рассмотрим подробнее условия сходимости


Введем обозначение



- .
h0

(4.1)

(x, h1 , h2 ) =

x

(h1 ),x (h2 ) . Заметим, что данная ковариационная матрица

является решением следующего уравнения:

(x, h1 , h2 ) = H (h1 )( (h2 ))T + (h1 )H T (h2 ) - H (h1 )H T (h2 ) (x) + +
X
или в операторной форме

K T (x, x , h1 )(x , h1 , h2 )K (x, x , h2 ) dx p(x , x)
при выполнении условия

(ћ, h1 , h2 ) = ч(ћ, h1 , h2 )+Kp (h1 , h2 )(ћ, h1 , h2 ),

(Kp (h1 , h2 )) < 1.
8


Нетрудно убедиться, что сходимость (4.1) эквивалентна существованию обобщенной второй 2 (x,h1 ,h2 ) . Таким образом, можно сформулировать следующее смешанной производной h1 h2 h1 =h2 =0 утверждение:
Утверждение 2.

Пусть выполнены следующие условия: 1. оператор K сохраняет мажорантное свойство первой компоненты вектор-функции Стокса, 2. в некотором интервале значений -1 < h1 < 1, -2 < h2 < 2 выполняется неравенство (h (h h (Kp (h1 , h2 )) qp < 1, причем операторы K h ,h ) , K h ,h ) , Kh ( h,h ) равномерно ограничены. Тогда существует сходимость
p 1 2 p 1 2 2 p 1 2 1 2 1 2

V

(h) - (0) - V (0) h0 h ћ

Таким образом, погрешность метода зависимых испытаний стремится к погрешности метода прямого дифференцирования при уменьшении шага

h.

В этом нет ничего удивительного, так

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

1 N

N n=1

n 1 (0) = ћ N C (N ) =

N n (h) + C (N )h + O(h2 ) n=1

Нетрудно заметить, что величина с математическим ожиданием

1 (0) Nn

является асимптотически нормально распределенной

I (0)

и дисперсией

V

2 (V )2 + (V )2 + O(h). hћN

5

Методы решения обратной задачи атмосферной оптики
n
слоев, в каждом из которых

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

=

( 1 , . . . , n ),

где

i ,

i = 1, . . . , n

значения коэффициентов взаимодействия со средой в

слоях разбиения. в

I (r, ) = ( I1 , I2 , I3 , I4 )T (r, ) значение вектор-функции Стокса, наблюдаемое точке r в направлении . Пусть fi (I (r, )), i = 1, . . . , m некоторые функции от компонент вектора I (r, ). Ставится
Обозначим через

задача: по заданным значениям

f1 (I ( ; r1 , 1 )), . . . , fm (I ( ; rm , m )) a a
вычислить значения коэффициентов аэрозольного рассеяния



a в слоях разбиения (

mn

).

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

f1 (I ( a ; r1 , 1 )) = f1 (I ( ; r1 , 1 )), a ... fm (I ( a ; rm , m )) = fm (I ( ; rm , m )) a

(5.1)

9


может быть решена итерационным методом Ньютона-Канторовича. Соответствующая линеаризованная (s ) система для нахождения s-го приближения a записывается, очевидно, в следующем виде:



n

4

i=1 ... j =1 n4 i=1 j =1

f1 Ij ( Ij a,i

(s-1) a

; r1 , 1 )

(s ) a

-

(s-1) a

= f1 (I ( ; r1 , 1 )) - f1 (I ( a

(s-1) a

; r1 , 1 ))
(s-1) a

fm Ij ( Ij a,i

(s-1) a

; rm , m )

(s ) a

-
( s)

(s-1) a

= fm (I ( ; rm , m )) - fm (I ( a

; rm , m )) (5.2)

или, в матричном виде,

D a = F

( s) (s )

(s )

. Здесь

(s) a =


(s ) a

-
. . .

(s-1) a

,
(s-1) a

f1 (I ( ; r1 , 1 )) - f1 (I ( a

; r1 , 1 )) ; rm , m ))



(s ) F

(s )

=

4 k=1 fi Ik Ik a,j

fm (I ( ; rm , m )) - fm (I ( a
Матрица производных

(s-1) a

j =1,...,n

является при m > n i=1,...,m прямоугольной, вследствие чего система (5.2) оказывается переопределенной. В этом случае для решения системы (5.2) применяется метод наименьших квадратов, приводящий, как известно, к необходимости решения следующей задачи:

D = {dij } =

(

(s-1) a

; rm , m )

D ( s) (s ) a - (s ) F

(s )

, D (s ) (s ) a - ( s ) F

(s )

- min

(s)



a

Отсюда получается квадратная система линейных уравнений относительно

( s)

a:

D

(s ) T

D (s ) (s ) a = D

(s ) T

( s) F

(s)

(5.3) Ax = b
с

Как известно, относительная погрешность решения системы линейных уравнений возмущенной матрицей формулой [5]:

~ A = A + A

и вектором правых частей

~ = b + b b ,

задается следующей

x (A) x 1 - (A) ћ
где

A A

b A + b A

(5.4)

(A) = A ћ A

-1

. При этом оценка

(5.4)

чем больше число обусловленности матрицы Выражение

является неулучшаемой. Отсюда видно, что, T D(s) D(s) (s) , тем точнее следует вычислять

элементы рассматриваемой системы линейных уравнений.

(5.4)

позволяет исследовать сходимость решения

x ~

системы линейных уравнений

~~ b Ax = ~,

элементы которой рассчитаны методом Монте-Карло, к решению

уравнений

Ax = b

, где

~ A = E A, b = E~. b x > , x x > , x

Действительно, для погрешности

x системы x = x - x ~
-1

линейных получаем:

P

x > =P x P

A < A A < A

-1

(A) + P

x > , x A > A (5.4)
-1

A > A (A)

(A)

-1

(A) + P

Рассмотрим подробнее второе слагаемое. С учетом формулы

имеем:

P

x > , x

A < A

-1

(A) P



(A)
A A

1 - (A)

A b + A b

> ,

A < A

-1

(A)



10


Так как функция зависящая от

(A c f (x) = 1-)(x+x A) A и b константа k

возрастает при любом положительном такая, что выполнены оценки:

c

, то существует не

1 - (A) (A) =P 1 - (A)

P



(A)
A A

A b + A b A b + A b

> ,

A < A

-1

(A)



=


A A

A > , < < k A > ,

-1

(A) +

+P 1 - (A) P A > A k P
Будем считать, что

(A)
A A

A b + A b

A < A k A < A k



+P

b k - (A)(1 + ) > , b k (A) +P

A > A k

b k - (A)(1 + ) > b k (A)
получаем:

< 1.

Тогда, выбирая

k = 3 (A),

P

x > P x

A > +P A 3 (A) 2P

b > +P b 3 (A) b > b 3 (A)

A > A

-1

(A)

A > +P A 3 (A)

Наши дальнейшие рассуждения будут зависеть от выбора векторной нормы с ней матричной нормы

b

и согласованной

. Возьмем, например, в качестве векторной нормы b = maxi |bi |, а n в качестве матричной нормы A = max1j n i=1 |aij |. Тогда, с учетом очевидного соотношения 2 A n maxi,j |aij |, получаем:

A

P max | xi | > x
i

2P max | aij | >
i,j

A b + P max | bi | > i 3n (A) 3 (A)
2

Это выражение позволяет нам утверждать, что порядок сходимости по вероятности решения системы линейных уравнений

~~ b Ax = ~

не меньше, чем порядок сходимости элементов матрицы

~ A

и вектора

~ b

к их математическим ожиданиям.

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

m

F ( a ) =
i=1

(fi (I ( ; ri , i )) - fi (I ( a ; ri , i )))2 - min a

a

для нахождения минимума функционала

F ( a )

использовался метод покоординатного спуска.

Этот метод имеет следующие преимущества. Для вычисления нескольких первых итераций можно использовать гораздо меньшее число траекторий и учитывать только знаки получающихся производных; другими словами, можно воспользоваться следующим итерационным процессом:

(n) a,i

F ( (n-1) ) (n-1) a = 1 - sgn a,i = 1 - sgn Fa,i ( (n-1) Fa,i ( a )







(n-1) a

)

a,

(n-1) i

(5.5)

При этом время счета можно существенно сократить.

11


6

Численные результаты

В качестве тестовой задачи рассмотрим модель земной атмосферы, в которой распределение по слоям коэффициентов поглощения, молекулярного и аэрозольного рассеяния задается таблицей 1. Альбедо подстилающей поверхности (Земли) равно нулю. Матрица аэрозольного рассеяния рассчитана по теории Ми на основе данных о микрофизических параметрах атмосферного аэрозоля.
Таблица 1. Коэффициенты взаимодействия со средой

Слой 02 25 5 11 11 25 25 38 38 55 55 70 70 80 80 100 100 120 120 150 150 500

m 1.106 0.906 0.594 0.258 0.264 0.384 0.248 0.430 0.743 0.176 0.998 0.254 Ч Ч Ч Ч Ч Ч Ч Ч Ч Ч Ч Ч 10-2 10-2 10-2 10-2 10-3 10-4 10-5 10-6 10-7 10-8 10-10 10-11

a 0.283 0.357 0.306 0.483 0.128 0.251 0.267 0.815 0.558 0.417 Ч Ч Ч Ч Ч Ч Ч Ч Ч Ч 0 0 10-3 10-3 10-3 10-3 10-2 10-3 10-5 10-7 10-8 10-10 0.456 0.695 0.474 0.131 0.187 0.387 0.300 0.518 0.964 0.219 0.107 0.697

c Ч Ч Ч Ч Ч Ч Ч Ч Ч Ч Ч Ч 10-1 10-2 10-3 10-3 10-4 10-5 10-6 10-7 10-8 10-9 10-10 10-12

Наблюдения проводятся с высоты 2.