CC BY
78
УДК 528.92
Н.М. Оскорбин, С. И. Суханов, Л.Ю. Федин Сглаживание массива данных рельефа
^ ^ * с использованием сплаиновых поверхностей
Ключевые слова: рельеф, топографический план, высотная отметка, GPS измерения, масштаб карты. Key words: relief, topographic plan, altitude mark, GPS measurements, map scale.
При представлении рельефа территории в процессе построения топографических карт и планов возникает проблема сглаживания исходных данных высотных отметок, полученных с использованием GPS измерений или известной базы данных НАСА. Необходимость выполнения указанной обработки данных состоит в том, что исходная информация содержит ошибки измерения, которые существенно искажают изолинии рельефа, а при малом разрешении неадекватно отображают реальную поверхность высот исследуемой территории. Проблема осложняется тем, что сглаживание необходимо проводить для массива данных большой размерности.
В представленной работе эта задача решается с использованием сплайновых поверхностей в среде электронных таблиц Excel.
Анализ литературы [1-4] показывает, что самым простым способом обработки трехмерных наборов данных является триангуляция, построенная на исходном массиве данных [4]. Полученная поверхность представляет собой систему плоских треугольников, «сшитых» значениями по координате г на линии граней и в узлах, в общем случае, нерегулярной сетки. Недостатком подхода являются, во-первых, отсутствие гладкости полученной поверхности и, во-вторых, он не решает проблему сглаживания, т.е. фильтрации хаотических ошибок измерения и не исключает выбросы, в том числе локально систематических.
Для построения сглаживающих поверхностей используют сплайны, в частности полный бикубический сплайн [2], взвешивающие полиномы высоких степей, ДММ сплайны, рестриктивные В-сплайны [3, 4]. В последнем методе построения сплайновых поверхностей допускаются ограниченные «скачки» производных соответствующих порядков. Заметим, что классические сплайны не содержат таких «скачков», т.е. при их построении обеспечивается «склеивание» в узлах сетки не только значений соседних полиномов, но и их первых производных. При выборе метода сглаживания необходимо учитывать главное требование к построенной поверхности и к его вычислительной сложности [2]:
* Работа выполнена при поддержке ведомственноаналитической программы «Развитие научного потенциала Высшей школы 2009-2010» (проект №2.2.2.4/4278).
во-первых, выбираемый класс функций должен описываться достаточно просто; во-вторых, поверхность не должна иметь особенностей, быть достаточно гладкой и иметь непрерывно изменяемую касательную или непрерывную кривизну; в-третьих, поиск нужной кривой в выделенном классе должен быть сравнительно легким, с вычислительной точки зрения. Следует отметить перспективный для развития геометрический подход к построению поверхностей [1], но его использование в качестве сглаживающих сплайнов требует специальных исследований.
Ориентируясь на возможности инструмента «Поиск решения» MS Excel, авторы применили неполный бикубический сплайн, который позволил на регулярной сетке обрабатывать массив высот размерностью (32 х 32). При построении сплайновой поверхности аналогично применению рестриктивных В-сплайнов [4] допускаются ограниченные и контролируемые «скачки» поверхностей по стыкам в координатах точек (Ху, у у) как по значениям, так
и по первой производной (по углам наклона рельефа). Используемый подход позволил разбить массив исходных данных на 16 классов, имеющих размеры (8 х 8) точек. Таким образом, приближение полигонов точек одного класса осуществляется по 64 точкам, что при числе свободных параметров, равном 11 для неполного бикубического сплайна, обеспечивает достаточное сглаживание в соответствии с рекомендациями работы [5, с. 70]. Опишем используемый подход.
Пусть (Ху, уу, гу) - координаты VGS-84 и отметки высот в метрах, i, j = 1,..., 32 . При обработке выделены полигоны (классы) (k,l), k,l = 1,...,4 и локальная для полигонов система координат точек (u, v) с прямоугольной на отметках (u, v) e D = = {-3.5;- 2.5;-1.5;...; 2.5; 3.5} сеткой. Приближение высот полигонов осуществляется следующей функцией (индексы k,l не приведены):
2 2 3 3
S(u, v) = X X a upvq + a30u +a03v , (1)
p=0q=0
где apq (p, q = 0,1,2); a03 , a30 - коэффициенты
в функции (1), которые необходимо определить из условий оптимального сглаживания.
Функция (1) имеет 11 свободных параметров и не учитывает пять членов полного бикубического сплайна.
Условие оптимального сглаживания отражает, прежде всего требование точности приближения
полигонов (/п , к,/ = 1,2,3,4). Кроме того, формализовано требование точности приближения стыков значений высот и углов наклона по линиям, разделяющим соседние полигоны, параллельных осям ох, оу . Соответствующие критерии обозначены как
1у/, Зух, Зку , 1 Гу . Для определения всех 176 коэффициентов функций (1) решается задача векторной оптимизации:
4 4,,,; 4
J = ЕЕ wkljn + Е
к =11 =1
d =0
е (jf + wjdx)+
+ i(dJy + w*j “)
. (2)
В формуле (2) индекс Л соответствует числу стыков по X и У, которые обозначаются соответственно (0 -1;1-2; 2- 3;3-4;4-0), Л = 0,1,2,3,4.
Задача (2) аналогична приближению функции обобщенным методом наименьших квадратов, в котором настраиваемые параметры типа Ш выступают в качестве весовых коэффициентов линейной свертки критериев.
Приведем выражения для частных критериев в задаче (2). Выражение для /Л имеет вид:
,, 1 8 8 г • -2
1п = Е I
64 т=1 н=11
т~\ к
где ит , Уи £ О ; 1т
ckll„kl ,kl\ kl
S ( (m , vn ) Zmn
(3)
отметка высоты соответст-
т и ’ ^ти
вующей координаты исходного массива.
При вычислении суммы квадратов значений и углов наклона стыков X и У начальная сетка удвоена и включает значения координат (и, у) как на сопрягаемых стыках, так и на внешней рамке рассматриваемого массива. Значение функций поверхностей и их углов наклона в клетках внешней рамки принято равными значениями прилегающих клеток. В том случае, когда обрабатываются несколько соседних массивов данных, эти значения могут служить дополнительным и эффективным средством согласования массивов поверхностей при независимом их построении.
Удвоенная сетка задается множеством Бр = {-4; -3,5; -3...;3; 3,5; 4} значений (и, у). Сопряжение поверхностей каждого полигона (к, /) осуществляется в точках у£,1 £ Ор ; ик е {-4,0; 4,0}
и иГ,1 е Ор ; Ур е {-4,0; 4,0}. При к < 4; / < 4 критерии согласования полигонов по значениям и углам наклона имеют вид (на примере ):
л = 1 17
17 т=1*-е £>
Jdl = — Е
т-
Л ,l ,.k+1,l
skl (,4.0)- sk+1,/ (, -4.0
u
, u„
(4)
где индекс d определяется по соответствующим значениям индекса k .
Аналогично выписываются критерии приближения соседних поверхностей по другим условиям. Если массив данных обрабатывается как изолированный, то значения критериев рассогласования «стыков» по клеткам рамки являются нулевыми, вычисление углов наклона осуществляется путем дифференцирования функции (1) в каждой точке расчетной сетки (u, v) e Dp ® Dp с последующим
преобразованием в градусную меру. При этом используется расстояние между узлами расширенной сетки в метрах (в нашем случае - 15 м).
Вычислительная технология сглаживания отметок рельефа и удвоения их числа включает следующие этапы:
1. Преобразование линейного массива точек (xt, yi, zi), i = 1,1024 в квадратный массив размера (32 х 32).
2. Получаемый массив точек (Ху, уу, Zj), i, у = 1,32 конвертируется для обработки с использованием элементов Ку, i, у = 1,32 матрицы К ,
состоящей из единиц. Эта матрица используется для корректировки отметок высот в случае обнаружения выбросов.
3. Производится визуализация исходных данных для задания начальных значений весовых коэффициентов критерия (2).
4. Формируются массивы полигонов, массив значений функции (1) на соседних стыках и массив углов наклона поверхности на полигонах и их границах.
5. Решается задача (2) определения оптимальных коэффициентов функции (1) и анализируется полученное решение.
6. Проводится поточечный анализ выбросов, по результатам которого корректируются значения коэффициентов функции (1) и исходных данных (путем изменения значений Ку матрицы К и весовых коэффициентов типа w).
7. Производится обратное конвертирование данных в линейный формат исходного массива и массива удвоенных точек.
Программная реализация метода сглаживания выполнена в двух книгах Excel. Первая книга служит для выполнения операций конвертирования исходных и расчетных данных. Обработка и визуализация данных выполняются во второй книге Excel, которая содержит 16 листов электронных таблиц. Выполняются расчетные операции по пунктам 3-6 рассмотренной вычислительной технологии.
Тестирование программы проведено при обработке массива в северо-западной части города Барнаула, на территории которого расположены участки с технологическими нарушениями, пахотные земли и природные ландшафты. Шаг сетки составляет 30 м, поэтому выбранный участок по размерам является квадратом со сторонами, равными 930 м.
Фрагменты расчетных таблиц и диаграмм приведены на рисунках 1-5. На рисунке 3 представлены исходные данные, анализ которых показывает необходимость выполнения процедур сглаживания. Макет расчетных таблиц на примере одного полигона ошибок сглаживания представлен на рисунке 1. Сверху и слева приведены координаты сетки. Заметим, что использование центрированной сетки, как показали исследования, значительно улучшает процесс приближения поверхности к исходным данным и их склеивания.
В С 0 Е Р С Н I Л
4 Квадр -3,5 -2,5 -1,5 -0,5 0,5 1,5 2,5 3,5
5 -3,5 0,22 0,78 -0,12 -0,35 0,21 -0,31 0,22 -0,08
6 -2,5 0,09 -0,33 -0,22 -0,44 0,12 0,595 0,11 -0,22
7 -1,5 0,01 -0,38 -0,25 -0,46 0,1 0,57 0,07 -0,27
8 -0,5 -0,04 0,59 -0,25 -0,46 0,12 0,584 0,08 -0,28
9 0,5 -0,11 0,56 -0,26 -0,45 0,13 0,605 0,09 -0,28
10 1,5 -0,23 0,48 -0,31 -0,47 0,12 0,602 0,09 -0,29
11 2,5 -0,41 0,33 -0,43 0,44 0,05 0,541 0,03 -0,35
12 3,5 0,29 0,08 0,36 0,25 -0,1 0,391 -0,11 -0,49
Рис. 1. Схема полигона и значения ошибок сглаживания (фрагмент)
Ж 7;:: Н І I ВШ к шш ншш N О І Р І О І Р I 8 I Т
Целевая функция Система ограничений по "стыку" значений X
79 1 2 3 4 Сумм 1 2 3 4
80 1 0,1182 0,28874 0,192402 0,163 0,763 Верх 0 0 0 0
81 2 0,1279 0,26532 0,170797 0,152 0,716 1-*2 0,0072 0.01026 0,0153 0.01963
82 3 0,239 0,34236 0,155969 0,125 0,862 2-*3 0,0004 0,00889 0,0232 0.01166
83 4 0,2837 0,16896 0.141441 0,236 0.83 3-М 0,0235 0.02042 0,0082 0,01237
84 нн Низ 0 0 0 0
85 Весовьк коэффии менты приближ ГК ЭЛИГОНО 0,02 0,16096
86 1 2 3 4 Вес приближ значений функции X
87 1 1 1 1 1 1 2 3 4
88 2 1 1 1 1 Верх 0 0 0 0
89 3 1 1 1 1 1_*2 1 1 1 1
90 4 1 1 1 1 2-*3 1 1 1 1
91 3-М 1 1 1 1
92 Приближение полигон-*-X ШШ Низ 0 0 0 0
93 I ]
94 Приближение полигон+У шш Система ого по "стыку" п роизводных X
95 I
96 Приближение ПОЛИГОН+Х+У шш 1 2 3 4
97 I ] Верх 0 0 0 0
98 Прибл полигон+Х+У+ПР X 1-*2 0,0025 0,00175 0,001 0,00199
99 2-*3 0,0111 0,01542 0,0067 0,00296
100 Прибл полигон+Х+У+ПР У ■■ 3-*4 0,0204 0,01163 0,0169 0,0098
101 I Низ 0 0 0 0
102 Прибл полигон+Х+У+ПР X, Ъ 0,10219
Рис. 2. Фрагмент таблиц вычисления критериев сглаживания и задания весовых коэффициентов
Рис. 3. Изображение исходных данных массива высот
Рис. 4. Изображение поверхности рельефа и данных массива высот (32x32)
Рис. 5. Изображение поверхности рельефа удвоенного массива данных высот (64 X 64)
На рисунке 2 представлена часть таблиц настройки весовых коэффициентов в задаче (2) и значения критериев. Вычислительный процесс по полному критерию задачи (2) является длительным при неконтролируемой погрешности. Рекомендуется расчеты проводить последовательно. Первоначально определяются коэффициенты неполного бикубического полинома по каждому из 16 полигонов раздельно. После оценки погрешностей и устранения выбросов следует выполнить операции приближенной сшивки по значениям и производным с контролем выбросов. Окончательное решение задачи (2) при хорошем начальном приближении для средней ЭВМ решается в пределах часа машинного времени. Контроль полученной поверхности можно провести визуально с использованием диаграмм, фрагменты которых представлены на рисунках 4 и 5, и анализом ошибок, максимальные и минимальные значения которых отображаются по полигонам в аналитических таблицах.
Построенная для тестовой территории поверхность рельефа приближает высоты исходного массива со средней по модулю погрешностью, равной 0,35 м, при размахе высот, равном 20 м. Эта оценка погрешности согласуется с условиями измерений, так как отсчеты рельефа представлены целыми числами (см. рис. 3).
Полученная в математической форме поверхность рельефа и сглаженные исходный и удвоенный массивы данных могут быть использованы при выполнении геодезических, картографических и проектных работ. Заметим, что использование предложенного подхода как метода математического моделирования пространственных «полей» и рельефа местности требует специальных исследований.
Библиографический список
1. Препората, Ф. Вычислительная геометрия / Ф. Пре-пората, М. Шеймос. - М., 1989.
2. Шикин, Е.В. Кривые и поверхности на экране компьютера: руководство по сплайнам для пользователей / Е.В. Шикин, Л.И. Плис. - М., 1996.
3. Поддубный, В.В. Методы инвариантного погружения и аппроксимации в рестриктивных задачах управления и фильтрации / В.В. Поддубный. - Томск, 1993.
4. Поддубный, В.В. Приближение сложных поверхностей с использованием аппарата рестриктивных В-сплай-нов / В.В. Поддубный, М.В. Черноусов // Геоинформатика. - Томск, 1999. - Вып. 1.
5. Максимов, А.В. Многопользовательские информационные системы: основы теории и методы исследования: монография / А.В. Максимов, Н.М. Оскорбин. - Барнаул, 2005.
