Документ взят из кэша поисковой машины. Адрес оригинального документа : 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.7 км в направлении в зенит при различных зенитных расстояниях Солнца. На рис. 23 приведены рассчитанные степени поляризации

P

и компоненты

I1 , I2

вектора Стокса. Рис. 47 иллюстрируют зависимость производных степени поляризации

по сечениям аэрозольного рассеяния от зенитных расстояний Солнца. Пунктиром показаны соответствующие характеристики однократно рассеянного излучения. Вычисления проводились на 1000000 траекториях. Погрешности соответствующих рассчитанных величин представлены в таблице 2:
Таблица 2. Погрешности производных степени поляризации.

Зенитн. расст. Солнца,

VP ,% (E P )2
0.19 0.17 0.24 0.34 0.47 0.91 1.32 1.58 1.67 2.00 3.66 4.79

VP (E P

a,2

a,2

)

2

,%

VP (E P

a,5

a,5

)

2

,%

90.5 91 92 93 94 95 96 97 98 99 100 101

1.25 1.34 1.48 1.97 1.71 6.73 14.02 10.73 9.26 12.63 42.40 12.51

1.17 1.02 1.27 1.58 2.40 5.22 5.06 9.16 12.38 21.21 34.49 70.08

12


Рис. 2 Степень поляризации

P

Рис. 3 Компоненты

ln(107 I1 )

(I) и

ln(107 I2 )

(I I)

Рис. 4 Производная

P

во втором слое

Рис. 5 Производная

I2

во втором слое

Рис. 6 Производная

P

в пятом слое

Рис. 7 Производная

I2

в пятом слое

13


D fi (I ( a ; ri , i )), i = 1, . . . , m
1. 2. 3. 4.

(s )

В Таблице 3 приводятся числа обусловленности (A) = max (A)/min (A) для матриц A = ~ T (s ) D , построенных в соответствии с п. 5 при помощи следующих наборов функций : ; ; ;

fi = I1 (ri , i ), i = 1, . . . , m fi = I2 (ri , i ), i = 1, . . . , m f2i = I1 (ri , i ), f2
i+1

= I2 (ri , i ), i = 1, . . . , m
.

fi = P (ri , i ), i = 1, . . . , m

Для сравнения приводятся также числа обусловленности матриц, построенных при помощи первого набора функций для излучения без учета поляризации (т. е. для уравнения вида (1.1) с ядром

Kij (x , x) = 0

при

i, j = 1). D
(s ) T

Таблица 3. Числа обусловленности матриц

D

(s)
Неполяриз. излучение

Зенит. расст. Солнца

Поляриз. излучение

I1 [90.5 , 91.7 ], [90.5 , 92.3 ], [90.5 , 95.5 ], [90.5 , 98.5 ], [90.5 , 100.5 ],
шаг шаг шаг шаг шаг

I2 106 106 106 106 107 2.1 1.4 4.1 8.4 1.0 Ч Ч Ч Ч Ч 105 105 105 105 106 1. 1. 3. 6. 8.

I1 , I2 0 2 3 9 9 Ч Ч Ч Ч Ч 106 106 106 106 106

I2 /I1 1.1 4.6 4.8 1.4 7.4 Ч Ч Ч Ч Ч 1015 1014 1016 1021 1021 4.0 4.5 1.2 2.1 2.7

I1 Ч Ч Ч Ч Ч 105 105 106 106 106

0.1 0.1 0.1 0.1 0.1

1. 1. 4. 8. 1.

8 5 6 6 1

Ч Ч Ч Ч Ч

Как известно, решение систем линейных уравнений с подобными плохо обусловленными матрицами сопряжено с определенными трудностями. В частности, такие задачи являются чувствительными к точности задания входных данных. В нашем случае ситуация усугубляется наличием вероятностных погрешностей компонент рассчитываемой матрицы решаемой системы. Очевидно, приведенные в Таблице 3 значения

(A) ~

являются оценками снизу для величин

(A).

Таким образом, видно, что при решении обратных задач необходимо весьма точно рассчитывать

производные характеристик поляризованного излучения. Однако и в этом случае при проведении численных экспериментов выяснилось, что алгоритм восстановления высотного хода коэффициента аэрозольного рассеяния в атмосфере, разбитой на большое число слоев, расходится. В Таблице 4 приведены значения коэффициентов аэрозольного рассеяния, вычисленные на 100000 траекториях методом Ньютона-Канторовича с использованием функций

fi = I2 (ri , i ).

Атмосфера толщины 100 км была разбита на пять сло?в толщиной 10,

5, 15, 40 и 20 км. Время счета на компьютере Pentium I I I 665 МГц составило 364 с.
Таблица 4. Значения коэффициентов

a

, рассчитанных методом Ньютона-Канторовича

Приближение Начальное 1-е Точное значение

3.962 6.768 3.741 1.141 7.811

Ч Ч Ч Ч Ч

10- 10- 10- 10- 10-

4 4 5 5 7

1. 4. 2. 9. 3.

981 725 496 564 905

Ч Ч Ч Ч Ч

10- 10- 10- 10- 10-

4 4 5 6 7

2.830 4.834 2.672 8.153 5.579

Ч Ч Ч Ч Ч

10- 10- 10- 10- 10-

4 4 5 6 7

Попытки рассчитать рассматриваемым методом параметры атмосферы, разбитой на число слоев большее или равное 6, с использованием любого из вышеприведенных наборов функций

fi

14


не увенчались успехом. Таким образом, становится понятно, что метод Ньютона-Канторовича не подходит для решения обратных задач атмосферной оптики для атмосферы, разбитой на большое число слоев. Итерационный процесс (5.5) в данном случае обнаруживает следующие преимущества. Вопервых, для его применения не требуется большая точность в расчетах производных поляризованного излучения. Во-вторых, результаты, приведенные в п. 3, показывают, что, по крайней мере для зенитных расстояний Солнца в интервале 90.5 95.5 , для определения знака производных можно пользоваться производными однократно рассеянного излучения, что также уменьшает время счета алгоритма. В качестве примера в Таблице 5 приведены значения коэффициентов аэрозольного рассеяния, вычисленные при помощи итерационного процесса (5.5) на 10000 траекториях. Коэффициент



был взят равным

0.03.

Время счета на компьютере Pentium I I I 665 МГц составило 1935 с.

Таблица 5.Значения коэффициентов

a

, рассчитанных методом (5.5)

Приближение Начальное 1-е 10-е Точное значение

3.962 5.001 4.291 6.768 1.788 3.507 2.672 1.141 7.811 5.845

Ч Ч Ч Ч Ч Ч Ч Ч Ч Ч

10 10-4 10-4 10-4 10-3 10-4 10-6 10-7 10-9 10-11

-4

3.843 4.851 4.162 6.565 1.734 3.402 3.629 1.107 7.576 5.670

Ч Ч Ч Ч Ч Ч Ч Ч Ч Ч

10 10-4 10-4 10-4 10-3 10-4 10-6 10-7 10-9 10-11

-4

2.922 3.688 3.164 4.991 1.400 2.586 3.110 1.008 6.494 4.577

Ч Ч Ч Ч Ч Ч Ч Ч Ч Ч

10 10-4 10-4 10-4 10-3 10-4 10-6 10-7 10-9 10-11

-4

2.830 3.572 3.065 4.834 1.277 2.505 2.672 8.153 5.579 4.175

Ч Ч Ч Ч Ч Ч Ч Ч Ч Ч

10-4 10-4 10-4 10-4 10-3 10-4 10-6 10-8 10-9 10-11

Время счета, сек.:

197

1935

7

Литература
1. G. I. Marchuk, G. A. Mikhailov, M. A. Nazaraliev, et al., The Monte Carlo Method in Atmospheric Optics. Springer Verlag, Berlin; Heidelb erg; New York, 1980

2. S. A. Ukhinov, D. I. Yurkov, Monte Carlo metho d of calculating the derivatives of p olarized radiation. Russian Journal of Numerical Analysis and Mathematical Model ling, 1998, Vol. 13, No. 5, 425-444. 3. G. A. Mikhailov, Minimization of Computational Costs of Non-analogue Monte Carlo Methods. World Scientic, Singap ore; New Jersey; London; Hong Kong, 1991 4. С. Уилкс, Математическая статистика. М., Наука, 1967 5. А. Н. Коновалов, Вычислительные методы линейной алгебры. Новосибирск, Наука, 1989

15