Сплайн-аппроксимация на основе триангуляции
Описан алгоритм полиномиальной сплайн-аппроксимации на основе триангуляции, используемый при построении карт. Метод обеспечивает более точное совпадение с заданными точками, чем B-сплайны. An algorithm of the polynomial spline approximation based on triangulation is described. The algorithm is used...
Saved in:
| Date: | 2008 |
|---|---|
| Main Author: | |
| Format: | Article |
| Language: | Russian |
| Published: |
Інститут програмних систем НАН України
2008
|
| Subjects: | |
| Online Access: | https://nasplib.isofts.kiev.ua/handle/123456789/1500 |
| Tags: |
Add Tag
No Tags, Be the first to tag this record!
|
| Journal Title: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| Cite this: | Сплайн-аппроксимация на основе триангуляции / Е.В. Назаренко // Пробл. програмув. — 2008. — N 2-3. — С. 657-664. — Бібліогр.: 11 назв. — рус. |
Institution
Digital Library of Periodicals of National Academy of Sciences of Ukraine| id |
nasplib_isofts_kiev_ua-123456789-1500 |
|---|---|
| record_format |
dspace |
| spelling |
Назаренко, Е.В. 2008-07-31T15:10:30Z 2008-07-31T15:10:30Z 2008 Сплайн-аппроксимация на основе триангуляции / Е.В. Назаренко // Пробл. програмув. — 2008. — N 2-3. — С. 657-664. — Бібліогр.: 11 назв. — рус. 1727-4907 https://nasplib.isofts.kiev.ua/handle/123456789/1500 680.3.06 Описан алгоритм полиномиальной сплайн-аппроксимации на основе триангуляции, используемый при построении карт. Метод обеспечивает более точное совпадение с заданными точками, чем B-сплайны. An algorithm of the polynomial spline approximation based on triangulation is described. The algorithm is used for building maps. The method provides more precise fitting to defined points than B-splines. ru Інститут програмних систем НАН України №2-3 С. 657-664 Прикладне програмне забезпечення Сплайн-аппроксимация на основе триангуляции Spline-approximation on the base of triangulation Article published earlier |
| institution |
Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| collection |
DSpace DC |
| title |
Сплайн-аппроксимация на основе триангуляции |
| spellingShingle |
Сплайн-аппроксимация на основе триангуляции Назаренко, Е.В. Прикладне програмне забезпечення |
| title_short |
Сплайн-аппроксимация на основе триангуляции |
| title_full |
Сплайн-аппроксимация на основе триангуляции |
| title_fullStr |
Сплайн-аппроксимация на основе триангуляции |
| title_full_unstemmed |
Сплайн-аппроксимация на основе триангуляции |
| title_sort |
сплайн-аппроксимация на основе триангуляции |
| author |
Назаренко, Е.В. |
| author_facet |
Назаренко, Е.В. |
| topic |
Прикладне програмне забезпечення |
| topic_facet |
Прикладне програмне забезпечення |
| publishDate |
2008 |
| language |
Russian |
| publisher |
Інститут програмних систем НАН України |
| format |
Article |
| title_alt |
Spline-approximation on the base of triangulation |
| description |
Описан алгоритм полиномиальной сплайн-аппроксимации на основе триангуляции, используемый при построении карт. Метод обеспечивает более точное совпадение с заданными точками, чем B-сплайны.
An algorithm of the polynomial spline approximation based on triangulation is described. The algorithm is used for building maps. The method provides more precise fitting to defined points than B-splines.
|
| issn |
1727-4907 |
| url |
https://nasplib.isofts.kiev.ua/handle/123456789/1500 |
| citation_txt |
Сплайн-аппроксимация на основе триангуляции / Е.В. Назаренко // Пробл. програмув. — 2008. — N 2-3. — С. 657-664. — Бібліогр.: 11 назв. — рус. |
| work_keys_str_mv |
AT nazarenkoev splainapproksimaciânaosnovetriangulâcii AT nazarenkoev splineapproximationonthebaseoftriangulation |
| first_indexed |
2025-11-26T23:47:27Z |
| last_indexed |
2025-11-26T23:47:27Z |
| _version_ |
1850783280808329216 |
| fulltext |
Прикладне програмне забезпечення
© Е.В. Назаренко, 2008
ISSN 1727-4907. Проблеми програмування. 2008. № 2-3. Спеціальний випуск 657
УДК 680.3.06
СПЛАЙН-АППРОКСИМАЦИЯ НА ОСНОВЕ ТРИАНГУЛЯЦИИ
Е.В. Назаренко
Институт кибернетики им. В.М. Глушкова НАН Украины,
03028, Киев, проспект Науки 43, кв. 70.
Тел.: 524 5330, eugn@ukr.net
Описан алгоритм полиномиальной сплайн-аппроксимации на основе триангуляции, используемый при построении карт. Метод
обеспечивает более точное совпадение с заданными точками, чем B-сплайны.
An algorithm of the polynomial spline approximation based on triangulation is described. The algorithm is used for building maps. The method
provides more precise fitting to defined points than B-splines.
Введение
Задача построения карт хорошо изучена и имеет широкое применение. Разработан ряд методов ее
решения: триангуляция, метод минимальной кривизны [1, 2], методы кригинга [3, 4], методы сплайн-
аппроксимации [5], метод Шепарда, скользящей плоскости, инверсных расстояний, N ближайших соседей,
натуральной окрестности, радиальных базисных функций и др.
Методы сплайн-аппроксимации имеют определенные преимущества: возможность использования
стабилизаторов для задания общих свойств картируемой поверхности, аддитивное включение в функционал
произвольного числа качественной и количественной информации с использованием весовых коэффициентов в
качестве управляющих параметров [6], устойчивость и надежность.
Недостатки сплайн-аппроксимации [5]:
- значения на карте могут выходить за диапазон известных значений в области картопостроения;
- вычислительная сложность;
- слабый по сравнению с кригингом контроль ошибки аппроксимации;
- кластерный эффект (приводит к тому, что точная аппроксимация внутри плотных скоплений известных
точек компенсирует большие отклонения в отдельно стоящих точках и на краях области аппроксимации).
Алгоритм полиномиальной сплайн-аппроксимации на основе триангуляции обеспечивает лучшую для
вычислений систему уравнений, чем аппроксимация на основе радиальных базисных функций, а также
отсутствие «эффекта бычьего глаза», характерного для большинства других методов.
Построение триангуляции
Диаграмма Вороного. Пусть S – конечное множество точек на плоскости. Многоугольником Вороного
)(* pV для точки Sp ∈ называется множество точек на плоскости, более близких к p , чем к любой другой
точке множества S [7]. Пусть ip и jp – две точки. Тогда множеством точек, более близких к ip , чем к jp ,
является содержащая точку ip полуплоскость ),( ji ppH , ограниченная прямой, проходящей через середину
отрезка ],[ ji pp перпендикулярно к нему. I
ji
jii ppHpV
≠
= ),()(* [8]. Границу многоугольника )(*
ipV будем
обозначать как )( ipV . Для заданного множества точек S диаграмма Вороного состоит из множества
многоугольников Вороного для каждой из точек множества S (рис. 1). Ребро диаграммы Вороного, лежащее
между двумя многоугольниками Вороного )(*
ipV и )(*
jpV , состоит из точек на плоскости, эквидистантных
точкам ip и jp , и расположено ближе к этим точкам, чем к остальным точкам множества S [7]. Вершины
диаграммы Вороного определяются вершинами многоугольников Вороного.
Рассмотрим алгоритм построения диаграммы Вороного, приведенный в [8]. Определим:
aP – массив точек множества S ;
lEdg – список ребер диаграммы;
aVor – массив многоугольников Вороного.
Список lEdg состоит из кортежей VOR_EDGE >∈∈∈∈< }2,1,0{,,, 21 oEeNnNn . E – множество отрезков
вида >∈∈< 22 , RdestRorg . Элементы 1n , 2n кортежа задают соответственно левую ( lp ) и правую ( rp )
относительно вектора e (направленного из org в dest) точки множества S такие, что )()( rl pVpVe ∩∈ . Ребро
e , в свою очередь, является ребром диаграммы Вороного. Ребра бывают трех типов: прямая, луч и отрезок. Им
отвечают значения 0, 1, 2 элемента o . Ребро является прямой тогда и только тогда, когда все точки диаграммы
Прикладне програмне забезпечення
658
лежат на одной прямой. Лучом ребро является лишь в случае, когда ))()(( 21 SCHnSCHn ∈⇔∈ . Обратное тоже
верно. Этот факт нам понадобится при вычислении SSSCH ⊂**),( . Отрезком ребро является во всех
остальных случаях.
Каждой точке aP[i] отвечает i-й элемент массива aVor,
который содержит многоугольник Вороного для данной
точки. Многоугольник aVor [i] организован в виде списка
ребер. Ребро имеет следующую структуру
(EDGE_FOR_NODE): >∈∈∈< ** ,, FpNENnEpE , где
*E – множество указателей на кортежи типа VOR_EDGE,
*F – множество указателей на кортежи типа
EDGE_FOR_NODE. Элемент pE кортежа задает
собственно ребро. n задает индекс соседа по ребру, на
которое указывает pE , в массиве aP. Элемент pNE
указывает на то же самое ребро у соседа. Взаимосвязь
структур данных показана на рис. 2. Отметим, что ребра в
многоугольнике отсортированы по часовой стрелке, что иллюстрирует, в частности, рис. 2.
Рис. 2. Представление диаграммы Вороного
Метод “разделяй и властвуй” состоит в следующем. Задача делится на несколько подзадач, которые
решаются рекурсивно и затем достигнутые результаты объединяются для получения решения исходной задачи.
По своему виду подзадачи сходны с исходной задачей, но они меньше по размеру и, естественно, сумма
размеров подзадач равна размеру исходной задачи. Если задача достаточно мала или проста, она может быть
решена непосредственно и нет необходимости ее дальнейшего деления – это условие является критерием для
окончания процесса деления [7].
В нашем случае задача каждый раз делится на две подзадачи, а условием окончания является случай, когда
на вход рекурсивной функции подана одна точка.
Алгоритм построения диаграммы Вороного:
1) sort( S );
2) удалить одинаковые точки;
3) если S > 1 VorD(0, S –1).
Функция N)lPointN,fPoint( ∈∈VorD имеет следующий вид:
если fPoint ≠ lPoint то
+
2
lPointfPoint
fPoint,VorD ;
+
+
lPoint,1
2
lPointfPoint
VorD ;
Рис. 1. Диаграмма Вороного
Прикладне програмне забезпечення
659
Объединить )( 1SVor и )( 2SVor , { }U
+
=
=
2
lPointfPoint
fPoint
1 ][
i
iaPS , { }U
lPoint
1
2
lPointfPoint
2 ][
+
+=
=
i
iaPS .
Первым этапом решения задачи является сортировка исходного множества точек S по координате x в
порядке возрастания. В случае совпадения значений, точки сортируются по координате y (тоже в порядке
возрастания).
Для разделения задачи на подзадачи в нашем случае используется функция VorD. В общих чертах
алгоритм объединения двух диаграмм описан далее.
Пусть два множества точек 1S и 2S ( SSS ⊂21, ) линейно разделимы вертикальной прямой const=x .
Пусть также имеем диаграммы Вороного )( 1SVor и )( 2SVor (рис. 1). Наша задача состоит в том, чтобы путем
объединения данных диаграмм за линейное время получить диаграмму Вороного для множества 212,1 SSS ∪= .
Введем понятие разделяющей цепи. Для заданного разбиения { }21, SS множества 2,1S разделяющая цепь
),( 21 SSσ обозначает множество ребер диаграммы Вороного, общих для пар многоугольников )( ipV и )( jpV
диаграммы )( 2,1SVor , где 1Spi ∈ и 2Sp j ∈ [8]. На рис. 1 изображена разделяющая цепь между множествами 1S
и 2S . Поскольку данные множества линейно разделимы вертикальной прямой, то σ однозначно разбивает
плоскость на левую lπ и правую rπ части. Пусть 1S лежит слева от 2S . Тогда диаграмма Вороного )( 2,1SVor
представляет собой объединение lSVor π∩)( 1 , rSVor π∩)( 2 и σ .
Получить σ можно следующим образом. Пусть ),( 21 eeI означает пересечение линий 1e и 2e , каждая из
которых – прямая, луч или отрезок. 21,vv – расстояние между точками 1v и 2v . λ – отсутствие пересечения.
Пусть также ось абсцисс направлена вправо, ось ординат – вниз.
1. Построить верхний опорный отрезок [ ]][],[ T
R
T
L paPpaP ; построить нижний опорный отрезок
[ ]][],[ B
R
B
L paPpaP ; lp = T
Lp ; rp = B
Rp ; opnd = 0; σ = 0;
2. l = прямая, перпендикулярная [ ]][],[ T
R
T
L paPpaP и делящая [ ]][],[ T
R
T
L paPpaP пополам; v = точка из l,
бесконечно удаленная вправо, если l – горизонтальная и точка из l, бесконечно удаленная вверх иначе;
3. e = прямая, перпендикулярная [ ]][],[ T
R
T
L paPpaP и делящая [ ]][],[ T
R
T
L paPpaP пополам;
4. *
lp = { ),( meI | )( lpVm ∈ и ),(, meIv – минимально, если opnd = 0 и vmeI ≠),( иначе};
*
rp = { ),( meI | )( rpVm ∈ и ),(, meIv – минимально, если opnd = 0 и vmeI ≠),( иначе };
5. если *
lp = λ и *
rp = λ перейти на 8;
если *
lp = λ перейти на 6.3;
если *
rp = λ перейти на 6.2;
6. если *, lpv = *, rpv
6.1. σ = { }],( *
lpv∪σ , если opnd = 0 и { }],[ *
lpv∪σ иначе;
v = *
lp ; opnd = opnd + 1;
lp = сосед по ребру *e (см. рис. 3.1, 3.2);
rp = сосед по ребру **e (см. рис. 3.1, рис. 3.2);
иначе если *, lpv < *, rpv
6.2. σ = { }],( *
lpv∪σ , если opnd = 0 и { }],( *
lpv∪σ иначе;
v = *
lp ; opnd = opnd + 1;
lp = сосед по ребру *e (см. рис. 3.2);
иначе
6.3. σ = { }],( *
rpv∪σ , если opnd = 0 и { }],( *
rpv∪σ иначе;
v = *
rp ; opnd = opnd + 1;
rp = сосед по ребру **e (см. рис. 3.2);
7. если lp ≠ B
Lp ∨ rp ≠ B
Rp перейти на 3;
Прикладне програмне забезпечення
660
8. e = прямая, перпендикулярная [ ]][],[ rl paPpaP и делящая [ ]][],[ rl paPpaP пополам; *v = точка из e,
бесконечно удаленная влево, если e – горизонтальная и точка из e, бесконечно удаленная вниз иначе;
σ = { }),( *vv∪σ , если opnd = 0 и { }),[ *vv∪σ иначе.
а
б
Рис. 3
Алгоритм построения верхнего и нижнего опорного отрезков подробно описан в [7]. Заметим только, что
он использует выпуклые оболочки множеств 1S и 2S , которые могут быть эффективно вычислены следующим
образом. Ранее мы установили связь между неограниченными многоугольниками и выпуклой оболочкой. Ребра
в неограниченном многоугольнике отсортированы таким образом, что лучами являются первое и последнее
ребро в списке ребер многоугольника. Для простоты предположим, что )( 2,1SVor не содержит прямых.
Следующая последовательность приводит к вычислению )( 2,1SCH :
1. если S = 1
A = {p}; //aP[p] – точка, принадлежащая выпуклой оболочке, например, точка с наибольшей
координатой x
конец вычисления;
2. A = ∅ ;
ip = p;
3. A = { }ipA ∪ ;
4. ip = элемент n первого элемента aVor[ ip ];
5. если ip ≠ p перейти на 3.
Таким образом, A – искомая выпуклая оболочка множества 2,1S . В приведенном алгоритме она
вычисляется в порядке обхода по часовой стрелке. )( 2,1SCH , упорядоченную против часовой стрелки, можно
вычислить, заменив шаг 4 шагом *4 :
*4 . ip = элемент n последнего элемента aVor[ ip ].
Случай, когда )( 2,1SVor состоит из прямых, требует внесения лишь незначительных уточнений в приведенный
алгоритм.
Шаг 4. алгоритма построения разделяющей цепи имеет следующую особенность. Цепь σ , оставаясь
внутри некоторого многоугольника )( ipV , может многократно менять свое направление, прежде чем она
пересечет какое-либо ребро многоугольника )( ipV . Однако, ввиду особенностей алгоритма, при поиске
пересечений (начиная со второго) следует рассматривать только ребра, идущие по часовой стрелке после
предыдущего пересечения для левой диаграммы и против часовой стрелки для правой диаграммы.
Чтобы завершить слияние диаграмм, остается удалить все ребра (и части ребер) диаграммы )( 2SVor ,
расположенные слева от σ , и все ребра (и части ребер) )( 1SVor , расположенные справа от σ . Обрезать и
удалять ребра можно одновременно с построением разделяющей цепи.
Триангуляция Делоне. Задача триангуляции множества точек может быть сформулирована следующим
образом. На плоскости заданы N точек. Соединить их непересекающимися отрезками так, чтобы каждая
область внутри выпуклой оболочки этого множества точек являлась треугольником [8].
Выпуклой триангуляцией называется такая триангуляция, для которой минимальный многоугольник,
охватывающий все треугольники, будет выпуклым. Говорят, что триангуляция удовлетворяет условию Делоне,
Прикладне програмне забезпечення
661
если внутрь окружности, описанной вокруг любого построенного треугольника, не попадает ни одна из
заданных точек триангуляции. Триангуляция называется триангуляцией Делоне, если она является выпуклой и
удовлетворяет условию Делоне [9].
Пусть S – множество точек на плоскости, причем все точки множества S не лежат на одной прямой. Пусть
также никакие четыре точки исходного множества не лежат на одной окружности. Тогда имеет место
утверждение.
Теорема Делоне. Граф, двойственный диаграмме Вороного (рис. 4), является триангуляцией (Делоне)
множества S [8].
Рис. 4. Диаграмма Вороного и двойственная ей триангуляция Делоне
Таким образом, триангуляция Делоне может быть получена из диаграммы Вороного. В случае, когда
четыре точки множества S лежат на одной окружности, после “обращения” ребер графа диаграммы, помимо
треугольников, образуются также многоугольники. Для получения триангуляции они должны быть разбиты на
треугольники произвольным образом.
Заметим, что вершины диаграммы Вороного являются центрами описанных окружностей вокруг
треугольников триангуляции Делоне того же множества, для которого построена диаграмма. Это свойство
может быть использовано для получения триангуляции Делоне из диаграммы Вороного:
1. Для каждого ребра, если оно не является лучом, сделать дубликат.
2. Отсортировать ребра по org (или по dest).
3. Сформировать триангуляцию.
В результате сортировки ребра, инцидентные одной вершине, будут расположены подряд. Располагая
информацией о точках множества S , между которыми построено ребро (элементы 1n , 2n кортежа
VOR_EDGE), триангуляция Делоне может быть построена за один проход по отсортированному списку ребер
диаграммы.
Рассмотрим шаг 2. алгоритма подробнее. Будем говорить, что 1e больше, чем 2e если функция ),( 21 eeless
возвращает true и 1e меньше, чем 2e – иначе. Функция ),( 21 eeless имеет следующий вид:
1. 1v = вершина, которой инцидентно ребро 1e ; 2v = вершина, которой инцидентно ребро 2e ;
2. если xv .1 = xv .2
если yv .1 = yv .2
α = угол, который образуется при переходе вектора (0,-1) в вектор 1e при вращении вектора
(0,-1) вокруг точки 1v по часовой стрелке;
β = угол, который образуется при переходе вектора (0,-1) в вектор 2e при вращении вектора
(0,-1) вокруг точки 2v по часовой стрелке;
если α < β , то less = true ; иначе less = false ;
иначе
если yv .1 < yv .2 , то less = true ; иначе less = false ;
3. если xv .1 < xv .2 , то less = true ;
4. less = false ;
Примечание. “.” – элемент кортежа.
Для задания триангуляции используем массив, состоящий из кортежей >∈∈∈< NpNpNp 321 ,, . Элементы
кортежа задают индексы вершин треугольника в массиве aP.
Прикладне програмне забезпечення
662
Построение карты
В общем случае задача построения карты может быть сформулирована так. На двумерной ограниченной
области аппроксимации по известной количественной и качественной информации необходимо отыскать
неизвестную функцию f , принадлежащую пространству H функций, непрерывных вместе с первыми
производными на всей области аппроксимации, за исключением заданного множества меры нуль [10]. Однако
на практике для представления карт используются значения f в узлах заранее заданной регулярной
прямоугольной сетки значений (заданной прямоугольником r и шагами x∆ , y∆ ).
В нашем случае количественная информация представлена известными значениями в некоторых точках и
описывается функционалом вида
( )2
),()( ∑ −=
i
iiiz zyxffJ .
Качественная информация чаще всего бывает следующих трех видов:
1. Искомая функция f напоминает другую известную функцию af – например, карту вышележащего
горизонта, крупномасштабную карту или карту, полученную на основании других способов измерения.
∑
∂
∂
−
∂
∂
+
∂
∂
−
∂
∂
=
ji
jiajijiaji
aa y
yxf
y
yxf
x
yxf
x
yxf
ffJ
,
22
),(),(),(),(
),( .
2. Искомая функция f должна изменяться плавно, с минимальными искривлениями. Физическая модель
плавно изменяющейся функции – плоская упругая пластина, искривляемая силами притяжения к точкам, чье
положение в пространстве определяется известными значениями:
∑
∂
∂
+
∂
∂
=
ji
jiji
y
yxf
x
yxf
fJ
,
2
2
22
2
2
2
),(),(
)( .
3. Искомая функция f соответствует стационарному процессу, ожидаемое значение в случайной точке
не зависит от горизонтального смещения:
∑
∂
∂
+
∂
∂
=
ji
jiji
y
yxf
x
yxf
fJ
,
22
1
),(),(
)( .
Таким образом подавляются всплески значений на краях области построения карты, где отсутствуют
известные значения.
С учетом вышеизложенного задача картопостроения сводится к минимизации функционала:
min))1(( 12
* →−+++=Φ JJJJ saazz γγλλλ .
Далее вместо функционала *Φ рассмотрим функционал (поскольку конечно-разностное представление aJ
аналогично конечно-разностному представлению 1J ) ))1(( 12 JJJ szz γγλλ −++=Φ .
Пусть S – множество точек на плоскости, в которых известно значение искомой функции f . Область
аппроксимации ограничим прямоугольником { }),),( maxminminmin yyyxxxyxr ≤≤≤≤= . Пусть
{ }),(),,(),,(),,( maxmaxminmaxmaxminminmin
* yxyxyxyxSS ∪= . Область аппроксимации разобъем сеткой триангуляции множества *S на
множество непересекающихся подобластей, каждая из которых представляет треугольник с вершинами в
точках множества *S . Тогда пространство H аппроксимируется конечномерным пространством Hs функций,
совпадающих на каждой подобласти с полиномом Лагранжа третьей степени:
yxbxybybxbxybybxbybxbbf L
2
9
2
8
3
7
3
65
2
4
2
3210 +++++++++=
и удовлетворяющих условиям непрерывности функции/производных на границах подобластей [5].
Для задачи восстановления f в узлах регулярной сетки производную целесообразно считать методом
конечных разностей:
( ) ( ) ( ) ( )
∑ ∑∑∑
= = =
−
=
−
∆
−
+
∆
−
=
x x yyn
i
n
i
n
j y
jiji
n
j x
jiji yxfyxfyxfyxf
fJ
2 1 2
2
1
1
2
1
1
,,,,
)( .
Прикладне програмне забезпечення
663
Точки, между которыми проходит разлом, не должны быть включены в 1J .
( ) ( ) ( ) ( ) ( ) ( ) 2
1 3
2
21
3 1
2
2
21
2
,,2,,,2,
)( ∑∑∑∑
= =
−−
= =
−−
∆
+−
+
∆
+−
=
x yx y n
i
n
j y
jijiji
n
i
n
j x
jijiji yxfyxfyxfyxfyxfyxf
fJ .
Точки, между которыми проходит разлом, не должны быть включены в 2J .
Для минимизации функционала Φ составляется система линейных уравнений
( )
0
,...,1 =
∂
Φ∂
k
n
a
aa
, nk ,1= .
Введем следующие обозначения:
yxaxyayaxaxyayaxayaxaaf iiiiiiiiii
i
L
2
910
2
810
3
710
3
610510
2
410
2
31021011010 +++++++++ +++++++++= ;
),( yxck – коэффициент функции ( )91010 ,..., +ii
i
L aaf при переменной ka (т.е.,
k
i
L
k a
f
yxc
∂
∂=),( ),
=
10
k
i ;
( ) ( )
+≤≤
=
иначе,0
91010,,
,
ikiеслиyxc
yxc ki
k ;
( )yxl , – индекс функции, построенной на треугольнике, содержащем точку ( )yx, .
Тогда
( ) +
−
λ=
∂
Φ∂
∑
=
3
1
1010101010101 ,,
,...,
2
1
j
k
j
k
j
k
j
k
L
k
j
k
jkz
k
n zyxfyxc
a
aa
( ) ( ) ( )( ) ( ) ( )( ) +
−−
∆
γ−λ+ ∑∑
= =
−−
−−
x y
jijijiji
n
i
n
j
ji
yxl
Lji
yxl
Lji
yxl
kji
yxl
k
x
s yxfyxfyxcyxc
2 1
1
),(),(
1
),(),(
2
,,,,
1
1 11
( ) ( ) ( )( ) ( ) ( )( ) +
−−
∆
γ−λ+ ∑∑
= =
−−
−−
x y
jijijiji
n
i
n
j
ji
yxl
Lji
yxl
Lji
yxl
kji
yxl
k
y
s yxfyxfyxcyxc
1 2
1
),(),(
1
),(),(
2
,,,,
1
1 11 (1)
( ) ( ) ( )( ) ( ) ( ) ( )( ))++−
+−
∆
γλ+ −−
= =
−−
−−−−∑∑ ji
yxl
Lji
yxl
Lji
yxl
L
n
i
n
j
ji
yxl
kji
yxl
kji
yxl
k
x
s yxfyxfyxfyxcyxcyxc jijiji
x y
jijiji ,,2,,,2,
1
2
),(
1
),(),(
3 1
2
),(
1
),(),(
4
2121
( ) ( ) ( )( ) ( ) ( ) ( )( ))2
),(
1
),(),(
1 3
2
),(
1
),(),(
4
,,2,,,2,
1
2121
−−
= =
−−
−−−− +−
+−
∆
γλ+ ∑∑ ji
yxl
Lji
yxl
Lji
yxl
L
n
i
n
j
ji
yxl
kji
yxl
kji
yxl
k
y
s yxfyxfyxfyxcyxcyxc jijiji
x y
jijiji .
Отметим, что ( )yxc yxl
k ,),( для большинства x и y равняется нулю, так, что при составлении (1) при
суммировании можно учитывать только те точки, которые отстоят от треугольника, которому соответствует
функция
10
k
Lf , не более чем на x∆2 и y∆2 соответственно.
Далее эта система решается, что приводит к вычислению неизвестных коэффициентов naa ,...1 и, таким
образом, к определению функции f . Для решения задачи остается лишь вычислить значения уже известной
функции f в узлах сетки. Заметим, что во избежание переполнений при составлении системы линейных
уравнений, значения x и y заменяются соответственными значениями с интервалов [ ]maxmin , xx и [ ]maxmin , yy ,
причем должно выполнятся равенство
*
*
y
x
y
x
∆
∆=
∆
∆
, где
1
minmax*
−
−
=∆
x
x n
xx
,
1
minmax*
−
−
=∆
y
y n
yy
. Интервалы могут быть
выбраны, например, такими: [ ]1,0 −∆ x , ( )
−
∆
∆
1,0 y
x
y n .
Реализация метода в пакете «ГеоПоиск»
Метод включен как один из методов картопостроения в пакет «ГеоПоиск» [11].
Рассмотрим применение метода для получения кровли пласта по имеющимся скважинам (рис. 5).
Положение скважин определяет разбиение области аппроксимации (рис. 6). Каждому треугольнику отвечает
функция – полином Лагранжа третьей степени, которая представлена в системе линейных уравнений десятью
коэффициентами. Таким образом, количество скважин определяет порядок системы.
Прикладне програмне забезпечення
664
Рис. 5
Рис. 6
Для расчета поверхности необходимо задать следующие параметры: размер сетки, вес известных значений,
вес поверхности-аналога, вес стабилизатора, доля 1-й производной в стабилизаторе. Размер сетки определяет
количество узлов, в которых будет вычислена поверхность. Кроме того, более густая сетка обеспечивает более
точную аппроксимацию поверхности в известных точках (за счет увеличения точности вычисления
производных), однако требует больше времени на составление системы линейных уравнений. Веса позволяют
управлять влиянием соответствующего фактора (упругость, поверхность-аналог и т.п.) на поверхность. На
рис. 7 показано вычисленную приведенным
методом карту.
Заключение
Рассмотренный метод обеспечивает
построение достаточно точной реалистичной
карты. Его эффективная реализация включена в
набор методов картопостроения пакета
«ГеоПоиск», предназначенного для обработки и
интерпретации данных геофизического
исследования нефте-газовых скважин с
привлечением смежной геологической и
геофизической информации. Разломы (нарушения)
с известной амплитудой, проходящие через
кровлю/подошву пласта, реализованы путем
предварительного вычитания из известных
значений вертикального сдвига в
соответствующих точках и последующего
прибавления к карте поверхности карты сдвига из-
за разломов. Усовершенствовать реализацию
разломов с неизвестной амплитудой можно, отдельно рассматривая части треугольника, разделенных
разломом. Также планируется реализовать правку сетки триангуляции вручную, так как триангуляция Делоне,
оптимизированная по форме треугольников, не отвечает требованию выстраивания длинного ребра
треугольника вдоль направления напластования пород, в соответствии геологическими представлениями
специалистов.
1. Smith W. H. F., Wessel P. Gridding with continuous curvature splines in tension // Geophysics. – 1990. – Vol. 55, N 3. – P. 293 – 305.
2. Красавчиков В.О. Математические методы, алгоритмы и технология геомоделирования для решения задач геологии нефти и газа.
Автореф… докт. дис. – Кемерово. – 2007.
3. Abdullah Arik. Area Influence Kriging // Mathematical Geology. – 2002. – Vol. 34, N 7. – P. 783 – 796.
4. Schenkova Z., Schenk V., Kalogeras I., Pichl R., Kottnauer P., Papatsimba C., Panopoulou G. Isoseismal maps drawing by the kriging method //
Journal of Seismology. – 2007. – Vol. 11, N 3. – P. 345 – 353.
5. Тульчинский В.Г. Технология картопостроения // Кибернетика и системный анализ. – 1999. – № 3. – С. 17 – 21.
6. Плавник А.Г. Алгоритмизация геоинформационных технологий в задачах, связанных с картопостроением. Автореф. канд. дис. –
Тюмень. – 2004.
7. Майкл Ласло. Вычислительная геометрия и компьютерная графика на C++. – 1997. – 304 с.
8. Препарата Ф., Шеймос М. Вычислительная геометрия: Введение. – 1989. – 478 с.
9. Скворцов А.В. Триангуляция Делоне и её применение. – 2002. – 128 с.
10. Стечкин С.В., Субботин Ю.Н. Сплайны в вычислительной математике. – М.: Наука, 1976. – 248 с.
11. www.geopoisk.com.
Рис. 7
|