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

где A и A₀ – заданные матрица и вектор. Заметим, что все полученные в предыдущем разделе соотношения МНК в конечном счете имеют именно такой вид, хотя здесь вектор X будем трактовать в более широком смысле, как любую линейно зависящую от Y величину.

Данные экспериментальных измерений Y содержат случайную погрешность, характеризуемую значением СКО компонентов y_i, i = 1,…,N. В более общем случае, хотя на практике этого стараются всячески избегать, ошибки измерений могут быть коррелированные, т.е. взаимосвязанные. Таким образом, в общем случае ошибки измерений описываются симметричной ковариационной матрицей S_Y, размерности N x N, определение которой, следуя [23], удобно схематически записать как

,

(4.3.2)

где – точное (неизвестное) значение измеряемого вектора, Y – измеренное значение вектора (отличающееся от точного вследствие наличия ошибок измерений), суммирование понимается как усреднение по всем статистическим реализациям измерений случайного вектора (по генеральной совокупности).

Для ковариационной матрицы ошибок S_X параметров X, размерности K x K, очевидно, следует записать соотношение, аналогичное (4.3.2). Но тогда, подставляя в него (4.3.1), получаем

(4.3.3)

Из соотношения (4.3.3) непосредственно вытекает ряд важных следствий.

Следствие 1. Формула (4.3.3) выражает связь между ковариационными матрицами ошибок измерений Y и линейно связанных с ними по (4.3.1) параметров X, т.е. позволяет по известным погрешностям измерений находить погрешности вычисляемых параметров. А именно: значения  есть СКО параметров x_k, значения есть коэффициенты корреляций между погрешностями параметров x_k и x_j. В частном, но часто встречающемся на практике случае некоррелированных ошибок измерений выражение (4.3.3) переходит в удобную для вычислений явную формулу:

k = 1,…,K , j = 1,…,K ,

(4.3.4)

где a_ki – элементы матрицы A, s_i – СКО параметра y_i. В случае равноточных измерений, т.е. s = s₁ = … = s_N, из (4.3.4) следует прямая пропорциональность СКО измерений и параметров:

Следствие 2. Из вывода формулы (4.3.3) очевидно, что генеральная совокупность может быть заменена в ней конечной выборкой из M измерений Y(m), m = 1,…,M, т.е. S_Y в (4.3.2) получена как оценка ковариационной матрицы по известным формулам

 i =1,…,N, j = 1,…,N .

Тогда по формуле (4.3.3) получаются аналогичные оценки и для матрицы S_X. Если речь идет именно о случайных ошибках измерений, то предполагается, что все M измерений относятся к одному реальному значению измеряемой величины. Но, с другой стороны, элементы матрицы S_Y можно трактовать и шире, как характеристику вариаций компонентов вектора Y, вызванных не только случайными погрешностями, но и любыми изменениями измеряемой величины. В таком случае (4.3.3) есть оценка вариаций параметров X по известным вариациям величин Y.

Следствие 3. Рассмотрим простейший случай соотношений вида (4.3.1) – вычисление среднего по всем компонентам вектора Y, т.е.  (здесь K = 1, поэтому X обозначен как скаляр). Тогда a_ki=1/N для всех номеров i и из (4.3.3) для СКО x получаем

.

(4.3.5)

Для некоррелированных ошибок измерений в сумме (4.3.5) остаются только диагональные члены матрицы, и она переходит в хорошо известное правило “суммирования дисперсии”

(4.3.6)

Из (4.3.6) следует, что СКО среднего уменьшается с ростом числа усредняемых значений как  (для равноточного случая s(x)=s(y)/). Поскольку под S_Y можно понимать не обязательно погрешности прямых измерений, свойства (4.3.5) и (4.3.6) часто используются при интерпретации решений обратных задач атмосферной оптики. Например, переход после решения обратной задачи от оптических характеристик тонких слоев атмосферы к достаточно толстым слоям или ко всему столбу атмосферы существенно уменьшает погрешность таких результатов [24]. Заметим также, что мы использовали соотношение (4.3.6) в разделе 2.1 при выводе выражения для дисперсии потоков (2.1.17) в методе Монте-Карло.

Следствие 4. Говоря о соотношении (4.3.6) нельзя не упомянуть еще одно обстоятельство. Оно написано для вещественных чисел, но реально любые представления результатов измерений носят дискретный характер, т.е., в конечном счете, соответствуют целым числам. Дискретность проявляется в наличии некоторой погрешности самого процесса снятия показаний с прибора. Она приводит к тому, что реальная дисперсия s(x) не может быть неограниченно уменьшена, даже если N -> (действительно, измеряя длину линейкой с миллиметровой шкалой, мы, даже проделав миллион измерений, очевидно не получим значение длины с точностью в 1 мкм, а ведь из (4.3.6) следует, что должно быть именно так). К сожалению, в литературе достаточно мало внимания уделяется вопросу влияния дискретности измерений на результаты их обработки. В качестве исключения можно назвать, например, монографию [25]. Хотя это явление хорошо известно в практике вычислений на компьютере, где разрядная сетка тоже (естественно) ограниченна. Это приводит к уже упоминавшемуся накоплению машинных погрешностей вычислений, для уменьшения влияния которых даже при простейших расчетах среднего арифметического (!), приходится прибегать к специальным вычислительным алгоритмам [25]. Как показывает приведенный выше краткий анализ, в общем случае дискретность приводит к занижению значений реальных погрешностей усредняемых величин.

Следствие 5. Помимо рассмотренного усреднения, другими часто встречающимися операциями вида (4.3.1) являются интерполяция, численное дифференцирование и интегрирование. Действительно, все они сводятся к некоторым линейным преобразованием величин y_i и без труда могут быть записаны в матричном виде (4.3.1). Таким образом, (4.3.3) есть решение задачи нахождения погрешности при операциях интерполяции, дифференцирования и интегрирования результатов измерений. Заметим, что в общем случае указанные погрешности становятся коррелированными, даже если погрешности исходных измерений независимы.

Следствие 6. Матрица S_X не зависит от значения вектора A₀ в уравнении (4.3.1). Полагая A₀ = AY₀, где Y₀ – некоторый вектор, состоящий из констант, получаем, что формула (4.3.3) справедлива не только для исходного вектора Y, но и для любого вектора вида Y + Y₀, т.е. ковариационная матрица ошибок вектора параметров X не зависит от прибавления к вектору измерений Y любой константы.

Следствие 7. Рассмотрим теперь нелинейную зависимость вида X = A(Y). Ее можно свести к описанной выше линейной (4.3.1) применяя линеаризацию, т.е., как и в предыдущем разделе, разлагая A(Y) в ряд Тейлора вблизи конкретного значения Y и ограничиваясь только линейными членами. Элементы матрицы A тогда будут частными производными , все постоянные члены по следствию 6 на оценку погрешности не влияют, и мы получаем все ту же формулу (4.3.3). Таким образом, например, вычислялись погрешности альбедо по полученным на втором этапе обработки результатов зондировок ковариационным матрицам погрешностей потоков в разделе 3.3, а также погрешности восстановленных параметров при решении обратной задачи для случая сплошной облачности, что будет рассмотрено в гл.6. Заметим, что для нелинейного случая соотношение (4.3.3) является лишь приближенной оценкой погрешности параметров, поскольку для точного ее получения следовало бы учесть все члены ряда Тейлора. Точность такой оценки, очевидно, тем выше, чем меньше погрешности измерений, т.е. элементы матрицы S_Y .

Вернемся теперь к решению обратной задачи и опять, для начала, ограничимся случаем линейной зависимости результатов измерений Y от искомых параметров X (4.2.3): =G₀ + GX. Пусть ошибки измерений подчиняются закону нормального распределения, плотность вероятности которого, как известно, зависит только от определенных выше , S_Y  и равна

Абстрагируемся от обсуждавшейся ранее неадекватности оператора решения прямой задачи и будем считать, что отличия результатов реальных измерений Y от вычисленных значений  обусловлены только наличием случайной погрешности. Тогда в качестве решения обратной задачи следует, очевидно, выбирать такой вектор X, которому соответствует истинное значение , т.е. =. Подставляя это условие в формулу для плотности вероятности, получим ее уже как функцию и измерений, и искомых параметров: (Y,X). Далее воспользуемся известным критерием максимального правдоподобия Фишера, согласно которому значениям искомых параметров должен соответствовать максимум совместной плотности вероятности (Y,X). Записав выражение в экспоненте явно через параметры x_k, находим максимум из уравнения , что дает систему линейных уравнений

(4.3.7)

Записывая, как в предыдущем разделе, (4.3.7) в матричной форме, получаем решение задачи

(4.3.8)

Обратим внимание, что, если положить W = S_Y^-1, (4.3.8)полностью совпадает с решением для “МНК с весами” (4.2.9). В частности, если случайные погрешности измерений некоррелированные и подчинены нормальному распределению, матрица SY диагональная и решение по МНК (4.2.9) есть одновременно и оценка максимального правдоподобия (4.3.8). Это утверждение составляет суть известной теоремы Гаусса-Маркова (см., например, [23]), она является строгим обоснованием выбора в качестве весов МНК обратных квадратов СКО измерений. Очевидно, что соотношение W = S_Y^-1 автоматически переносится и на все последующие алгоритмы МНК, описываемые выражениями (4.2.14), (4.2.17) ‑ (4.2.19), (4.2.22), (4.2.24) и (4.2.26).

Поскольку соотношение (4.3.8) имеет вид линейной связи (4.3.1) между Y и X, ковариационная матрица погрешностей восстановленных параметров S_X находится по (4.3.3). Подставляя из (4.3.8) A=(G⁺S_Y^-1G)^-1G⁺S_Y^-1 в (4.3.3) и учитывая симметрию матрицы (G⁺S_Y^-1G)^-1, получаем соотношение

Выражение (4.3.9) позволяет находить оценки погрешности восстанавливаемых параметров по известным оценкам погрешности измерений, т.е. полностью решает поставленную задачу учета последних. Для нелинейных алгоритмов (4.3.9), очевидно, сохраняет свой вид, только матрицу G следует брать на последней итерации. Заметим, что (4.3.9) относится и к методу штрафных функций (4.2.24), (4.2.26), поскольку для него, на последней итерации, дополнительный вклад в невязку (штраф), по крайней мере теоретически, должен быть нулевым, следовательно, матрица системы (4.3.1) будет такая же, как выше.

Основным этапом получения решения по МНК и методу максимального правдоподобия (4.3.8) является решение некой системы линейных уравнений, т.е. обращение ее матрицы. Однако в общем случае указанная матрица может быть очень близка к вырожденной. При реальных вычислениях на компьютере это приведет к тому, что матрица (G⁺S_Y^-1G)^-1 либо вообще “не будет обращаться”, либо при операции обращения будет внесена значительная вычислительная погрешность. Причина этого явления связана с общим свойством большинства обратных задач атмосферной оптики – их некорректностью. Глубокий теоретический анализ понятия некорректности обратных задач и многочисленные примеры подобных задач приведены в монографии [2]. В прикладном стиле, достаточно простое толкование этому явлению мы уже дали выше, при обсуждении в предыдущем разделе явления сильного разброса искомых значений при итерациях. В плане техническом некорректность проявляется в виде отмеченных трудностей обращения матрицы (G⁺S_Y^-1G)^-1, т.е. близости к нулю ее определителя. Заметим, что, вообще говоря, не все конкретные обратные задачи атмосферной оптики некорректны, однако, если их корректность не следует из теоретических соображений, следует всегда применять к ним методы решения некорректных задач, поскольку анализ корректности не выгоден в чисто техническом плане, потому что он требует значительно большего объема вычислений [2]. Таким образом, дальше мы договоримся рассматривать задачу нахождения параметров X по измерениям Y как некорректную. Ограничимся, для краткости формул, линейным случаем, а затем автоматически перенесем полученные результаты на алгоритм, рекомендованный для нелинейных обратных задач.

Методом решения некорректных обратных задач является их регуляризация – прием, в конкретном случае решения системы линейных алгебраических уравнений, сводящийся к замене исходной системы, близкой в некотором смысле к ней, для которой матрица системы всегда невырожденная [2]. Ниже мы рассмотрим два метода регуляризации, стандартно применяющихся при решении обратных задач атмосферной оптики.

Простейшим приемом регуляризации является прибавление к матрице исходной системы некоторой заведомо невырожденной матрицы. Рассмотрим вместо решения (4.3.8) решение вида

(4.3.10)

где I – единичная матрица, h – некоторый числовой параметр. Очевидно, что при h -> 0, решение (4.3.10) стремится к “настоящему” решению (4.3.8). Из этого следует простой алгоритм: находится последовательность решений (4.3.10) при уменьшении параметра h и за решение принимается значение X, которому соответствует минимум невязки. Этот прием называется регуляризация по Тихонову (хотя как эмпирический метод она была известна давно, А.Н. Тихонову принадлежит ее строгое математическое обоснование [2]).

Регуляризацию по Тихонову легко связать с рассмотренным в предыдущем параграфе методом штрафных функций. Действительно, потребуем для решения выполнения условий x_k = 0, тогда решение по методу штрафных функций (4.2.22) непосредственно переходит в (4.3.10). Конечно, нельзя требовать строгого соблюдения x_k = 0, поэтому и множитель h выбирается как можно меньшим. Таким образом, регуляризация по Тихонову соответствует наложению на решение определенного ограничения, конкретно – требование минимального расстояния решения от нуля, т.е. сужению множества возможных решений обратной задачи. В плане теоретическом к этому сводятся все приемы регуляризации. Требование x_k = 0 означает, что компоненты вектора X не должны слишком сильно различаться друг от друга, т.е. отсекает возможность получения сильно осциллирующих решений. Но ведь это фактически и есть прием борьбы с сильным разбросом решений на итерациях нелинейных задач. Действительно, в этом случае регуляризация по Тихонову применяется во всех стандартных алгоритмах современных нелинейных МНК, см., например, [15].

В рассматриваемой постановке обратных задач атмосферной оптики все искомые параметры X имеют физический смысл. Следовательно, определенная информация о них известна до проведения измерений Y, она называется априорной. Примем, что параметры X характеризуются априорным средним значением и априорной ковариационной матрицей D. Предположим, что ошибки параметров подчинены нормальному распределению, т.е.

Подчеркнем, что указанные априорные характеристики и D – это информация о параметрах, известная заранее, без рассматриваемых измерений, в частности это относится и к априорным СКО параметров X. Учитывая теперь полученную выше плотность вероятности ошибок измерений (Y,X) и полагая отсутствие корреляции между ошибками измерений и параметров потребуем выполнения для их совместной плотности (Y,X)(X) критерия максимального правдоподобия. Рассматривая для удобства в качестве независимой переменной X – и проделав все выкладки, аналогичные выводу (4.3.8), получаем :

(4.3.11)

Решение (4.3.11) известно как метод статистической регуляризации [26-28]. Регуляризация здесь достигается за счет прибавления к матрице системы обратной априорной ковариационной матрицы D^-1, действительно, легко проверить, что решение (4.3.11) существует даже в самом наихудшем случае G⁺S_Y^-1G. С другой стороны, чем больше априорное СКО параметров, тем меньше вклад D^-1 в (4.3.11) и в пределе, при D^-1 = 0 (4.3.11) переходит в решение без регуляризации (4.3.8). Статистическая регуляризация (4.3.11) значительно удобнее, чем (4.3.10), поскольку не требует никакого итерационного подбора параметра h (хотя, с другой стороны, требует знания априорной информации), поэтому в обратных задачах атмосферной оптики используется в основном она. Заметим, что для нелинейных задач, где при разложении в ряд Тейлора рассматривается именно разность параметров, зависимость решения от исчезает, т.е. статистическая регуляризация эквивалентна просто добавлению D^-1 в подлежащую обращению матрицу. При этом стандартно выбирается за нулевое приближение. Используя тождество

(4.3.12)

которое элементарно проверяется умножением обеих частей слева на G⁺S_Y^-1G + D^-1 и справа на GDG⁺ + S_Y, решение (4.3.11) для некоторых типов задач удобнее переписать в эквивалентной форме, не требующей обращения ковариационных матриц

(4.3.13)

Вычислим погрешности полученных параметров X по погрешностям измерений S_Y, т.е. апостериорную ковариационную матрицу ошибок параметров X. Здесь, по определению , где X –решение (4.3.12), а  – случайное отклонение от него, вызванное погрешностью измерений. Подставляя (4.3.13) в определение S_X, учитывая , после элементарных выкладок, получаем Заметим, что в этом выражении из априорной ковариационной матрицы параметров вычитается некоторая положительно определенная матрица, таким образом, измерения приводят к уменьшению априорного СКО параметров, что имеет простой физический смысл: измерения уточняют априорно известные значения параметров. Для дальнейшего преобразования S_X, докажем соотношение ,

для чего воспользуемся тождеством , и учтем (4.3.12). Окончательно получаем

.

(4.3.14)

Интересно заметить, что, несмотря на весьма сложный путь получения, (4.3.14) имеет ту же форму, что и (4.3.9): ковариационная матрица ошибок искомых параметров есть просто обратная матрица решаемой системы алгебраических уравнений, т.е. она “автоматически” получается в процессе вычислений.

Как отмечалось выше, полученные по (4.3.14) апостериорные СКО  всегда не превосходят априорные . Отношение этих СКО характеризует информативность измерений относительно данного параметра – чем оно меньше, тем больше в проведенных измерениях содержится информации о данном параметре. Любопытно, что для вычисления апостериорного СКО (4.3.14) в линейном случае не требуются собственно результаты измерений, достаточно иметь лишь алгоритм решения прямой задачи (матрицу G). Поэтому вычисление возможной точности восстановления параметров и оценка информативности может проводиться еще на предварительном этапе, до проведения измерений. Для нелинейных задач, где матрица производных G зависит от решения X, такие утверждения, строго говоря, неверны, тем не менее, и в этом случае часто используют (4.3.14) для анализа информативности задачи до измерений.

Выбор априорной ковариационной матрицы D представляет некоторые сложности при использовании метода статистической регуляризации. Если имеется достаточная статистика прямых измерений искомых параметров, то D просто вычисляется по ней. Иначе приходится прибегать к различным физическим и эмпирическим оценкам и моделям. Для конкретной задачи обработки результатов зондировок априорные модели будут обсуждаться в главе 5. Заметим, что при необходимости интерполяции D она элементарно пересчитывается по (4.3.3) – см. следствие 5. Следует сказать, что при расчетах ковариационной матрицы необходимо приводить результаты с достаточно высокой точностью, не округляя полученные коэффициенты корреляции. Иначе ошибки округления (следствие 4 из соотношения (4.3.3)) приводят к искажениям структуры матрицы, в результате чего возникают проблемы с ее использованием. В частности, все литературные данные по коэффициентам корреляции атмосферных параметров приводятся с точностью 2-3 знака, поэтому при попытке работы с такими матрицами они, как правило, “не обращаются”. С другой стороны, сложности с обращением D могут носить и принципиальный характер, она может быть близка к вырожденной, если искомые параметры имеют сильные корреляции между собой.

Для преодоления указанных трудностей, а также в плане общей оптимизации скорости работы алгоритма, необходимо перейти от искомых параметров к независимым, для которых отсутствуют корреляции и матрица D диагональная. Такое преобразование, как известно [11], осуществляется матрицей P из собственных векторов D, при этом матрица D перейдет в диагональную матрицу L по известным формулам преобразования координат [11] L = P^-1DP. После вычисления апостериорной ковариационной матрицы следует осуществить обратный переход к искомым параметрам P^-1S_XP, что с учетом ортогональности собственных векторов (P^-1 = P⁺) для решения (4.3.11) дает

(4.3.15)

При вычислении собственных чисел и векторов матрицы D следует использовать метод вращения [11], который, хотя и медленный, но успешно работает при близких (кратных) собственных числах. Чтобы при вычислении собственных чисел не происходила потеря точности, рекомендуется следующий прием нормировки. Примем за единицу измерения параметра x_k его априорное СКО, т.е. введем вектор  и перейдем к величинам , где матрица D' – корреляционная. Решаем обратную задачу со штрихованными переменными, а после ее решения переходим к исходным единицам измерения x_k = d_k, (S_K)_ik=(S'_X)_ikd_id_k. Заметим еще, что в силу отмеченных выше искажений ковариационных матриц при округлениях, их собственные числа могут становиться значимо отрицательными. Для борьбы с таким явлением рекомендуется регуляризация по Тихонову: вместо матрицы D' используется D'+h²I с последовательным увеличением h до исчезновения отрицательных собственных чисел.

В часто встречающемся на практике случае сильной корреляции между искомыми параметрами существенно отличны от нуля только несколько максимальных собственных чисел D, обозначим их число m. Тогда все вычисления значительно ускоряются, если в матрице L оставить только m указанных собственных чисел (она станет размерности m x m), а в P – только m соответствующих им столбцов (размерность m x K). Этот прием составляет суть известного метода главных компонент.

Обозначая полученные матрицы как Lm и Pm из (4.3.15) получаем

(4.3.16)

Используя формулу (4.3.16) вместо (4.3.15) в ряде задач иногда удается сократить объем вычислений на порядок и более.

Критерии выбора значения m в (4.3.16) могут быть различными. Математические основаны на сравнении исходной матрицы D и матрицы P_m⁺L_mP_m, которые теоретически при m=K должны совпадать. Соответственно, m выбирается исходя из допустимой величины такого несовпадения. При этом поэлементное сравнение указанных матриц вряд ли имеет смысл, обычно сравнивают диагональные элементы (дисперсии) или суммы этих элементов (инвариант при преобразованиях координат [11]). Объективный физический выбор m предложен в информационном подходе В.П. Козлова [27], хотя, заметим, он удобен не для всех типов обратных задач, поскольку требует очень громоздких вычислений. Согласно следствию 2 из (4.3.3), вариация измерений, вызванная априорными вариациями параметров, есть GDG⁺. Перейдем к собственному базису этой матрицы, т.е. к независимым вариациям измерений. Тогда собственные числа GDG⁺ есть “полезный сигнал”, который следует сравнивать с шумом – СКО измерений. Если измерения равноточные и некоррелированные с СКО равным s, то число m есть число собственных чисел, больших s². Более сложным является случай неравноточных и коррелированных измерений (именно он имеет место при обработке зондировок). Тогда измерения следует предварительно привести к независимости и единой точности s = 1. Такой переход основан на теореме об одновременном приведении двух квадратичных форм к диагональному виду [11] и осуществляется матрицей P_YL_Y^-1/2, где P_Y – матрица собственных векторов S_Y, а L_Y– соответствующая им диагональная матрица из собственных чисел S_Y. Таким образом, согласно (4.3.3) выбор m определяется числом собственных чисел матрицы P_YL_Y^-1/2GDG⁺L_Y^-1/2P_Y⁺ , больших единицы. Заметим, что для нелинейных задач матрица G меняется от итерации к итерации, но проделывать каждый раз столь громоздкие вычисления значения m нереально. Поэтому его вычисляют предварительно по матрице G₀, усиливая, для гарантии, условия выбора, т.е. сравнивая собственные числа не с единицей, а с меньшим значением.

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

(4.3.17)

Алгоритм с улучшенной сходимостью (4.2.26), конкретно использовавшийся при обработке зондировок переходит в

(4.3.18)

В обоих случаях апостериорная ковариационная матрица параметров вычисляется по формуле:

S_X= P_m⁺(P_mG_n⁺S_Y^-1G_nPm⁺ + Lm^-1)^-1 P_m .

.