Сингулярный анализ решения некоторых задач спутниковой геодезии Текст научной статьи по специальности «Математика»

Геодезия

УДК 528.1:629.7

В.И. Дударев СГГ А, Новосибирск

СИНГУЛЯРНЫЙ АНАЛИЗ РЕШЕНИЯ НЕКОТОРЫХ ЗАДАЧ СПУТНИКОВОЙ ГЕОДЕЗИИ

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

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

V.I. Dudarev SSGA, Novosibirsk

SYNGULAR ANALYSIS OF SOME SPACE GEODESY TASKS SOLUTION

The questions of the way of linear equation solution choice are considered, which provide with regular and correct problem statement of complex nonlinear dynamics state parameters evaluation in space geodesy.

dynamic system, technique of least squares, the system of nonlinear equations, condition number, singular decomposition, solution stability.

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

Решение нелинейной задачи МНК включает процедуру замены нелинейной задачи оценивания на линейную, т. е. переход от нелинейных моделей измерений к линейным. В частности, решение системы нелинейных уравнений методом Г аусса - Ньютона предусматривает последовательное выполнение ряда линейных задач до достижения принятого критерия окончания итеративного процесса. Иными словами, неоднократно решается линейная задача наименьших квадратов, суть которой заключается в следующем [2-6]. Задана несовместная система линейных алгебраических уравнений (СЛАУ)

A • X - F = V, (1)

35

Геодезия

в которой A - прямоугольная m х n матрица коэффициентов (m - число измерений; n - число оцениваемых параметров); X - n-мерный вектор-столбец неизвестных, состоящий из поправок к приближенному значению вектора оцениваемых параметров; F - m-мерный вектор-столбец правой части; V - m-мерный вектор-столбец поправок к результатам измерений.

Требуется найти такой действительный вектор X, для которого евклидова норма Ц-Ц Е вектора A • X - F будет минимальной:

min

A • X - F

Б’

(2)

где X - оценка вектора X (или вектор решений).

По отношению к несовместной системе уравнений

A • X = F (3)

V будет вектором невязок [7], значение которого определяется равенством

V = A • X - F. (4)

Г еометрически условие выражения (2) означает, что необходимо найти такую точку с координатами A • X в n системе m-мерных базисных векторов

A = [ai a2 ... an], (5)

которая ближе всего находится к вектору F по евклидовой норме. Этому требованию удовлетворяет точка, для которой вектор невязок V будет ортогонален каждому базисному вектору a i (i = 1, ..., n), т. е.

AT • (A• X-F) = 0. (6)

Матричная запись (6) представляет собой совместную СЛАУ, известную под названием линейной системы нормальных уравнений с квадратной n х n

т т

матрицей коэффициентов A A и вектором правой части A F. Решением этой системы является вектор

X = (Aт • A)-1 • Aт • F.

(7)

Обобщенным решением (псевдорешением) системы (3) в смысле наименьших квадратов называется любой вектор X , для которого евклидова норма вектора невязок V из (4) достигает наименьшего значения (2). Решение X,

имеющее наименьшую евклидову норму min

X

Б

называется нормальным

обобщенным решением (нормальным псевдорешением) [8]. Обобщенные решения и только они удовлетворяют СЛАУ (6). Если матрица A - неполного

36

Геодезия

ранга, то существует бесчисленное множество обобщенных решений. Среди них имеется нормальное псевдорешение, которое всегда существует и единственно. Как будет показано ниже, единственное решение СЛАУ (3) находится из выражения [3, 7-9]

X = A + • F, (8)

где A + - псевдообратная (либо эффективная псевдообратная) матрица для матрицы A.

На заключительном этапе линейной задачи МНК можно находить вектор X из решения СЛАУ с квадратной (7) либо прямоугольной (8) матрицей коэффициентов. Зачастую такие матрицы определены с ошибками. Кроме того, элементы этих матриц в оперативной памяти вычислительной машины представляются также с ошибками. Вполне возможно, что матрица коэффициентов будет неполного ранга (факт нерегулярной постановки измерительной задачи) или близка к этому (факт некорректной постановки измерительной задачи). А это приводит либо к неопределенности решения, либо к повышению чувствительности решения СЛАУ к ошибкам входных данных. В связи с перечисленными выше положениями на этапе постановки задачи оценивания необходимо выбирать оптимальную стратегию решения обратной задачи и правильно интерпретировать полученные результаты. Поэтому следует привлекать современные методы вычислительной алгебры, чтобы в этой части обеспечить регулярную и корректную постановку задачи оценивания расширенного вектора параметров состояния конкретной нелинейной динамической системы.

Заключительный этап выполнения линейной задачи МНК - решение СЛАУ (1) - представим в виде решающей подсистемы, у которой входом служат матрица коэффициентов и вектор правой части, а выходом - векторы X и V, и рассмотрим методику исследования ее устойчивости.

По определению А.Н. Тихонова и В.Я. Арсенина [10], решение СЛАУ будет устойчивым, если малым изменениям входных данных соответствуют малые изменения вектора решения. Решение неустойчиво, если малые изменения входных данных могут привести к большим изменениям вектора решения. Линейные системы, чувствительные к изменению входных данных, называются плохо обусловленными или неустойчивыми [8, 11]. Такие СЛАУ содержат плохо обусловленные матрицы коэффициентов [12].

При наличии возмущений ДА и AF матрицы коэффициентов и вектора правой части фактически решается система уравнений

(А + ДА) • X = F + AF, (9)

где X - n-мерный вектор-столбец решения этой возмущенной системы.

37

Геодезия

Верхний предел его относительной ошибки

X - X

-х

определяется неравенством [3, 11]

с параметрами

а =

HL <JVL: у =<

II-AIe -IIXE IIA• XIE

llE II llE

'll X E (10)

+ Y • sF]/(1 - MA) • sA ) (11)

F = |^F||E '1 lFllE ; (12)

fL .. ..

I, ; ^(A) = a L • A • XL y 11 llE IIE A+ E’ (13)

где ^(А) - число обусловленности матрицы A.

Из неравенства (11) следует, что устойчивость решения СЛАУ к изменению входных данных существенно зависит от числа обусловленности ^(А). Это число выступает в роли коэффициента увеличения совместного влияния относительных ошибок матрицы коэффициентов еА и вектора правой части sF на результат: с увеличением числа обусловленности соответственно увеличивается относительная ошибка решения sX. Если число обусловленности мало (близко к единице), то матрица A будет хорошо обусловленной. Если число обусловленности велико, то данная матрица будет плохо обусловленной [4, 8, 12].

Из анализа выражений (11) и (12) следует, что на число итераций существенное влияние оказывает точность задания матрицы A. Если эта матрица плохо обусловлена (^(А) >> 1), то необходимо принять меры к повышению точности ее расчета так, чтобы выполнялось условие: ^(А) ■ sA < 1. В этом случае будет соблюдаться неравенство sX < 1 и итеративный процесс решения системы нелинейных уравнений будет сходиться.

Переходя к методам решения линейных систем, отметим, что выражение (8) здесь и в целом ряде работ используется только для краткости записи. Обычно при решении СЛАУ псевдообратная (либо обратная) матрица не вычисляется, так как более эффективно находить решение с применением методов разложения матрицы коэффициентов на матрицы-сомножители (обычно на две или три). Для квадратных матриц чаще всего используются разложения Холесского, Грама-Шмидта, Гаусса, а для прямоугольных - Хаусхолде-ра, SVD (сингулярное разложение). Достоинство этих методов разложения состоит в том, что они позволяют путем ортогональных преобразований привести исходную матрицу к более простому виду: треугольному, диагональному и т. д. При этом ортогональные преобразования не вносят дополнительных ошибок в вычисления. Упомянутые методы разложений достаточно

38

Геодезия

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

Мощным средством вычислительной алгебры совместно с линейной задачей МНК является сингулярное разложение [5, 7-9, 11, 12]. Сингулярным разложением прямоугольной m х n матрицы A называется разложение вида

A = U • L • VT, (14)

где U и V - соответственно левая m х m и правая n х n ортогональные матрицы; L - прямоугольная m х n матрица с элементами аji = 0 при i ф j и aii = ai > 0

(i = 1, ..., n, j = 1, ..., m). Величины аг называются сингулярными числами матрицы A.

Разложение (14) можно использовать для установления степени устойчивости вектора решения к изменению исходных данных. С этой целью применяются зависимости (10)—(12), для которых число обусловленности следует вычислять по формуле [11]

^(A) а max/а min, (15)

где а max, amin - максимальное и минимальное сингулярные числа.

Сингулярное разложение сводит общее решение линейной задачи МНК к задаче с диагональной матрицей, что позволяет достаточно просто определить ранг матрицы коэффициентов и получить единственное решение СЛАУ. Кратко изложим суть этого вопроса, основываясь на работах [3, 5, 9, 12].

Применив к матрице коэффициентов A в СЛАУ (3) сингулярное разложение, получим преобразованную систему

L • Z = G. (16)

В ней n-мерный вектор-столбец преобразованных неизвестных Z и m-мерный вектор-столбец правой части G находятся из равенств

Z = VT • X, (17)

G = U1 • F

Определив псевдообратную матрицу L+ как

L'

а+ = а-1, i = j (i = 1,..., n,j = 1,..., m). o+ji = 0 (i ф j) ,

(18)

(19)

получим оценку Z вектора Z в (16):

Z = L+ • G. (20)

39

Геодезия

Теперь обобщенное решение СЛАУ (3) находится из равенства

X = V • Z. (21)

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

A + = V • £+•UT, (22)

после чего обобщенное решение можно найти по формуле (8).

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

X =

S

0 ’

S = diag(oi, ..., о п),

(23)

то удобно представить

gi

g 2 _ ’

(24)

где g ь g2 - векторы-столбцы размерности n и n-m соответственно.

Тогда вместо формулы (20) рациональнее вычислять вектор Z по формуле

Z = S-1 • gb (25)

в которой

S-1 = diag(o-1, ..., оП1). (26)

Затем уже обобщенное решение определяется из равенства (21). Очевидно, что эффективность процедуры по формуле (25) возрастает с ростом числа m уравнений относительно числа n н еизвестных.

Для обобщенного решения X из формулы (21) вектор невязок получается подстановкой выражений (14), (16), (18), (23)-(25) в (4):

V = U • (UT • F - X • Z) = U •

0 g 2

(27)

Евклидова норма этого вектора равна

(28)

40

Геодезия

Если матрица A - неполного ранга, г < n (г = rankA), то псевдообратная матрица Х+ рассчитывается следующим образом:

Е+

Г + -1 • ^

оii = оi , 1 = J , оi ^ т,

< о+ = 0, 1 = J, оi <т,

о +i = 0, 1 ф J.

(29)

Здесь параметр т отражает точность задания матрицы A и точность представления чисел в памяти конкретной вычислительной машины [5, 7, 9]:

т = оmax ' ^A • (30)

Отбрасывание сингулярных чисел, меньших т, приводит к уменьшению числа обусловленности, значение которого находится как в выражении [5]:

MA) = отах /т • (31)

С введением величины т вводится также понятие эффективного ранга, равного количеству сингулярных чисел, больших т [5].

Если псевдообратная матрица Е+ вычисляется по формуле (29), то вектор решения Z в формуле (20) лучше представить состоящим из векторов-столбцов Z1 и Z 2 размерности г и n-r соответственно. Матрицу G будут составлять векторы-столбцы g 1 и g2 размерности г и m-r. Очевидно, что преобразованной системе уравнений (16) удовлетворяет единственный вектор Z1 и произвольный вектор Z2. Такой выбор вектора решений не увеличивает евклидовой нормы вектора невязок. Отмеченная множественность решений для преобразованной системы обусловливает множественность решений X (21) исходной СЛАУ. Единственным решением задачи (2) является решение минимальной евклидовой нормы. Очевидно, что в этом случае вектор Z 2 должен содержать нулевые компоненты. Поэтому нормальное обобщенное решение будет иметь вид

X = V'

(32)

Для решения линейной системы уравнений с матрицей коэффициентов неполного ранга г < n также применима экономичная схема вычислений. При этом решение СЛАУ (16) осуществляется с использованием формулы (25), в которой матрица

S-1 = diag(o-1, ..., оГ1)

(33)

41

Геодезия

состоит из г обратных сингулярных чисел, больших величины т. Для нормального обобщенного решения (32) евклидова норма вектора невязок (28) также будет минимальной.

Заметим, что нормальное обобщенное решение может быть получено из выражения (8), в котором эффективная псевдообратная матрица A + рассчитывается по формуле (22) с использованием псевдообратной матрицы E+ из формулы (29) [5].

В работе [3] отмечается, что при соотношении в СЛАУ числа строк и столбцов p = m / n больше 2 целесообразнее к матрице коэффициентов применять ортогональное разложение Хаусхолдера. А затем уже полученную треугольную матрицу подвергать сингулярному разложению. Тогда при значении p «10 экономия времени при вычислениях может достигать 50 %.

Применение сингулярного анализа показывает, что решение систем линейных уравнений с прямоугольной матрицей гораздо предпочтительнее, чем решение этих же систем с использованием нормальной матрицы коэффициентов. Для доказательства этого утверждения выполним сингулярное разложение (14) нормальной матрицы коэффициентов в (6). Можно записать

AT • A = V• ET • UT • U• E• VT = V• ET • E• VT. (34)

Здесь матрица E • E будет состоять из квадратов сингулярных чисел матрицы E. Отсюда имеем

P(eT • E) = О2max/о2min = й2 (A), (35)

т. е. число обусловленности нормальной матрицы коэффициентов A • A в квадрат раз больше числа обусловленности прямоугольной матрицы А.

Возвращаясь к оценке относительной ошибки решения СЛАУ, подчеркнем, что вектор V не может служить надежной мерой ее решения. Анализ формулы

ех = KA) -| |V| |Е /| |F|| Е, (36)

приведенной в работе [12], позволяет сделать следующий вывод: сильно отличающийся от своего точного значения вектор решения может давать весьма малые невязки. Это отличие увеличивается с ростом числа обусловленности р(А) матрицы коэффициентов СЛАУ.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Дударев, В.И. Постановка задачи оценивания состояния динамической системы [Текст] / В.И. Дударев // Математическая обработка результатов геодезических наблюдений: межвуз. сб. науч. тр. - Новосибирск: НИИГАиК, 1993. - Т. 51(91). - С.19-30.

42

Геодезия

2. Брандин, В.Н. Экспериментальная баллистика космических аппаратов [Текст] /

В.Н. Брандин, А.А. Васильев, А.А. Куницкий. - М.: Машиностроение, 1984. - 264 с.

3. Лоусон, Ч. Численное решение задач метода наименьших квадратов [Текст] /

Ч. Лоусон, Р. Хенсон; пер. с англ. X^. Икрамова. - М.: Наука, 1986. - 232 с.

4. Форсайт, Д. Численное решение систем линейных алгебраических уравнений [Текст] / Д. Форсайт, К. Молер; пер. с англ. В.П. Ильина, Ю.И. Кузнецова. - М.: Мир, 1969. -166 с.

5. Форсайт, Д. Машинные методы математических вычислений [Текст] / Д. Форсайт,

М. Малькольм, К. Молер; пер. с англ. X^. Икрамова. - М.: Мир, 1980. - 279 с.

6. Эльясберг, П.Е. Определение движения по результатам измерений [Текст] / П.Е. Эльясберг. - М.: Наука, 1976. - 416 с.

7. Беклемишев, Д.В. Дополнительные главы линейной алгебры [Текст] / Д.В. Беклемишев. - М.: Наука, 1983. - 336 с.

8. Молчанов, И.Н. Машинные методы решения прикладных задач. Алгебра, приближение функций [Текст] / И.Н. Молчанов. - Киев: Наукова думка, 1987. - 288 с.

9. Дэннис, Д. Численные методы безусловной оптимизации и решения нелинейных уравнений [Текст] / Д. Дэннис, Р. Шнабель. - М.: Мир, 1988. - 440 с.

10. Тихонов, А.Н. Методы решения некорректных задач [Текст] / А.Н. Тихонов,

В.Я. Арсенин. - М.: Наука, 1979. - 288 с.

11. Воеводин, В.В. Матрицы и вычисления [Текст] / В.В. Воеводин, Ю.А. Кузнецов. -М.: Наука, 1984. - 320 с.

12. Хорн, Р. Матричный анализ [Текст] / Р. Хорн, Ч. Джонсон; пер. с англ. Х.Д. Икра-мова, А.В. Князева, Е.Е. Тырышникова. - М.: Мир, 1989. - 655 с.

Получено 12.05.2010

© В.И. Дударев, 2010

43