Авторефераты по всем темам >>
Авторефераты по земле
На правах рукописи
ПОПАДЬЕВ ВИКТОР ВАЛЕРЬЕВИЧ
СОВЕРШЕНСТВОВАНИЕ КОЛЛОКАЦИОННЫХ МЕТОДОВ РЕШЕНИЯ ЗАДАЧ ФИЗИЧЕСКОЙ ГЕОДЕЗИИ
25.00.32 Геодезия
Автореферат диссертации на соискание ученой степени кандидата технических наук
Москва 2012
Работа выполнена на кафедре высшей математики Московского государственного университета геодезии и картографии (МИИГАиК).
Научный руководитель доктор технических наук, профессор Нейман Юрий Михайлович Официальные Яшкин Станислав Николаевич, Московоппоненты: ский государственный университет геодезии и картографии (МИИГАиК), доктор технических наук, профессор, кафедра астрономии и космической геодезии, профессор;
Зуева Анастасия Николаевна, кандидат технических наук, ФБУ ЦНИИ Минобороны России, научноисследовательский центр топографогеодезического и навигационного обеспечения, ведущий научный сотрудник.
Ведущая организация Институт астрономии РАН
Защита состоится 18 декабря 2012 года в 12:00 ч. на заседании диссертационного совета Д 212.143.03 при Московском государственном университете геодезии и картографии (МИИГАиК) по адресу: 105064, Москва, Гороховский пер., д. 4, зал заседаний Учёного совета.
С диссертацией можно ознакомиться в библиотеке университета.
С авторефератом диссертации можно ознакомиться на сайте www.miigaik.ru.
Автореферат разослан ноября 2012 года.
Ученый секретарь диссертационного совета Климков Юрий Михайлович
Общая характеристика работы
Обоснование актуальности темы диссертации.
Благодаря трудам М. С. Молоденского, В. В. Бровара, Л. П. Пеллинена, T. Krarup, H. Moritz и других отечественных и зарубежных учёных линейные методы определения физической поверхности Земли и её внешнего гравитационного поля базируются в настоящее время на строгой теории, позволяющей решать основную проблему геодезии принципиально с любой точностью и обладающей, таким образом, практически неограниченными потенциальными возможностями. Однако вопросы устойчивой численной реализации этой теории представляют интерес и в настоящее время, поскольку неизбежные на практике дискретность измерений и наличие разного рода погрешностей вносят существенную специфику. К тому же за последние десятилетия возникла многообещающая возможность существенно расширить состав исходных данных за счёт альтиметрии и различных методов спутниковой гравиметрии, в частности градиентометрии. Одним из наиболее эффективных методов, позволяющих совместно использовать такую разнородную информацию для решения как основной задачи геодезии, так и других более частных задач физической геодезии, является метод среднеквадратической коллокации. Характеристической чертой этого метода является тот факт, что любой объект геодезических измерений трактуется как определённый функционал на одном и том же потенциале силы тяжести Земли. Это и позволяет рассматривать любые разнородные по составу геодезические измерения с единых позиций: измерены с определённой точностью значения некоторых функционалов на геопотенциале, требуется оценить значения некоторых других функционалов и их точность или восстановить потенциал. По-видимому, такая уникальная ситуация все измерения суть функционалы на единой функции имеет место только в случае геодезических измерений, и было бы неоправданно этим не пользоваться. К тому же, упомянутая единая функция обладает целым рядом полезных свойств, поддающихся теоретическому изучению в рамках специальной теории потенциала. Коллокация использует всё это полностью и даёт практический метод оценки одних функционалов на геопотенциале по результатам измерения других (локальная задача коллокации) или метод восстановления самого потенциала по результатам измерения определённых функционалов на нём (глобальная задача коллокации).
Но, конечно, присущи методу коллокации и недостатки. Поэтому совершенствование этого метода и исследование возможностей его развития является актуальной задачей современных численных методов физической геодезии.
Цель и основные задачи исследования Основной целью диссертационной работы является исследование и разработка таких модификаций метода коллокации, которые частично или полностью преодолевают следующие основные недостатки этого метода:
необходимость решать количество уравнений, равное количеству исходных измерений;
необходимость опираться на гипотезу о стационарности и даже изотропности гравитационного поля Земли (ГПЗ).
Для достижения этой цели решались следующие задачи:
1. Анализ и экспериментальные исследования новейших зарубежных разработок последовательной и быстрой коллокации.
2. Разработка и исследование собственного двухэтапного гармонического анализа ГПЗ.
3. Исследование возможностей фурье-анализа и вейвлет-анализа для выявления признаков нестационарности ГПЗ.
4. Разработка и исследование методов коллокации в условиях нестационарного ГПЗ.
5. Разработка и исследование возможностей определения функций влияния исходных функционалов.
6. Адаптация программного пакета Gravsoft и использование возможностей программgeocolиsphgricдля уточнения гравитационного поля по результатам спутниковой градиентометрии.
Научная новизна и практическая значимость работы.
1. Совмещение различных идей, связанных, с одной стороны, с аппроксимацией функционалов на геопотенциале, а с другой стороны с быстрым преобразованием Фурье, позволило C. C. Tscherning (Дания) и F. Sanso (Италия) создать совсем новый алгоритм быстрой коллокации. Проделанные нами численные эксперименты с быстрой коллокацией показали существенное увеличение быстродействия при решении ресурсоёмких задач космической гравиметрии, что позволяет рекомендовать алгоритм быстрой коллокации для практического использования и в нашей стране.
2. К этому же направлению исследований принадлежит и наша разработка в виде двухэтапного гармонического анализа. Вычисления вдоль параллелей методом быстрого преобразования Фурье позволяют обеспечить некоррелированность коэффициентов Фурье, соответствующих различным частотам.
Поэтому на втором этапе вместо решения одной, но очень большой системы mp уравнений, доказана возможность решать 2N+1 систем, но каждая из них содержит лишь p уравнений. Здесь m и p число меридианов и параллелей, соответственно, определяющих сетку с исходными данными, а N наивысшая степень искомых гармонических коэффициентов. Разработаны методы контроля вычислений, составлены необходимые программы для ПК.
3. Ковариационный анализ гравитационного поля, обычно выполняемый в процессе решения задач физической геодезии методом коллокации, рекомендуется, для выявления анизотропности, выполнять с помощью двумерного преобразования Фурье и дополнять построением ковариационных карт и зависимостей радиуса корреляции от азимута по образцу рис. 7 12 (см. цвет.
вклейку на стр. 11 12).
4. Предложенные операции вейвлет-анализа позволяют эффективно выявлять локальные нестационарности и обоснованно делить обширную область с нестационарным гравитационным полем на такие подобласти, в которых поле можно считать стационарным, см. рис. 9, 10 и 11 (см. цвет. вклейку на стр. 13).
5. Доказано, что отличие выборочного спектра от непрерывного имеет определённый порядок малости относительно шага сетки, и поэтому при вычислении выборочного спектра можно пользоваться известным в теории аппроксимации правилом Рунге.
6. В тех случаях, когда исходными данными в задаче коллокации служат измеренные значения одноимённых функционалов, разработана процедура непосредственного вычисления функций влияния исходных функционалов, без решения систем уравнений. Показано, что функции влияния не зависят от вида функции, на которой заданы функционалы, а зависят только от используемого гильбертова пространства. Если поле стационарно, то достаточно найти только одну функцию влияния, а все остальные функции получаются путём смещения аргумента. Область с нестационарным полем трактуется как объединение подобластей с полем стационарным (сегментация), и для каждой такой подобласти материнскую функцию влияния надо определять отдельно.
7. Использование специального вида ковариационной функции нестационарного ГПЗ обеспечивает выполнение непрерывной сегментации ГПЗ и таким образом практически снимает ограничительное требование метода коллокации о стационарности поля.
8. Алгоритм непосредственного вычисления функций влияния исходных функционалов, снимающий проблему решения большого количества уравнений.
Основные результаты, выносимые на защиту:
1. Алгоритм двухэтапного гармонического анализа ГПЗ при большом количестве исходных данных.
2. Методы ковариационного анализа локального ГПЗ, позволяющие выявлять, графически интерпретировать и учитывать в вычислениях анизотропность ГПЗ.
3. Методы вейвлет-анализа локального ГПЗ, позволяющие выявлять и графически интерпретировать участки нестационарности поля.
4. Метод создания непрерывной сегментации ГПЗ, позволяющий при использовании коллокации отбросить гипотезу о стационарности ГПЗ.
5. Алгоритм непосредственного вычисления функций влияния исходных функционалов, снимающий проблему решения большого количества уравнений.
Вклад автора в исследование.
Основные результаты диссертации получены автором самостоятельно. Отдельные результаты получены совместно с к. ф.-м. н. Сугаиповой Л. С.
Структура диссертации и объем работы.
Диссертация состоит из введения, трёх глав основного текста и заключения общим объёмом 143 стр. машинописного текста, имеется 41 иллюстрация и 11 таблиц. Список литературы насчитывает 75 наименований, в том числе 48 на английском и других языках.
Публикации и апробации работы.
Результаты исследований представлялись на научно-технических конференциях студентов, аспирантов и молодых учёных МИИГАиК 4 апреля 20г., 6 апреля 2011 г., 3 апреля 2012 г. и опубликованы в трёх статьях в журналах, включённых в перечень ВАК.
Основное содержание работы
Во введении обоснована актуальность темы, сформулированы цель и основные задачи исследования, указаны научная новизна и практическая значимость ожидаемых результатов.
Первая глава Основы среднеквадратической коллокации начинается с краткого исторического обзора развития метода коллокации.
Благодаря работам датского математика T.Krarup, шведского геодезиста A. Bjerhammar, советского геофизика В. И. Аронова, австрийского геодезиста H. Moritz, датского геодезиста C. C. Tscherning и других ученых, среднеквадратическая коллокация в настоящее время является одним из самых эффективных численных методов физической геодезии. Этому способствуют следующие факты:
возмущающий потенциал, определяемый в результате решения глобальной задачи коллокации, является гармонической функцией, то есть автоматически решается уравнение Лапласа, обусловленное физикой явления;
решениям локальных задач коллокации обеспечены наилучшие линейные несмещённые оценки;
имеется теоретически обоснованная возможность совместно обрабатывать разнородные исходные данные, полученные как на поверхности Земли, так и вне её;
имеется теоретически обоснованная возможность приписывать адекватную меру точности получаемым оценкам.
Но, конечно, присущи методу коллокации и недостатки. Основной из них состоит в необходимости решать такое количество линейных алгебраических уравнений, которое при строго формальном подходе равно количеству измерений. В современных задачах спутниковой гравиметрии это количество, как известно, достигает десятков миллионов и решить совместно такое количество уравнений пока невозможно. К тому же при увеличении числа исходных измерений обычно не только растёт подлежащая решению система линейных уравнений, но и ухудшается её обусловленность. Однако совсем необязательно составлять все нужные уравнения одновременно и решать их все сразу.
Целесообразно пользоваться различными многогрупповыми методами и алгоритмом так называемой быстрой коллокации (по аналогии с быстрым преобразованием Фурье). Кроме того, основной путь решения данной проблемы усреднение измерительной информации по ячейкам подобранных размеров и формы.
Другим серьёзным недостатком метода коллокации является предположение о том, что ковариация между точечными значениями трансформант ГПЗ не зависит от координат точек, а зависит только от вектора, соединяющего пару точек или даже только от модуля этого вектора. Это предположение является следствием гипотезы о стационарности или даже изотропности гравитационного поля. Конечно, в общем случае это не так, и использование в коллокации единой ковариационной структуры приводит к излишнему заглаживанию ГПЗ в горах и, наоборот, к неоправданным скачкам результатов в равнинных районах. И хотя известны доказательства того, что эти неприятные влияния уменьшаются с увеличением плотности исходных данных и вообще вариации воспроизводящего ядра не сильно влияют на оценки искомых функционалов, а больше сказываются на результатах оценивания точности, тем не менее, гипотезу о стационарности поля следует считать принципиальным недостатком метода коллокации.
Наконец, некоторые исследователи считают недостатком тот факт, что перед вычислениями методом коллокации обычно необходимо выполнять ковариационный анализ гравитационного поля, с которым предстоит работать, и определять подходящую ковариационную функцию. В связи с этим стоит вспомнить, что трансформанты ГПЗ представляют собой непрерывные функции, а измерения выполняются в дискретных точках. Поэтому все практические методы физической геодезии в той или иной мере используют интерполяцию. Но при отсутствии пространственной корреляции между данными любая интерполяция просто не имеет смысла.
С математической точки зрения метод коллокации самодостаточен в том смысле, что вполне обеспечивает согласованность между всевозможными трансформантами ГПЗ при любом воспроизводящем ядре. Но результаты ковариационного анализа позволяют уточнить то гильбертово пространство, в котором наиболее целесообразно искать решения, и, следовательно, задать наиболее подходящую меру аппроксимации (то есть правило измерения расстояния между функциями аппроксимируемой и аппроксимирующей) именно в имеющихся конкретных обстоятельствах. Подобная априорная оценка гладкости искомого решения свойственна многим методам аппроксимации, поскольку невозможно действовать разумно без информации о том, что решение принадлежит определённому классу и, следовательно, удовлетворяет определённым ограничениям.
Изложены основы теории среднеквадратической коллокации в рамках традиционной гипотезы об изотропности ГПЗ. Описаны различные трактовки метода: с функциональной и вероятностно-статистической точек зрения.
Разобран смысл часто смешиваемых терминов коллокация и кригинг.
Собраны и систематизированы используемые в геодезии ковариационные функции и спектральные плотности различных трансформант стационарного геопотенциала.
В заключение отмечается, что в разных науках о Земле часто используются похожие на коллокацию численные методы под разными названиями. Во избежание путаницы, представляется целесообразным пользоваться термином среднеквадратическая коллокация только по существу, то есть только в тех случаях, когда обработке подлежат измеренные значения разнородных функционалов на некоторой функции, принадлежащей определённому гильбертову пространству, и задача состоит или в восстановлении этой функции, или в оценке других функционалов на ней. В остальных случаях лучше говорить об интерполяции, экстраполяции, разного рода кригинге, регрессии или об уравнивании методом наименьших квадратов, поскольку если объекты исходных измерений не являются функционалами на одной и той же функции, то между ними нет зависимостей, обусловленных физикой соответствующего явления и, следовательно, дополнительно возникает непростая задача разумного назначения весов разнородной информации, различающейся не только точностью измерений, но и единицами измерений. Игнорирование физической размерности при этом может легко привести к ошибочным выводам.
Вторая глава Развитие теории коллокации является основной в диссертации и посвящена развитию метода коллокации и преодолению её недостатков.
Начнём с обсуждения стационарности и нестационарности гравитационного поля в задачах коллокации. Давно известно, что реальное ГПЗ редко бывает стационарным и к тому же изотропным. Поэтому в диссертации выполнен обзор известных способов ослабления влияния нестационарности поля: удаление тренда, разбиение на блоки (кусочная сегментация), пространственное преобразование данных. В данной работе мы не стремимся устранить неизотропность поля, а пытаемся её учесть. Для этого традиционный для коллокации предварительный ковариационный анализ гравитационного поля мы рекомендуем проводить по профилям различных направлений и отражать в виде карт ковариаций (термин заимствован из литературы по геостатистике). Рис. 1 4 содержат один из результатов такого анализа, наглядно показывающего структуру анизотропности.
Рис. 5, 6 (см. цвет. вклейку на стр. 11) демонстрируют зависимость радиуса корреляции от азимута.
Для того чтобы практически учесть выявленную анизотропность, необходимо определить ковариацию как функцию не только расстояния между точками, но и направления между ними. Это можно сделать с помощью двумерного преобразования Фурье имеющегося массива исходных данных g(x, y). Полученную матрицу G(x, y) надо умножить на такую же матрицу с комплексно сопряжёнными элементами и разделить на общее количество исходных данных. Соответствующая двумерная ковариационная функция получается обратным двумерным преобразованием Фурье: C(x, y) = -F S(x, y). Пример результата таких вычислений показан на рис. 7 для того же миллионного листа карты N45 с аномалиями Фая.
То же для листа карты O41 показано на рис. 8.
Полученные ковариационные функции являются функциями разности координат x, y двух точек, то есть функциями вектора, соединяющего эти две точки, и, следовательно, описывают неизотропное поле.
Если поле не является стационарным, то есть его спектральный состав изменяется с местоположением, то для изучения структуры поля естественно отказаться от фурье-анализа и воспользоваться вейвлет-анализом.
На рис. 9 показан результат одномерного непрерывного вейвлет-преобразования по параллелям листа карты N45 с аномалиями Фая (с запада на восток) на базе вейвлета sym2 (семейство symlets почти симметричная модификациия классических вейвлетов Daubechies).
Исходной информацией послужила матрица размерностью 200300. Преобразование выполнено построчно с запада на восток. Вверху все 200 строк матрицы с данными, а внизу в увеличенном виде (как под микроскопом) только первые две строки. Подробная вейвлет-спектрограмма демонстрирует мельчайшие детали частотного образа: в её нижней части хорошо просматривается структура высокочастотных компонент, а в верхней низкочастотных (изменения яркости менее частые, чем в нижней части); при этом отчетливо фиксируются начало и конец импульсов, темный тон соответствует переходам аномалии через нуль, а светлый тон экстремумам.
На рис. 11 показаны вейвлет-спектрограммы, соответствующие коэффициентам cA1, cH1, cD1, cV1 одного уровня двумерного дискретного вейвлетпреобразования миллионного листа карты N45 с аномалиями Фая. Cпектрограмма, например, для преобразования cA1 вычисляется по следующему правилу: 1) оценивается спектр, то есть каждый коэффициент массива возводится в квадрат, 2) оценивается общая мощность спектра, то есть составляется сумма квадратов всех коэффициентов, 3) вычисляется процентное содержание квадрата каждой амплитуды, то есть квадрата каждого коэффициента, в полученной общей мощности спектра. Аналогично получены вейвлетспектрограммы для коэффициентов cH1, cD1, cV1. Cпектрограмма cA1 отражает низкочастотную структуру поля, а спектрограмма cD1 высокочастотную.
Наглядно видно, что изучаемое поле можно трактовать как стационарное, за исключением небольшой подобласти, спектр которой заметно отличается от среднего. Ориентация этой подобласти соответствует наклону анизотропности, выявленному ранее, см. рис. 1 и 5.
Таким образом, даже простейшие методы вейвлет-анализа работают как своеобразный математический микроскоп и позволяют эффективно выявлять тонкую структуру имеющейся информации (скачки, резкие переходы производных через нуль и т. п.). Подобный анализ очень полезен для изучения локальных нестационарностей и обоснованного деления большой области с нестационарным гравитационным полем на такие подобласти, в которых поле можно считать стационарным (сегментация поля).
Итак, пусть ГПЗ в изучаемом районе D не только неизотропное, но и нестационарное вовсе. С математической точки зрения это означает, что ковариационная функция Rj(p, p) представляет собой функцию четырёх переменных x, y, x, y (на плоскости) и, таким образом, полностью меняется с изменением местоположения точек p, p. Однако естественно ожидать, что эти изменения происходят достаточно плавно и, следовательно, всегда можно разбить изучаемую область на некоторое количество таких подобластей D1, D2,..., в которых ГПЗ достаточно обоснованно можно полагать стационарным (хотя, быть может, и неизотропным). Количество таких подобластей зависит, конечно, от тщательности предварительного учёта многих факторов, в том числе от структуры и неоднородности рельефа, но в любом случае эти подобласти не должны пересекаться, а их объединение должно совпадать с D.
Для каждой такой подобласти с номером i = 1, 2,... надо отдельно определить параметры анизотропности i, d1i, d2i и составить симметричную положительно определённую матрицу анизотропности подобласти Di вида d2 cos2 i + d2 sin2 i (d2 - d2 ) cos i sin i 1i 2i 1i 2i i =. (1) (d2 - d2 ) cos i sin i d2 cos2 i + d2 sin2 i 1i 2i 2i 1i Здесь азимут i определяет направление анизотропности, см. рис. 1 8, а d1i, d2i радиусы корреляции в двух взаимно ортогональных направлениях, соответствующих собственным векторам. Ковариационная функция нестационарного поля для произвольной пары точек pi(xi, yi) Di и pj(xj, yj) Di COV(,s) 2.1871.65400 320.18 -100 -1.5 -1 -0.5 0 0.5 1 1.50 100 150 2, rad Рис. 1: Ковариационная карта (мГал2) мил- Рис. 2: Изолинии ковариационной карты лионного листа карты номенклатуры N45 с (мГал2) аномалий Фая на листе N45 в сианомалиями Фая стеме координат азимут , расстояние s.
COV(,s) 2.1551.70 443320.21-1.5 -1 -0.5 0 0.5 1 1.20 40 60 80 100 120 140 160 180 2, rad Рис. 3: Ковариационная карта (мГал2) мил- Рис. 4: Изолинии ковариационной карты лионного листа карты номенклатуры O41 с (мГал2) аномалий Фая на листе O41 в сианомалиями Фая стеме координат азимут , расстояние s.
90 0.6 120 60 120 0.4 150 30 150 0.2 180 0 180 210 330 210 3240 300 240 3270 2Рис. 5: Зависимость радиуса корреляции от Рис. 6: Зависимость радиуса корреляции от азимута на листе карты N45. азимута на листе карты O41.
s, s, 12212132132214332434554556676677777722232224342232234334445555Рис. 7: Ковариация неизотропного поля аномалий Фая на листе карты N45.
Рис. 8: Ковариация неизотропного поля аномалий Фая на листе карты O41.
Analyzed all rows (points 1 - 60000).
---0 1 2 3 4 5 x Coefficients of the Continuous Transform 240 110 1 2 3 4 5 Space x Analyzed 1 and 2 rows (points 1 - 600) ---0 50 100 150 200 250 300 350 400 450 500 550 6Coefficients of the Continuous Transform 0 100 200 300 400 500 6Space Рис. 9: Результаты одномерного непрерывного вейвлет-преобразования (спектрограммы) по параллелям листа карты N45 с аномалиями Фая (с запада на восток). Вверху: все 2строк, внизу: первые две строки.
Вейвлет-спектрограмма. Вейвлет-спектрограмма.
Доля энергии каждого коэффициента cA1 Доля энергии каждого коэффициента cH0.Approximation A1 Horizontal Detail H1 0.0.0.120 10.40 50 100 50 10.0.80 60 100 10.06 0.60 0.80 1140 0.0.20 100 12250 100 150 50 100 1З - В З - В 100 200 3100 200 3Вейвлет-спектрограмма. Вейвлет-спектрограмма.
Доля энергии каждого коэффициента cV1 Доля энергии каждого коэффициента cDVertical Detail VDiagonal Detail D0.0.0.0.0.40 0.0.100 0.10.0.1180 0.0.0.200 12150 100 150 100 1100 200 3100 200 300 З - В З - В Рис. 10: Результат двумерного дискретного Рис. 11: Вейвлет-спектрограммы. Доля (в вейвлет-преобразования миллионного листа процентах) квадрата амплитуды каждого карты с аномалиями Фая. коэффициента в общей мощности спектра результатов двумерного дискретного вейвлет-преобразования первого уровня миллионного листа карты N45 с аномалиями Фая.
Scale Scale С - Ю С - Ю С - Ю С - Ю Q = 3.1725x2 + 1.2384x1 x2 + 2.4575x1 Vi0.0.V0.0.i-0.-0.-0.-0.--1 -0.5 0 0.5 alpha = Рис. 12: Линии уровней 0,25, 1,00 и 3,00 квадратичной формы (3) с 1 = 3, 53, 2 = 2, 10, = /6, x = x1, y = x2. Также показаны координатные орты i1, i2 и собстенные векторы V1, V2 матрицы -1.
Рис. 13: Дисперсия аномалии силы тяжести в единицах мгал2, вычисленной на поверхности Земли и сглаженной по одноградусным трапециям.
x0.0.0.области D имеет вид -1/ - i + j i + j C(pi, pj) = C0 |i|1/4|j|1/4 exp -(-)T (-), p- ppj pj i i 2 (2) где C0 дисперсия, а индексы i, j соответствуют номерам подобластей D1, D2,.... Важнейшую роль при этом играет квадратичная форма -Qij = (-)T (i + j) (-). (3) p- ppj pj i i Квадратичную форму (3) удобно геометрически интерпретировать соответствующими линиями уровня. На рис. 12 показаны три такие линии уровней 0,25, 1,00 и 3,00 квадратичной формы (3) с 1 = 3, 53, 2 = 2, 10, = /6, x = x1, y = x2. Там же показаны координатные орты i1, i2 и собственные векторы V1, V2 матрицы -1.
Если положить Qij = (-)T -1(-), (4) p- ppj pj i i где произвольная положительно определённая матрица, то получим известное в математической статистике понятие расстояние Mahalanobis rij = Qij. (5) Соответствующая ковариационная функция остаётся стационарной, но уже неизотропной. Если матрица = E единичная, то получается обычное расстояние в евклидовом пространстве:
rij = |-|. (6) ppj i Если = 0, то ковариационная функция способна описывать неизотропное поле при условии, что направление анизотропности совпадает с направлением координатных осей. Если 1 = 2, то Qij не зависит от и соответствующая ковариационная функция предназначена для характеристики только изотропных полей.
Интересно, что с помощью замены r2 = (-, p- на Qij можно получить p- -) pj pj i i целое семейство нестационарных ковариационных функций.
Таким образом удаётся избегать тех проблем с плохой стыковкой результатов коллокации, которые обычно возникают на границах подобластей.
Сама идея такого деления поля, конечно, не нова и обычно называется сегментацией поля. Однако традиционная сегментация кусочна в том смысле, что на границах подобластей полученные результаты коллокации не согласуются. В диссертации описана непрерывная сегментация, при которой указанные проблемы не возникают.
Далее речь идёт об ослаблении недостатка коллокации, связанного с необходимостью решать большое количество уравнений. Направления возможных действий при этом условно можно разделить на два пути. Основы первого пути состоят в разработке оптимального порядка использования исходных данных. Хорошим примером успеха на этом пути являются практически неизвестные в нашей стране последовательная коллокация и особенно развитая F. Sanso (Италия) и C. Tscherning (Дания) быстрая коллокация, названная по аналогии с известным быстрым преобразованием Фурье.
Совмещение различных идей, связанных, с одной стороны, с аппроксимацией функционалов на геопотенциале, а с другой стороны с быстрым преобразованием Фурье, позволило заметно упростить метод коллокации, в частности, процесс гармонического анализа ГПЗ.
Примерно этому же направлению исследований принадлежит и наша разработка в виде двухэтапного гармонического анализа, сокращающая количество подлежащих решению уравнений. Задача состоит в определении безразмерных гармонических коэффициентов cnk, snk разложения возмущающего потенциала T в ряд по полностью нормированным шаровым функциям (квадрат нормы = 4). Решение основано на следующих рассуждениях.
Пусть произвольная функция f(, ), квадратично интегрируемая на единичной сфере , представлена рядом по полностью нормированным сферическим функциям. Соответствующие гармонические коэффициенты определяются известными формулами 2 ank cos k = Pnk(cos ) f(, ) d sin d. (7) bnk sin k 4 0 Здесь внутренний интеграл имеет форму коэффициента Фурье при разложении функции f(, ) по параллели в тригонометрический ряд Фурье, то есть, если 2 cos k Ak() f(, ) d =, (8) sin k Bk() k то f(, ) = [Ak() cos k + Bk() sin k]. (9) k=Подставим (9) в (7) и учтем, что, благодаря ортогональности функций тригонометрической системы, 2 cos k Ak() [Ak () cos k + Bk () sin k] d = k. (10) sin k Bk() k=В результате получается связь гармонических коэффициентов и широтных коэффициентов Фурье:
k ank Ak() = Pnk(cos ) sin d; (11) bnk Bk() порядок k гармонического коэффициента совпадает с частотой коэффициентов Фурье.
Если вместо произвольной функции f(, ) взять, для определённости, вторую радиальную производную Trr(, ) возмущающего потенциала T, то гармонические коэффициенты ank, bnk для Trr можно вычислить в два этапа:
1) интегрированием по меридианам получаем широтные коэффициенты Фурье для всевозможных частот k = 0, 1, 2,...
2 Ak() cos k = Trr(, ) d; (12) Bk() sin k k 2) интегрированием по параллелям получаем гармонические коэффициенты k ank Ak() = Pnk(cos ) sin d, (13) bnk Bk() тогда n+3 ank R2 r cnk =. (14) snk bnk (n + 1)(n + 2) R Пользуясь правилом преобразования ковариаций и учитывая ортогональность синусов и косинусов, мы получили ковариации cnk с Al() и snk с Bl():
n+(n + 1)(n + 2) R cnk Al() cov, = n(T ) Pnk(cos ). (15) snk Bl() (2n + 1)R2 r Таким образом, гармонические коэффициенты tnk различных степеней n коррелируют только с такими коэффициентами Фурье Al(), Bl(), частота которых l совпадает с порядком гармонического коэффициента k.
В результате вычисления вдоль параллелей методом быстрого преобразования Фурье позволяют обеспечить некоррелированность коэффициентов Фурье, соответствующих различным частотам. Поэтому на втором этапе вместо решения одной, но очень большой системы mp уравнений, доказана возможность решать 2N + 1 систем, но каждая из них содержит лишь p уравнений. Здесь m и p число меридианов и параллелей соответственно, определяющих сетку с исходными данными, а N наивысшая степень искомых гармонических коэффициентов. Разработаны методы контроля вычислений, составлены необходимые программы для ПК.
Второй наш путь избегать большое количество уравнений принципиально отличается от первого. При его реализации необходимость в решении систем уравнений вообще не возникает и заменяется отысканием функций влияния для исходных функционалов. Основные вычисления выполняются в частотной области, а исходной информацией о гравитационном поле служит не ковариационная функция, а спектральная плотность.
Идея такого подхода вытекает из связи коллокации с кратномасштабным вейвлет-анализом, указанной C. Kotsakis (Греция). В наших исследованиях мы не пользуемся кратномасштабным анализом, но заимствуем критерий оптимальности коллокации в частотной области и показываем возможность определения под этим условием функций влияния для исходных функционалов.
Приведём рассуждения в рамках функций одной переменной, поскольку обобщение для функций нескольких переменных носит чисто технический характер.
Наилучшие оценки методом коллокации в частотной области при заданном шаге h предложено искать под условием минимума среднего интегрального значения квадрата модуля разности образов Фурье функций аппроксимируемой g(x) и аппрокимирующей (x). В терминах функций одной перменной это выглядит следующим образом:
h/ D(, h) |G() - (, h, x)|2 dx min для . (16) h -h/Здесь |G()|2 = G() G() спектральная плотность функции g(x), а (, h, x)|2 спектральная плотность функции (x). Усреднение спектральной плотности разности образов Фурье функций аппроксимируемой и аппрокимирующей по промежутку длиной в период h объясняется тем фактом, что при заданном шаге h сетку можно выбрать бесчисленным количеством способов, смещая начало отсчёта переменной x по промежутку длиной в период h. Здесь (, h, x) обозначает образ Фурье аппроксимирующей функции, сдвинутой по оси x на величину x, 2k (, h, x) = (, h) exp -i x, -h/2 < x < h/2. (17) h Соотношение (16) определяет среднеквадратический принцип минимума в частотной области. Величина D(, h) обусловлена усреднённой спектральной плотностью погрешностей аппроксимации. Однако усреднение производится не по ансамблю из гипотетических повторных измерений, а по всевозможным дискретизациям при заданном уровне разрешения h. Тем самым обеспечивается чисто детерминированный подход к задачам коллокации, несмотря на общность терминологии с теорией случайных функций.
Если исходными данными служат значения одноимённых функционалов (например значения аномалии силы тяжести, значения каких-нибудь производных геопотенциала и т. п.), то глобальную задачу коллокации можно решать в виде n- h(x) = g(kh) vh(x - kh), (18) k=где g(kh) исходные измерения изучаемой функции g(x) в узлах kh регулярной сетки с шагом h, а vh(x - kh) соответствующие функции влияния.
Практически надо вычислить только одну функцию с помощью одного (обратного) преобразования Фурье - vh(x) = F [V (, h)], (19) так как все остальные функции влияния для каждого узла kh можно получить аналитически или численно путём сдвига найденной материнской функ ции vh(x) по оси абсцисс на kh. Образ Фурье Vh() требуемой функции vh(x) определяется отношением |G()|2 |G()| Vh() = = h. (20) + |Gs(, h)|2k |G( + )|h k=- Таким образом, единственной необходимой информацией об изучаемой функции g(x) служит её спектральная плотность |G()|2 = G() G(), то есть образ Фурье соответствующей автоковариационной функции cov(x, x) = cov(x - x) = cov(x):
+ |G()|2 = cov(x)exp(-ix)dx, (21) 2 - а |Gs(, h)|2 спектральная плотность выборочной функции gs(x) + + 1 Gs(, h) = gs(x)exp(-ix)dx = G( - ks). (22) 2 h k=- - Отметим, что искомые функции влияния не зависят от вида изучаемой функции g(x), а зависят только от используемого гильбертова пространства H.
Точность решения задачи коллокации определяется дисперсией + Vh() 2(h) = |G()|2 1 - d, (23) h - которая зависит только от шага h, то есть от разрешающей способности выборки.
Чем ближе спектр выборочной функции к спектру реальной исходной функции, тем меньше дисперсия аппроксимации и тем оптимальнее используемая величина шага.
Предположим теперь, что нас интересует не сама восстанавливаемая функция g(x), а её какое-нибудь линейное преобразование типа g(x) = L(g)(x). В таком случае n-1 n- h(x) = g(xk) L(vk,h)(x) = g(xk) vk,h(x). (24) k=0 k=Если преобразование L представляет собой свёртку с какой-нибудь весовой функцией s(x) (например g(x) аномалия силы тяжести, а нужно определять высоты геоида), то есть g(x) = L(g)(x) = g(x)s(x - x) dx = g s, (25) - то образ Фурье новой функции влияния vk,h(x) находится по формуле |G()|2 |S()|2 |G()|2 |S()| V () = = h. (26) h + |Gs(, h)|2k |G( + )|h k=- Здесь S() обозначает образ Фурье весовой функции s(x) свёртки (например, известный образ Фурье функции Стокса).
Третья глава описывает численные эксперименты, выполненные в процессе подготовки диссертации, и содержит следующие разделы: 3.1 Численный эксперимент с двухэтапным гармоническим анализом; 3.2 Численные эксперименты с определением функций влияния функционалов; 3.3 Опыт работы с современным программным обеспечением (речь идёт о пакетах GravSoft и GOCE User Toolbox (GUT)); 3.3.1 Численные эксперименты с результатами спутниковой градиентометрии; 3.3.2 Создание регулярной сетки точечных и усреднённых значений различных трансформант геопотенциала.
Последний раздел отражает практический опыт диссертанта работы с гравиметром.
Проделанные эксперименты подтвердили теоретические результаты второй главы диссертации. Кроме того, установлено следующее.
1. Гладкость изучаемой функции оказывает заметное влияние на точность определения выборочного спектра.
2. Отличие выборочного спектра от непрерывного имеет определённый порядок малости относительно шага сетки и поэтому при вычислении выборочного спектра можно пользоваться известным в теории аппроксимации правилом Рунге.
3. Точность спектрального синтеза заметно повышается, если предварительно выполнить кусочно-линейную интерполяцию исходных значений функции.
4. Модельные вычисления с различными комбинациями данных вторых производных возмущающего потенциала Trr, Tee, Tnn, Tnr, Tn, Tr, T показали пригодность и достаточное удобство используемых алгоритмов и программного обеспечения. Алгоритм быстрой коллокации снижает временные затраты на вычисления примерно на порядок. Теоретически возможно совместно использовать все 6 видов измеренных производных и определять коэффициенты до 360-ой степени, но мощности доступного компьютера явно не хватает необходимо распараллеливание вычислительных потоков.
5. Некоторые типы возможных исходных данных так сильно физически коррелированы, что приводят к практически идентичным результатам; из трёх диагональных элементов Trr, Tee, Tnn матрицы Гессе наилучшие по точности результаты дают Trr; исходные значения Tee и Tnn по точности результатов примерно одинаковы, но их совместное использование оказывается полезным.
6. Хотя низкочастотная часть спектра (2 24 или 2 36) удалялась из исходных данных, оценки соответствующих коэффициентов по результатам гармонического анализа значимо отличались от нуля. Этот эффект проявляется тем заметнее, чем выше порядок исследуемой производной потенциала.
По-видимому, это является следствием наложения более высоких частот.
7. Процедура определения гармонических коэффициентов геопотенциала по результатам спутниковой градиентометрии целесообразно разделить на две части: создание на некоторой сфере (радиус которой выбирается с учётом усреднённой высоты полёта спутника) регулярной сетки с усреднёнными по ячейкам результатами измерений вторых производных геопотенциала (первичная обработка); гармонический анализ на этой сфере с использованием известных преимуществ за счёт симметрии узлов регулярной сетки.
8. Для строгого оценивания точности значений различных трансформант геопотенциала, вычисленных гармоническим синтезом по глобальным моделям GOCE, теперь доступны ковариационные матрицы ошибок соответствующих гармонических коэффициентов. Вычисления удобно выполнять с помощью программcovhsmp(вычисляет дисперсии заданных значений на сетке, регулярной по широте, долготе и постоянной высоте над эллипсоидом) иcovhs2p(оценивает и ковариации заданных пар). На рис. 13 показан результат одного из экспериментов по определению точности аномалии силы тяжести, сглаженной по одноградусным трапециям на Земле, по результатам гармонического синтеза вторых производных геопотенциала до 150-ой степени и порядка по данным проекта GOCE.
9. С помощью преобразования ковариаций и теоремы о выборке можно аналитически найти ковариационные функции между заданными точечными и искомыми усредненными по сферическим ячейкам значениями вторых производных геопотенциала.
Заключение Ранее отмечалось, что среднеквадратическая коллокация является в настоящее время одним из самых эффективных численных методов физической геодезии. Однако необходимость опираться на гипотезу о стационарности и даже изотропности гравитационного поля Земли (ГПЗ) существенно ослабляет теорию этого метода, а необходимость решать количество уравнений, равное количеству исходных измерений, заметно затрудняет его практическое использование. В данной работе сделана попытка преодолеть или, по крайней мере, ослабить указанные недостатки.
Исследования подтвердили давно известный факт о том, что гравитационное поле заметно зависит от рельефа и может обоснованно трактоваться как стационарное только в равнинных районах. Но даже в равнинных районах возможность пользоваться изотропностью поля требует практического обоснования. В диссертации предложены подробные процедуры структурного анализа ГПЗ в виде построения разного рода ковариационных карт и вычисления ковариационных функций. Наиболее целесообразно, на наш взгляд, определять ковариационную функцию для исходных данных через двойное дискретное преобразование Фурье. Полученные таким способом ковариационные функции являются функциями разности координат двух точек, то есть функциями вектора, соединяющего эти две точки, и, следовательно, описывают неизотропное поле.
Если поле не является стационарным, то есть его спектральный состав изменяется с местоположением, то для изучения структуры поля приходится отказаться от фурье-анализа и естественно воспользоваться вейвлетанализом. Именно так и сделано в диссертации. Построение одномерных вейвлет-преобразований отдельных профилей поля и двумерных вейвлетспектрограмм позволяет выделить явно отличающиеся друг от друга участки поля и разделить, таким образом, изучаемую область на неперекрывающиеся подобласти, внутри каждой из которых поле можно считать стационарным.
Для каждой такой подобласти составляется матрица анизотропности, что позволяет построить ковариационную функцию для нестационарного поля в исходной области. Дальнейшие вычисления выполняются по стандартным формулам коллокации, но евклидовы расстояния между точками заменяются расстояниями Mahalanobis. Таким образом, известная идея о сегментации поля в данном случае реализуется непрерывно, в отличие от традиционного подхода, приводящего к проблемам нестыковки на границах подобластей.
Что касается ослабления практического недостатка коллокации, связанного с необходимостью решать большое количество уравнений, то вкладом диссертанта являюся алгоритмы двухэтапного гармонического анализа и определения функций влияния. Первый из них уменьшает количество уравнений, подлежащих решению, примерно вдвое, а второй алгоритм вообще не требует решения уравнений. Основные вычисления при этом выполняются в частотной области, а исходной информацией о гравитационном поле служит не ковариационная функция, а спектральная плотность.
Если поле стационарно, то достаточно найти только одну функцию влияния с помощью одного (обратного) преобразования Фурье, так как все остальные функции влияния для каждого узла kh можно получить аналитически или численно путём сдвига найденной материнской функции vh(x) по оси абсцисс на kh. Область с нестационарным полем трактуется как объединение подобластей с полем стационарным, и для каждой такой подобласти материнскую функцию влияния надо определять отдельно.
Все теоретические разработки автора обеспечены программами для персонального компьютера, позволяющими оперативно выполнять нужные вычисления и, при необходимости, сопровождать результаты графической интерпретацией.
Список публикаций по теме диссертации.
1. Попадьев В. В. Современное состояние метода коллокации. // Известия вузов. Геодезия и аэрофотосъемка, № 5, 2011. С. 10 15.
2. Басманов А. В., Попадьёв В. В., Сермягин Р. А. Развитие государственной гравиметрической сети Вьетнама // Геодезия и картография, № 5, 2011. С. 16 19.
3. Нейман Ю. М., Сугаипова Л. С., Попадьёв В. В. Эксперименты со спутниковой градиентометрией // Геодезия и картография, специальный выпуск, посвящённый конференции Научно-технические разработки в области геодезии и картографии и их применение в хозяйственной и оборонной деятельности страны, 2012.