Помимо собственно расчета потоков солнечной радиации, для решения обратной задачи, методами, изложенными в гл. 4., необходим и расчет производных от них по всем восстанавливаемым параметрам.

Вычисление производных в используемом методе Монте-Карло основано на дифференцировании формального ряда Неймана (2.1.22) [9], который есть общий вид решения задачи переноса излучения [7,9]. Чтобы не повторяться, договоримся сохранять ниже все обозначения, что и в разделе 2.1. Пусть требуется найти производную от решения прямой задачи (2.1.22)

B = q + Kq + K²q + K³q + …

по некоторому параметру a, т.е. величину , где индекс a у интегральных операторов символически обозначает их зависимость от параметра. Продифференцируем ряд Неймана, выражаемый формулой (2.1.22):

(5.2.1)

Используя для краткости обозначения производных запись вида :

(5.2.2)

и учитывая символическую запись интегральных операторов (2.1.20), получаем для производной левой части ряда (5.2.1)

(5.2.3)

где

.

(5.2.4)

Формула для производной преобразована к виду (5.2.3), (5.2.4) специально. Записанная таким образом, она имеет вид интеграла (2.1.20), который непосредственно вычисляется по методу Монте‑Карло, согласно выражению (2.1.21):

.

(5.2.5)

То есть, по соотношению (5.2.5) вычисление производных в методе Монте‑Карло сводится просто к умножению записываемой в счетчик величины на некоторую дополнительную “весовую” функцию W_a() [9].

Чтобы построить конкретный алгоритм вычисления W_a(), найдем явный вид производных правой части ряда (5.2.1), для чего воспользуемся известным выражением производной произведения через сумму логарифмических производных (xyz…)'=(xyz…)(x'/x + y'/y + z'/z +...). Получаем:

(5.2.6)

Записывая выражения (5.2.6) в виде (5.2.5), более удобном для осуществления метода Монте-Карло, окончательно имеем

(5.2.7)

Из (5.2.7) следует, что производные в методе Монте-Карло можно вычислять при помощи тех же алгоритмов, что и сами искомые величины, умножая при каждой записи в счетчик величину () на специальный “вес” , причем если () зависит только от текущего значения случайной величины , т.е. текущих координат фотона, то  – сумма зависит от всей “истории” его траектории.

Возвращаясь к задаче расчета потоков, получаем, что для вычисления их производных достаточно продифференцировать по восстанавливаемым параметрам явные выражения для функций _a(u), q_a(u) и K_a(u,u'). После чего в описанный в разделе 2.1 вычислительный алгоритм расчета потоков вносятся элементарные изменения: подсчет при каждом моделировании элемента траектории фотона величин Wa (для всего набора параметров) и запись в специальные счетчики производных одновременно с записью в счетчики потоков. Хотя в нашем случае потоки вычисляются как интегралы по длине волны (5.1.4), но при моделировании каждой отдельной траектории длина волны остается фиксированной. Следовательно, при дифференцировании достаточно ограничиться рассмотрением монохроматического случая, производная же от интеграла (5.1.4) при этом получится “автоматически”. Также следует учесть, что оптическая толщина слоя сама является функцией дифференцируемых параметров. Поэтому при вычислении производных в качестве вертикальной координаты следует использовать атмосферное давление. При реальном моделировании от этого ничего не изменится, но при выводе (2.1.8) вероятность свободного пробега фотона от уровня высоты по давлению P1 до уровня P теперь запишется через явное выражение для оптической толщины слоя как :

где (P) коэффициент ослабления, соответственно плотность вероятности (2.1.8) примет вид :

(5.2.8)

Именно формулу (5.2.8) следует использовать в качестве плотности вероятности свободного пробега фотона при дифференцировании.

Теперь, конкретно для алгоритма вычисления величин потоков, описанного в разделе 2.1, учитывая приведенный там же явный вид функций, входящих в выражения (5.2.7), получаем следующий алгоритм вычисления производных.

Вводим счетчики W_a для всего набора параметров. В начале каждой траектории фотона присваиваем W_a:= 0. При моделировании каждого свободного пробега фотона с учетом выражения (5.2.8) присваиваем счетчику значение:

(5.2.9)

где '(P₁,P₂) – длина свободного пробега фотона от уровня P₁ до уровня P₂ (2.1.7). Если фотон попал на подстилающую поверхность, то слагаемое c величиной (P₂), очевидно, будет отсутствовать. При моделировании каждого акта взаимодействия фотона с атмосферой, т.е. при умножении веса фотона на ₀('), присваиваем счетчику значение :

(5.2.10)

где P' – текущая координата фотона (по атмосферному давлению), соответствующая оптической глубине '. Аналогично для взаимодействия фотона с подстилающей поверхностью в соответствии с выражением (2.1.23) счетчику присваиваем значение:

(5.2.11)

На каждом шаге моделирования процесса рассеяния фотона в атмосфере согласно выражению (2.1.9) присваиваем счетчику значение:

(5.2.12)

Наконец, при каждой записи в счетчик потоков веса по формуле (2.1.18) одновременно в счетчики производных записываются величины :

,

где P – координата счетчика.

Полученный алгоритм существенно упрощается, если учесть, что фактически в точке моделирования рассеяния P'=P₂ одновременно вычисляется сумма

(5.2.13)

которая, при подстановке выражений для входящих в нее оптических параметров через молекулярную и аэрозольную составляющие (1.2.12), (1.2.13), с помощью элементарных преобразований приводится к виду :

(5.2.14)

где _m, x_m, _a, x_a – объемные коэффициенты и индикатрисы молекулярного и аэрозольного рассеяния. При этом напомним, что индикатриса молекулярного рассеяния, определяемая по формуле (1.2.14) от оптических параметров не зависит. Окончательно, в алгоритме при моделировании длины свободного пробега фотона счетчику присваивается только значение:

(5.2.15)

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

(5.2.16)

Далее приведем явные выражения для указанных выше производных по восстанавливаемым параметрам обратной задачи. Общий набор восстанавливаемых параметров был определен в предыдущем разделе. Это – вертикальные профили температуры воздуха T(P_i) и профили концентрации четырех поглощающих излучение газов Q_H2O(P_i), Q_O3(P_i), Q_NO2(P_i),  Q_NO3(P_i) (концентрация O₂ постоянна), объемные коэффициенты аэрозольного рассеяния и поглощения _a(P_i,_j), _a(P_i,_j) и альбедо поверхности A(_j). Концентрации атмосферных газов договоримся выражать в объемном отношении смеси, что дает по формуле (1.2.15) для их счетных концентраций простое соотношение

(5.2.17)

Наборы высотных уровней P_i и длин волн _j пока будем считать заданными в общем виде, их конкретные значения будут получены на основе анализа производных (см. раздел 4.4). Заметим также, что на практике для упрощения вычислений производных (и избежания ошибок в процессе программирования) следует записывать их в виде цепочки простейших формул, максимально используя правило дифференцирования сложной функции. Причем это имеет смысл делать, даже если при подстановке производных в общие выражения происходят упрощения формул. Другим приемом, эффективно упрощающим вычисления, является применение выражения производной произведения через логарифмические производные (см выше).

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

найдем производную от функции

Определив номер n из условия u_n u u_n+1, имеем равенства:

 если i < n или i > n+1;

 если i = n,

 если i = n+1.

Поскольку производная зависит только от аргумента, введем для нее обозначение

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

Длина свободного пробега фотона '(P₁P₂) по выражениям (2.1.1)-(2.1.4) является квадратичной функцией объемного коэффициента ослабления (P_i). Отсюда для вычисления производной

,

где для определенности положим P₂ < P₁, получаем следующий алгоритм:
1) Находим номера n₁ и n₂ из условий P_n P₁ P_n+1, P_n P₂ P_n+1.
2) Далее, в зависимости от величины разности n₂ – n₁ рассматриваем три случая.
n₂ > n₁ + 1

В этом случае производная , если i < n₁ или i > n₂+1;

 если i = n1;

  если i = n₂+1;

, если n1+2 i n2-1;

(5.2.18)

 если i = n₂;

 если i = n₂+1.

n₂ = n₁ + 1. Этот случай отличается от предыдущего тем, что при i = n₁ + 1 = n₂ производная равна

(5.2.19)

n₂ = n₁. При этом условии производная равна:

, если i < n1 или i > n1 + 1; , если i = n1 = n2, , если i = n1+1=n2+1.

(5.2.20)

Учтем, что объемный коэффициент ослабления при использовании в описанном алгоритме вычисляется на единицу давления aP(Pi), хотя как исходный вычисляется на единицу высоты az(Pi). Тогда, дифференцируя соотношение (5.1.1), находим

(5.2.21)

В силу “правил сложения” (1.2.13) аналогичные соотношения получатся и для объемного коэффициента аэрозольного рассеяния.

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

Производные по концентрациям поглощающих излучение атмосферных газов, (кроме водяного пара). От этих концентраций зависит объемный коэффициент молекулярного поглощения _m(P_i), а от него, в свою очередь, по формуле (1.2.13) объемный коэффициент ослабления. Тогда, обозначая индексом k конкретный газ, получаем:

(5.2.22)

где:

  согласно выражениям (1.2.12) и (1.2.13),

 согласно (1.2.11) и

 согласно (5.2.9).

Сечения Сa молекулярного поглощения в зависимости от длины волны и (для озона) температуры вычисляются линейной интерполяцией по формулам (1.2.17), (5.1.7). Разумеется, производные по концентрациям газов не равны нулю только в спектральных диапазонах учета поглощения этих газов (табл.5.1).

Производная по концентрации водяного пара. От концентрации H₂O помимо объемного коэффициента молекулярного поглощения зависит по формуле (1.2.16) объемный коэффициент молекулярного рассеяния. Это дает для производной от длины свободного пробега

(5.2.23)

и для производной от объемного коэффициента молекулярного рассеяния

.

(5.2.24)

Производные, связанные с объемным коэффициентом поглощения, вычислены выше, сечение поглощения для H₂O вычисляется по формуле (5.1.8).

Для производной от объемного коэффициента молекулярного рассеяния согласно (1.2.16) имеем :

,

(5.2.25)

где

.

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

,

(5.2.26)

где  с учетом соотношений (1.2.12), (1.2.13).

Производная по объемному коэффициенту аэрозольного рассеяния. От этого параметра зависят объемные коэффициенты поглощения, рассеяния и по (5.1.9) индикатриса аэрозольного рассеяния. Получаем :

,

(5.2.27)

где  с учетом соотношений (1.2.12), (1.2.13).

Далее получаем:

.

(5.2.28)

Наконец, для производной от индикатрисы рассеяния справедливы соотношения

(5.2.29)

и с учетом формулы (5.1.9), после несложных преобразований получаем:

,

(5.2.30)

где

 .

Производная по температуре воздуха. От температуры зависят достаточно много величин. Начнем с длины свободного пробега фотона. Для нее получаем:

(5.2.31)

и для объемного коэффициента молекулярного рассеяния

.

(5.2.32)

Важной особенностью расчета производных по температуре является необходимость учета температурной зависимости в формуле пересчета объемных коэффициентов ослабления в терминах атмосферного давления (5.1.1). Получаем

.

(5.2.33)

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

.

(5.2.34)

В конечном счете, задача сводится к дифференцированию объемных коэффициентов молекулярного рассеяния и поглощения. Первый коэффициент по формуле (1.2.11) равен сумме коэффициентов поглощающих излучение газов (всех, включая O₂). Соответствующая сумма получится и для производных. Обозначая конкретный газ индексом k, имеем с учетом (5.2.9):

.

(5.2.35)

Используемые сечения поглощения газов NO₂, NO₃ и O₃ в диапазоне 426‑848 нм от температуры не зависят, для них справедливо равенство . Для O₃ в диапазоне 330‑356 нм из (5.1.7) получаем

.

(5.2.36)

Для O₂ и H₂O соотношение (5.1.8), с учетом линейной интерполяции сечений по длинам волн, дает выражение:

.

(5.2.37)

Для производной от объемного коэффициента молекулярного рассеяния из соотношения (1.2.14) с учетом (1.2.15) получаем :

,

(5.2.38)

а формула (1.2.16) дает следующее выражение:

(5.2.39)

По результатам расчета производных, методами, описанными в разделе 4.4, были выбраны конкретные наборы высот и длин волн для восстанавливаемых параметров атмосферы, а именно: Сетка по длинам волн: от 325 до 685 нм с шагом 20 нм и от 725 до 985 нм с шагом 40 нм (всего 28 точек). Сетка по высоте: от 1000 мбар до 800 мбар с шагом 10 мбар, от 800 – 500 мбар через 20 мбар, 500 – 110 мбар через 30 мбар, 90 – 10 мбар с шагом 10 мбар и уровни 5, 2 и 0,5 мбар (всего 61 точка).

Выбор столь подробной сетки в нижних слоях атмосферы связан с тем, что сами измерения потоков имеют шаг 100 мбар по высоте. Заметим, что высота верхней границы атмосферы – 0,5 мбар (примерно 55 км) хорошо согласуется со стандартно используемой в расчетах переноса излучения в коротковолновой области спектра верхней границей атмосферы 60 км [3,4].

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

Зависимость величины восходящего потока от альбедо подстилающей поверхности хорошо изучена [30,31]. Она имеет вид неоднородной линейной функции (y = ax + b), где мультипликативный член есть прямо пропорциональная альбедо доля отраженного от поверхности нисходящего потока, а аддитивный член связан с рассеянием излучения в атмосфере. Соответственно, чем больше альбедо, тем сильнее зависит от него восходящий поток, и тем менее он информативен относительно атмосферных параметров. Зависимость от альбедо поверхности проявляется также и в нисходящих потоках (см. раздел 3.4). Соответствующая производная тем больше, чем больше альбедо и сильнее рассеяние в атмосфере, и, как следует из расчетов, для светлых поверхностей типа снега величина производной может достигать нескольких десятых процента вариации потока на один процент вариации альбедо. То есть влияние альбедо поверхности на нисходящий поток вполне может превышать погрешность измерений величины этого потока.

Вне полос поглощения O2 и H2O зависимость величин потоков от температуры проявляется очень слабо, даже при максимальных априорных вариациях температуры вариации потоков оказываются на уровне погрешности измерений. Это относится и к полосам поглощения O₃. Таким образом, вне полос поглощения кислорода и водяного пара температурная зависимость потоков может не учитываться, и соответствующие производные могут считаться равными нулю. Вместе с тем в полосах поглощения O₂ и H₂O, включая даже слабые полосы, температурная зависимость существенна, а в некоторых участках спектра, например, на длине волны 932 нм в центре полосы H2O – сильная, достигающая процента вариации потока при вариации всего профиля температуры на один градус.

Производные по концентрации водяного пара также существенны только в полосах его поглощения, т.е. зависимостью объемного коэффициента молекулярного рассеяния от концентрации H₂O можно пренебрегать. Значения указанных производных максимальны в полосе поглощения 910 – 980нм, где вариация потока для априорных вариаций всего вертикального профиля концентрации H₂O достигает 40%.

Максимум величин производных по концентрации O₃ находится в стратосферном озоновом слое. Заметим, что выбор верхней границы атмосферы для расчетов (0,5 мбар) обусловлен именно влиянием стратосферного озона на величину потоков солнечной радиации, поскольку влияние всех прочих компонентов (включая аэрозоли), как показывает анализ производных, на больших высотах пренебрежимо мало. Максимальная вариация величины потока для диапазона априорных вариаций озона на длине волны 330 нм составляет около 5%.

Значение производных по концентрациям N₂O очень малы и, даже учитывая возможный широкий диапазон априорных вариаций его концентрации, вряд ли позволяет надеяться на восстановление концентрации N₂O. Этот вывод не противоречит результатам, полученным в предыдущем разделе, поскольку там мы намеренно использовали экстремально высокие значения концентрации поглощающих газов, а производные рассчитывали для средней модели [3,4]. Аналогичная картина получается и для NO₃, хотя производные по его концентрации в максимумах поглощения (524 и 662 нм) существенно больше и, в принципе, при высоких концентрациях NO₃ позволяют получать определенную информацию о содержании этого газа.

Достаточно сложную высотную зависимость имеют производные по объемному коэффициенту аэрозольного рассеяния. Величина объемного коэффициента аэрозольного рассеяния влияет на потоки солнечной радиации за счет двух противоположных процессов: уменьшения потоков за счет роста оптической толщины слоя атмосферы и их увеличения за счет рассеяния. Поэтому профили указанных производных знакопеременны: вблизи точки измерения они имеют положительный максимум, который уменьшается по мере удаления от этой точки и затем значения производной переходят в отрицательную область. По-видимому, это связано с более “локальным” характером влияния рассеяния излучения: оно дает максимальный вклад в потоки непосредственно вблизи точки их измерения. По абсолютной величине производные от объемного коэффициента аэрозольного рассеяния весьма велики: вариации коэффициента даже в отдельных слоях могут вызывать вариации потока до 10% и более. Спектральный ход указанных производных выражен слабо. Существует методика восстановления высотного хода аэрозольных параметров по дистанционным измерениям радиационных характеристик в полосе поглощения кислорода 760 нм [32,33]. Действительно, некоторое различие вида вертикальных профилей производных внутри и вне этой полосы имеется, однако оно весьма слабое, что собственно подтверждают и выводы, сделанные в работе [33]. Впрочем, для рассматриваемых результатов измерений потоков вертикальный ход восстанавливаемых параметров “автоматически” получается из самолетных измерений на разных уровнях атмосферы.

Величина производных по объемному коэффициенту аэрозольного поглощения существенно зависит от типа выбранной аэрозольной модели. Значения этих производных тем больше, чем больше собственно аэрозольное поглощение. Это явилось причиной того, что по данным измерений над Ладожским озером восстановление объемных коэффициентов аэрозольного поглощения оказалось проблематичным, а для измерений над пустыней – вполне возможным. Этот вывод, собственно, следует уже из проведенного в разделе 3.3 анализа спектральных притоков солнечного излучения.