CC BY
27
УДК 528.2, 004:528
ПРИМЕНЕНИЕ ТЕОРИИ ВЕЙВЛЕТОВ ПРИ СЖАТИИ И ФИЛЬТРАЦИИ ГЕОИНФОРМАЦИИ
А.С.ЯРМОЛЕНКО1, О.В.СКОБЕНКО2
1 Новгородский государственный университет им. Ярослава Мудрого, Великий Новгород, Россия
2 Белорусская государственная сельскохозяйственная академия, г. Горки, Могилевская область, Республика Беларусь
Целью статьи является разработка детальной и доступной технологии применения вейвлетов в обработке геоинформации, предметом исследований - вейвлет-методы фильтрации и сжатия геоинформации. Методология исследований базируется на современной теории вейвлетов в свете линейной алгебры. Методы исследований - изучение и обобщение, абстрагирование, формализация, математическое моделирование с помощью составленных авторами программ для ЭВМ.
После введения и постановки задачи приведены базовые положения линейной алгебры, на которых основывается содержание статьи при построении ортонормированных базисов в одно- и двумерном случаях. Вначале приводится применение общей теории к разложению вектора исходных данных в базисах Хаара и Шеннона. Далее на основе базиса Хаара строятся ортонормированные базисы вейвлет-преобразований и фильтрации информации. Рассмотрен порядок создания вейвлет-фильтров последовательностью сверток, применение КМА-анализа для построения ортонормированного базиса вейвлет-преобразования. Реализована практическая возможность вейвлет-фильтрации на основе составленных конкретных программ моделирования полей данных геоинформации и изображений, сжатия данных и их фильтрации. Результатом работы являются методики построения ортонормированных базисов различными методами вейвлет-преобразования, на основе которых составлены алгоритмы и соответствующие программы для ЭВМ по сжатию геоинформации на примере рельефа местности и фотоизображений. Исследована эффективность сжатия геоинформации и фильтрации шумов с помощью вейвлетов. Разработана методика определения значения величины фильтра в зависимости от точности исходной геоинформации, проиллюстрированная на примере расчета значения фильтра для сжатия информации о высотах рельефа местности. Такая же методика рекомендована и для фильтрации изображений.
Ключевые слова: базис вейвлет-преобразования; свертка; КМА-анализ; фильтрация; сжатие; точность; программа; изображение; поле данных
Как цитировать эту статью: Ярмоленко А.С. Применение теории вейвлетов при сжатии и фильтрации геоинформации / А.С.Ярмоленко, О.В.Скобенко // Записки Горного института. 2018. Т.234. С. 612-623. DOI: 10.31897/PMI.2018.6.612
Введение и постановка задачи. Сжатие геоинформации и ее фильтрация являются актуальной проблемой теории математической обработки как геодезической информации (измерений), так и изображений. При этом необходимо максимально использовать полученную информацию и получать конечные результаты с достаточной точностью и минимальными затратами при ее хранении, что связано со сжатием информации. Это важно при работе с геоинформацией в геодезии, землеустройстве, природообустройстве, мониторинге земель, при ведении точного сельского хозяйства [20, с.200]. При этом алгоритмы должны обладать простотой и точностью вычислений.
В настоящее время для моделирования объектов сжатия и фильтрации информации широко применяются преобразования Фурье [6, 9-11, 16-19, 23]. Однако даже в своем быстром варианте (быстрое преобразование Фурье - БПФ) оно сопряжено с большим количеством вычислений. В свою очередь в работах [11, 18] отмечается, что в отличие от преобразований Фурье возможны преобразования с другими базисами, восстанавливающие дискретные и непрерывные функции, но значительно сокращающие вычисления. Одним из таких базисов является вейвлетный. В настоящее время он постепенно находит практическое применение [7-15]. Следует при этом отметить, что работы [5-7, 21, 22] носят лишь ознакомительный характер, в основательной работе [21] о применении вейвлетов ничего не сказано, но детально описывается большинство фильтров подавления шумов, которые могут использоваться при моделировании процессов на основе вейвлетного базиса. В трудах [3, 17] приводятся основные теоретические положения по теории вейвлетов, базирующиеся на основе зарубежных исследований. Для создания технологий обработки геоинформации на основе этих работ необходимы дополнительные исследования. В работах [1, 11] выполнены значительные исследования по описанию вейвлетами гравитационного поля Земли. Показано, без детального описания технологии, что применение сжатия информации не всегда эффективно при работе с вейвлетами, хотя при использовании всей информации осуществляется быстрое и полное восстановление поля. В [24, 25] детально
Ж АС.Ярмоленко, О.В.Скобенко 001: 10.31897/РМ1.2018.6.612
:,и ■ Применение теории вейвлетов...
изложена теория лишь первого этапа вейвлет-преобразования, при этом она не доведена до универсального алгоритма. В [14, 15] выдвигается оригинальная идея комбинирования базиса Фурье и вейвлет-базиса при моделировании гравитационного поля Земли (ГПЗ), приводятся существующие вейвлеты, представлена лишь общая схема вейвлет-анализа-синтеза. В [4] в вейвлет-преобразовании используются зарубежные пакеты при неизвестном алгоритме.
Перечисленные труды имеют практическое и теоретическое значение при обработке геоинформации, но в них теория вейвлетов не доведена до инженерного уровня их применения. На основе отмеченного в данной статье с целью разработки детальной и доступной технологии применения вейвлетов в обработке геоинформации ставятся и решаются следующие задачи: 1) построение ортонормированных базисов вейвлет-преобразований и фильтрация; 2) создание вейвлет-фильтров; 3) построение фильтров последовательностью сверток; 4) применение КМА-анализа для построения ортонормированного базиса вейвлет-преобразования; 5) вейвлет-фильтрация и сжатие данных на основе составленных конкретных программ моделирования полей данных геонформации.
Базовые положения. В соответствии с [18] приведем следующие определения. Оператором циклического сдвига последовательности Z [18, с.125] на к-позиций вправо является оператор Як, осуществляющий создание новой последовательности Як2 по формуле
(Як2)(п) = 2(п - k), (1)
где п - номер элемента в создаваемой последовательности. Пусть имеется ряд (вектор значений)
2 = (1,0, 0,1)г (2)
при N = 4. Тогда при к = 1 (Я2)(0) = 2(0 - 1) = 2(-1) = - 1) = Z(4 - 1) = 2(3) = 1; (Я2)(1) = = 2(1 - 1) = 2(0) = 1; (Я2)(2) = 2(2 - 1) = 2(1) = 0; (Я2)(3) = 2(3 - 1) = 2(2) = 0. Таким образом получим новую последовательность-вектор
Я2 = (1, 1, 0, 0)г (3)
При к = 4 получаем снова исходный вектор
Я2 = (1,1,0,0)г. (4)
Согласно определению 3.7 и теоремы 3.8 [18, с.176] при некоторых заданных векторах и и V, принадлежащих тому же пространству элементов, что и вектор 2 , например (2), возможно построение ортонормированного базиса вида
В = {VМ1 и {ВДМ-1 = {V, Я?, Я^,..., RN-V, и, Я2и, Я4и,..., RN-2и}, (5)
где и, согласно [18], можно назвать отцовским вейвлетом, а V- материнским; символ и - логическое объединение множеств, М = N/2.
В данной статье составляющие базиса, построенные по вектору и, будем называть отцовскими, а по вектору V - материнскими.
Ортонормированность (5) возможна тогда лишь [13, с.250-258, лемма 7.1; 18, теорема 3.8], когда система матриц
( Т-7/„Ч т?/„ч
А(п)=т?
и(п) У(п) Щп + М) N+М)
(6)
для п = 0, 1, 2, ..., М- - унитарна.
Под унитарной понимается такая матрица [18, с.100], для которой
А-1 = А*, (7)
где А-1 - обратная к А матрица, а А* - сопряженная к А матрица, получаемая комплексно-сопряженной значением от всех элементов матрицы Аг, транспонированной к А. Комплексно-сопряженным к числу г = х + ¡у есть число г = х + ¡у. Для матриц вещественных чисел А = Аг,
и(п), V (п) _ соответственно преобразования Фурье заданных векторов и и V. Считается, что пара и и Vпорождает ортонормированный базис. В [18] ^мерное пространство элементов вектора 2 обозначено как ¡2(2м) пространство значений, в том числе и комплексных, на котором определена /2-норма [18, с.112]:
А = [£ I2 (к )|:
Множество индексов последовательностей обозначается как = {0,1, ..., N-1} [18, с.111]. Обобщением для двумерного случая является МД^мерное пространство ^(АдуАд^). Если принять, что {В0, В1, ..., ВД1_1} есть ортонормированный базис [18, с.134] ¡2(2Д]) и {С0, С\, ..., Сд-1} есть ортонормированный базис /2(2дт2), то ортонормированный базис пространства I2(2Д1-2Д2) будет Dm1,m2(nl, п2) = Вщ(п1)Ст2(п2), где 0 < Ш1 < N1 _ 1; 0 < щ < N1 _ 1; 0< Ш2 < N2 _ 1, 0 < П2 < N1 _ 1. Например, при
1
1 1 ^ V1 -
1 Г1 1 ^
V1 - Ъ
получим
Dm1,m2(nl, П2) --
Г1 1 1 1 ^
1 -11 -1
1 1 -1 -1
V1 -1 -1 1 У
Скалярным произведением (комплексным) векторов 7 и w называется выражение [18, с.91]
п
1-1
где есть число, комплексно-сопряженное Wj. Сверткой z*w называется вектор с компонентами
N-1
г * w(m) - ^ г(т - п^(п).
п-0
Для использования аппарата свертки в [18, с.172] вводится сопряженное отражение ю вектора ю:
ю(к - п) - ю(п - к).
тогда
N-1
N-1
г *ю(к) - £ 2 (п) *5(к - п)2 (п) * ю(п - к).
п-0
п-0
Построение ортонормированных базисов вейвлет-преобразований и фильтрация. В качестве примера построения ортонормированного базиса примем векторы Хаара [13, с.270; 18, с.190]:
и-| * , * ,0,...,0 | ;
л/2 л/2
(8)
V-|-^,0,...,0 | .
л/2 л/2
(9)
Поскольку в [18] установлено, что матрицы А(п) (6) для векторов и и V унитарны, то пользуясь правилом (5), построим ортонормированный базис для вейвлет-преобразования вектора 2, например (2). Здесь N = 4 и М = N/2 = 2.
2
А.С.Ярмоленко, О.В.Скобенко
Применение теории вейвлетов...
Тогда в соответствии с (1), (5) на основании (8), (9) можно составить ортонормированный вейвлет-базис:
В - (V, Я2V ,и, Я2и)
( 1 0 10 ^ -10 10 0 10 1
ч 0 -10 1 у
(10)
Тогда вейвлет-преобразование на данном этапе будет
В нашем примере
= В 2 .
( 1 0 10 у 1 ^
2 - -Г
л/2
-10 10 0 10 1
V0 -10 1А1 у
1
72
(1 ^ -1 1
V 1 у
(11)
(12)
а обратное преобразование, следуя (11), имеет вид
В примере это
2 - В2 - В*2 - ВТ2.
(1 -10 0 у 1 ^
2 - -Г
л/2
0 0 1 -1 110 0
V0 0 1 1 Л 1 у
-1 1
1
И
(1 ^ 0 0
V1 у
(13)
Теперь осуществим фильтрацию сигнала, т.е. его разложение по отцовским и материнским составляющим вейвлет-базиса. Согласно [13, с.26; 18, с.99] в пространстве с комплексным скалярным произведением <,> и ортонормированным базисом Я={ы\, и2, ..., ип} для любого V этого пространства справедливо
V-2
1 -1
V,и]и].
(14)
Здесь V соответствует вектору данных 2, а базисные векторы иу - всем столбцам матрицы (10). После подстановки их в (14) получим
где
2 - Q(г) + Р(2),
^2) - 2
(1 У 1 ^
(1001)
-1
V 0 у
-1
V 0 у
( 0 у 0 У ( ( 1 ^ ( 0 У ( 1 ^
+ (1001)
0 1
v- 1у
0 1
V- 1уу
-1
V 0 у
+(-1)
0 1
V- 1уу
-1 -1
V 1 у
(15)
р( 2) - ^
(1У1 ^
(1001)
(0 У 0 у
+ (1001)
0 1
V1 уууу
( (1 ^ 1 0
V V °у
+1"
( 0 У 0 1
V1 уу
(Г| 1 1
ч1У
Из примера видно, что Q(Z) и Р(2 составляющие сигнала высоких и низких частот соответственно. Следуя [18, с.180], составляющую высоких частот сигнала представим в виде формулы
*
1
0
0
0
А.С.Ярмоленко, О.В.Скобенко
Применение теории вейвлетов...
( N/2)—1
Q( 2) = £ (2, я2у)я2у = ((2, яу) )я0У
+
к=0
+
((2, яу) )Я2У +... +
( ,
\
2, V
\\ 2 /у
Ям
2
(16)
а составляющую низких частот в виде
(N/2)—1
Р(2) = £(2, Я^У = ((2,Щ)яи+
к=0
-((2, Я2и) +... +
Л
2, и
и.
2
(17)
Примеры фильтров на основе базисов Шеннона с вещественными и комплексными числами. Если воспользоваться составляющими векторов и и V (преобразований Фурье) в [18, с.179, 181], то при N = 4 с использованием обратных преобразований Фурье можно по аналогии с (14)-(17) определить высоко- и низкочастотные составляющие для данных примеров и в этих базисах (табл.1).
Таблица 1
Разложение вектора ^ по отцовскому и материнскому вейвлетам в различных базисах
Исходный Базисы
вектор 2 Хаара Шеннона 1 Шеннона 2
Р(2)
1 1 3 1 . 1
2 1 4 4 2
0 1 2 11 — - — г 4 4 0
0 1 11 — + — г 44 1
2 2
1 1 2 3 1 — + — г 44 1
6(2)
1 1 11 — + — г 44 1
2 2
0 1 2 11 - — + — г 44 0
0 1 1 1 1
2 — — г 44 2
1 1 2 11 — — — г 44 0
Из сравнения базисов, приведенных в табл.1, а также базисов Мейера, Баттла - Лемарье, Добеши [13, с.207] преимущество следует отдавать базису Хаара по следующим причинам:
• базис Хаара прост в вычислении;
• понятна на его основании и фильтрация сигнала; так, низкочастотная часть на первом этапе равна нулевому коэффициенту преобразования Фурье, высокочастотная - соответствует отклонениям сигнала от его середины;
• в базисах, отличных от Хаара, задаются дополнительные требования к числу N. Например, в базисах Шеннона оно должно быть кратным 4, а в базисе Хаара кратности 2р при р = числом N определяется и вид базиса Добеши [18].
Ж АС.Ярмоленко, О.В.Скобенко DOI: 10.31897/PMI.2018.6.612
:,и ■ Применение теории вейвлетов...
Построение фильтров последовательностью сверток. В вейвлет-разложениях (16), (17) замечено, что коэффициенты при R2kV, R2kU являются соответствующими скалярными произведениями, которые в виде сверток можно записать так [13, 18, 19]:
Z, R2kv) = ZV~(2k), Z, RiU) = ZU~(2k).
Настоящие выражения являются теоретической основой быстрого вейвлет-преобразования. Тогда фильтр на основе такой свертки строится следующим образом [13, с.279; 18, с.187]: 1) образуются векторы U, V (описаны в конце базовых положений); 2) осуществляются свертки Z*U, Z*V; 3) вводится оператор децимации - удаления составляющих с нечетными номерами D(ZU), D(ZV); 4) реализуется оператор разрежающей выборки - удваивания размера вектора вставкой нуля между двумя смежными значениями U(D(Z*U)), U(D(Z*V)); 5) осуществляется фильтрация -строится вектор низкочастотной составляющей в виде свертки P(Z) = U*U(D(Z*U)) и высокочастотной Q(Z) = V*U(D(Z*V)); 6) восстанавливается сигнал Z = P(Z) + Q(Z).
На этом заканчивается первый этап разложения сигнала на высоко- и низкочастотные составляющие. Число всех этапов определяется по формуле р = log2N. Каждый последующий этап состоит из фазы анализа и фазы синтеза. В фазе анализа на этапе n осуществляется: 1) ввод вектора ZTn_x *Un_j и его децимация Zn = D(ZTn_x *Un_x); 2) ввод векторов Un, Vn нормированного базиса размерности N\ = N/2"-1; 3) осуществляются децимация и разреживание сверток ZTn *Un, ZTn *Vn; UDUn = UD(ZTnUn), UDVn = UD(ZTnVn); 4) вычисляются высокочастотная и низкочастотная составляющие вектора Z„: Q(Z„) = Vn*U(D(Zn* Vn)); P(Z„) = U*U(D(Zn*U„)).
В фазе синтеза осуществляется: 1) разреживание векторов Q(Zn), P(Zn) (при этом разреживание выполняется числом n - 1-раз до достижения размерности исходного сигнала): U(Q(Zn)), U(P(Zn)); 2) операциями свертки получают высоко- и низкочастотную составляющие сигнала на этапе n: Qn(Z) = U(Q(Z„))*Ub P„(Z) = U(P(Zn))*Ul.
Окончательный результат анализа синтеза:
Z = P( Z) + 2 Qp_x( Z).
i=0
Настоящий алгоритм реализован в специально составленной авторами программе Sub Вейв-лет Анализ Синтез ().
Сигнал Z можно представить в виде разложения
Z = c0b0 +C1b1 + ... + CN _1bN (18)
по базису B, который представляет собой совокупность ортонормированных базисов [18, с.176, с.185-186] B = {bi }N=01 .Очевидно, что в этом случае ci = (ZbT). Задача заключается в том, чтобы определить ортонормированные векторы в разложении (18). В теории вейвлетов [18, с.209-225, определение 3.28] вместо (18) принимается следующая запись:
( N/2 p_1 _1) ( N/2 p_( p _1) _1)
+ 2 C2,k p,k + 2 С3
k=0 k=0 k=0 k=0
Z = 2 c1,k Ф_ p,k + 2 c2,k p,k + 2 c3,k V_( p_1),k +... + 2 CN/2,k V_( p_( p_1),k . (19)
Если принять р = то (19) перепишется в виде
Z = С1,0ф_р,0 + С2,0ф-р,0 + С2,оф-р,1 + • • • + + С#/2,0ф-1,0 + сМ2,1ф-1,1 + — + СМ2,М2-1ф-1,М2-1. (20)
Здесь все векторы базиса В ф_д записаны слева направо по степени детализации вектора 2. В [18] они записаны наоборот справа налево. Построение базисных векторов ф_д осуществляется в следующем порядке [18, с.225]: 1) применяется последовательность вейвлет-фильтров и1, У\, и2, У2; ...; ир, Ур; при этом и1, У1 е 12( например, исходя из (8), (9) при 1 = 2 будет
U2 1 l V2 j
_ 1
V2
Г1 1 0 0^ 1 _ 10 0
; 2) каждый из базисов строится по формуле [18, с.225]
Ж АС.Ярмоленко, О.В.Скобенко 001: 10.31897/РМ1.2018.6.612
:,и ■ Применение теории вейвлетов...
V— ;,к = Я2 >к/1 ,
где/ = gl-lUl-1(Vl), а gl = g/-1U/-1(Ul) при начальных значениях, равных/1 = VI, gl = Ц\. В приведенных выражениях оператор Ul-1(Vl) означает /-1-кратное разреживание вектора Vl или и1. Элементы сверток/ (п), g/ (п) для п = 0, 1, ..., N-1 вычисляются по формулам
N-1 N-1
fx(«) = Z(gi-i(m)Ul-1(Vl)(n - m)); gl(n) = £(gl-l(m)Ul-1(Ul)(n - m)) .
m=0 m=0
КМА-анализ для построения ортонормированного базиса вейвлет-преобразования. В фундаментальных работах [13, 18, 19] вейвлет-разложение предполагает наличие масштабирующей и уточняющей функций. В [18, определение 5.30] масштабирующая функция еще называется отцовским вейвлетом. В [13, с.103; 19, с.44] термин «отцовский вейвлет» в определении масштабирующей функции не применяется. В тех же работах [18, с.198; 19, с.39-42, гл.5] уточняющая функция (материнский вейвлет) просто названа вейвлетом.
В системе Хаара в соответствии с [18, (5.49)] отцовский вейвлет записан в виде
Г1, 0 < х < 1;
Ф( х) = Г (21)
[0, иначе.
Материнский вейвлет по формуле [18, с.380]
у( х) = ф(2 х — 1) — ф(2 х), (22)
и, что легко показать, имеет вид
у( х) =
-1, 0 < х <-;
2
1, 1 < х < 1; (23)
2
0, иначе.
Отцовский и материнский вейвлеты позволяют построить ортонормированный базис вейв-лет-преобразования для представления дискретно-заданной геоинформации с целью ее обработки. Для этого в вейвлет-теории разработан так называемый кратно-масштабный анализ (КМА-анализ). На его основе создается удобный в использовании алгоритм построения ортонормированного базиса вейвлет-преобразования. КМА-анализ базируется на функциях
Фм = 2-2 ф(2-jx - k); (24)
у;,k = 2-* у(2-'х - k). (25)
Эти формулы приведены в [5, с.113, с.128; 13, с.241; 18, с.368; 19, с.193]. При этом в приведенных работах показателям степени в (24), (25) приписывается как положительный, так и отрицательный знак, как в нашем случае. В случае отрицательной степени график функции растягивается по оси х, а при положительной - сжимается. В работе нас интересует растяжение по оси х, поэтому принята запись степени с отрицательным знаком. При такой записи уточняются значения функций в зависимости от числа ортонормированных векторов базиса вейвлет-разложения, т.е. осуществляется увеличение деталей анализируемой информации или увеличение разрешения. Поэтому в зарубежной литературе [5, с.113] КМА-анализ справедливо называется многоразрешающим анализом (multiresolution analysis). Построение вейвлет-базисов в системе Хаара будем вести на основе (21)-(25) в следующем порядке:
1) построение базисного вектора нулевого приближения ф;
2) построение последующих уточняющих базисных векторов.
1. Построение базисного вектора нулевого приближения. В основу построения всех базисных векторов положим вейвлет-базис вида [18, (3.67), (3.84)], применяемый нами в разложении (20). Хотя в теории вейвлетов [3, 5, 13, 17-19] допускается возможность построения нескольких
А.С.Ярмоленко, ОЯ.Скобенко
Применение теории вейвлетов...
базисных векторов нулевого приближения, в данной работе мы ограничимся лишь одним в соответствии с данным разложением. При этом такое разложение наиболее часто применимо в различных рядах, в том числе и при разложении в ряд Фурье.
В соответствии с разложением (20), принятым здесь за основу, базисный вектор ф-р, с единствен, так как С «пробегает» значения от 0 до ((N2) - 1) [18, (3.67)]. (В КМА-анализе обозначениям ф_р, с, у-р, с разложения (20) соответствуют ф-р, с, у-р, с.) Поскольку N = 2Р, то С = 0. Следовательно, в (24) > = р, С = 0. В соответствии с (21) 0 < 2-рх - С < 1, или С < х < (1 + С)/2-р. При С = 0 будет 0 < х < 1/2-р и составляющие этого вектора определятся по формуле
Фх)=2—р я,0 <х <2 р;
.0,
Таким образом, на основе (26) получим вектор
иначе.
(26)
—р_
фр,0 = 2 2(1,1,...,1)
гр,0 - (27)
с количеством одинаковых членов 2р.
2. Построение последующих уточняющих базисных векторов. Последующие уточняющие базисные векторы вычисляются по (25) с учетом (23). Порядок множества непересекающихся базис. векторов [18, (3.67)] и{у—> С 2 ^ определяется значением>. Это значение изменяется отр - 1
ных
до 0 с шагом - 1. Тогда число базисных уточняющих векторов порядка > составит величину N/21. Значение С будет изменяться от 0 до N/21 — 1. Например, для > = р будем иметь следующий уточняющий вектор (единственный): у-р, 0. Для> = / уточняющими векторами порядка / будут: у^, 0; у^, 1 ; у—/ ^21 )-1. Так, при N = 8 р = 3 уточняющими векторами порядка / = 1 при верхнем пределе С, равном
N/21 - 1 = 2р/2-1= 3, будут: у-1,0; у-1,1; у-1,2; Уу. Каждый из этих векторов при определенном С также определяется по (25) с учетом материнского вейвлета (23). Тогда можно записать
У/с = 2 2
— 1, 0 < 2—/х — С < 1 2
1, 1 < 2—/х — С < 1; 2
0, иначе
или
У/с = 2 2
— 1, С21 < х <|1 + С |21;
1, ^1 + С|21 < х < (1 + С)21;
2
0, иначе.
(28)
В качестве примера возьмем р = 3, / = 1, С = 2. Тогда
— 1, 4 < х < 5;
1
= 21
1, 5 < х < 6; 0, иначе.
В табл.2 дан пример вейвлет-разложения как методом сверток, так и КМА-методом.
Таблица 2
Разложение вектора ^ по составляющим ф и у
2
Исходный Составляющие
вектор с1,0-ф-3,0 с1,0-У-3,0 с1,0-У-3,0 с3,ГУ-2,1 04ЖУ-1,0 с4,ГУ-1,1 с4,2-У-1,2 с4,3-У-1,3
4 7,25 -3,25 -1 0 1 0 0 0
2 7,25 -3,25 -1 0 -1 0 0 0
3 7,25 -3,25 1 0 0 -2 0 0
7 7,25 -3,25 1 0 0 2 0 0
10 7,25 3,25 0 -1,5 0 0 1 0
8 7,25 3,25 0 -1,5 0 0 -1 0
10 7,25 3,25 0 1,5 0 0 0 -2
14 7,25 3,25 0 1,5 0 0 0 2
Ж АС.Ярмоленко, О.В.Скобенко DOI: 10.31897/PMI.2018.6.612
:,и ■ Применение теории вейвлетов...
Каждая составляющая вычислялась по формуле:
для низкочастотной части
N/2 p - j
Pu = 1/2p-j ZP,j-1 i=1
и высокочастотной
Qij = Py-1 - Pij
Здесь j - номер этапа; i - номер составляющей сигнала в группе (группа состоит соответственно из 2, 4,... 2р-элементов в зависимости от порядкового номера этапа): p - число этапов. Окончатель-
p-i
ный результат: Z = P1 + ZJQp-i , где P1 - составляющая нулевой частоты (вектор у-3,0 в примере
i=0
табл.1); Qp-,i - уточняющие составляющие. В приведенных преобразованиях на каждом этапе при нецелом p возможен остаток элементов, число которых меньше числа 2ip при ip = 1, 2, ..., 2p. По этому остатку находится также среднее значение TZ, которое записывается в старшей строке pip на место этого остатка элементов. В последней строке, следующей сразу после строки с номером целой части p, обозначаемой через pf, находится среднее по всем элементам предыдущей строки Ppf. Это среднее является одним и тем же для строки Ppf + 1. Подстрока qpf + 1 вычисляется в общем порядке:
qpf+1 = Ppf- Ppf + 1
Исследование эффективности сжатия геоинформации и фильтрации шумов с помощью вейвлетов. В основу исследований положены специально составленные авторами программы на языках VISUAL BASIC Excel (VBE) и IDL системы ENVI по сжатию и фильтрации геоинформации. В качестве первого объекта исследований принята модель рельефа, приведенная в [23]. В VBE-программе Sub МакросВЕЙВЛсИстОш() истинные высоты точек представлены массивом сс1(), а высоты, отягощенные случайными ошибками, - массивом cm(). Так же, как и в [23], связь настоящих массивов определяется формулой cm(i) = cc1(i) + delta, где i изменяется от 0 до N - 1, а delta = Randbetwen(-t, +t)Std - функция языка VBE-генерирования случайного числа в интервале значений квантиля от -t до +t; Std - задаваемый стандарт случайных ошибок (шумов). Именно высоты массива cm() подвержены вейвлет-разложению по частотам. При вейвлет-сжатии и соответственно фильтрации оставались самые большие по амплитуде члены разложения на всех частотах. Фильтрационный вектор представлен массивом Filt = Array(1; 0,8; 0,6; 0,5; 0,25; 0,15; 0,1; 0). Например, при значении фильтра Filter = Filt(1) = 0,8 на всех частотах разложения оставались значения более 0,8, остальные значения обнулялись. На выходе формировался отфильтрованный массив высот Tw(). Среднее квадратическое отклонение (СКО) отфильтрованных высот определялось по формуле
CKO =
N
Z (Tw(i) - cc(i))2
i =1_
N
Таблица 3
Средние квадратические отклонения (СКО) по каждому фильтру (в метрах) в зависимости от стандартов (Std) распределения случайных ошибок высот
Filter Оставшаяся информация, % Значения Std
0,05 0,1 0,3 0,4 0,5
1 25 0,66 0,88 0,71 0,93 0,94
0,8 33 0,55 0,35 0,63 0,86 0,88
0,6 40 0,23 0,26 0,63 0,86 0,94
0,5 48 0,23 0,26 0,64 0,86 0,94
0,25 51 0,19 0,24 0,63 0,86 0,99
0,15 66 0,15 0,20 0,64 0,86 0,99
0,1 70 0,11 0,17 0,63 0,86 0,99
0 100 0,09 0,18 0,63 0,86 0,99
В табл.3 приведены значения СКО в зависимости от значения фильтра - Filter и точности отмоделированных высот - величины Std. Оставшаяся информация (в процентах) вычислена как процентное отношение числа наибольших оставшихся коэффициентов к их общему числу до сжатия, равному 15.
Влияние систематической части вейвлет-сжатия на точность конечного результата после обнуления деталей приведено ниже.
СКО по каждому фильтру при высотах, не отягощенных случайными ошибками:
Filter 0,55 0,34 0,25 0,23 0,15 0,10 0,06 0 СКО 1 0,8 0,6 0,5 0,25 0,15 0,1 0
Из экспериментальных исследований следует:
• для сохранения максимально высокой точности исходных высот при их вейвлет-разложении сжатие информации недопустимо; даже при отсутствии сжатия (Filter = 0) и наличии случайных ошибок среднее квадратическое отклонение результирующих высот больше их стандартного отклонения на входе;
• сжатие информации можно допускать, но при этом следует учитывать порог понижения точности изображения рельефа; например, при трехкратном сжатии (Filter = 0,8) и Std = 0,1 среднее квадратическое отклонение высот на выходе СКО = 0,35 больше Std в 3,5 раза, в остальных случаях оно больше в 2 раза. На рис.1 показан рельеф, построенный по истинным высотам, а на рис.2 - рельеф, полученный при трехкратном вейвлет-сжатии (Filter = 0,8) и стандарте ошибок высот 0,3 м (Std = 0,3 м);
• закономерности влияния вейвлет-сжатия на точность получаемой на выходе информации такие же, как и при фильтрации информации в рядах Фурье [23], но здесь объем вычислений ничтожно мал по сравнению с преобразованиями Фурье.
Методика расчета значения фильтра для сжатия информации о высотах рельефа может быть принята следующей:
1. По известной методике определяется оценка стандарта высот рельефа (Std) в виде средней квадратической ошибки съемки рельефа.
2. Для данного объекта определяются средние квадратические отклонения по каждому фильтру.
3. Если значение средней ошибки рельефа, полученное как 0,8Std [2], меньше трети высоты сечения рельефа, то по полученной экспериментально можно выбрать порог (величину СКО) понижения точности сжатой информации о рельефе при определенном значении фильтра (Filter). Если средняя ошибка рельефа больше трети высоты сечения рельефа и приближается к величине СКО, получаемой при 2530 %-ном сжатии информации (т.е. трех- и четырехкратном сжатии), то в таком случае значение порога (величины СКО) понижения точности сжатой информации о рельефе не суть важно.
Одновременно авторами составлена программа pro oroi_data_ corr24bitWAVE на алгоритмическом языке IDL системы ENVI вейвлет-
Рис. 1.Рельеф в горизонталях, построенный по высотам, принятым в качестве исходных [23] (сечение рельефа 0,25 м)
Рис.2.Отфильтрованный рельеф при трехкратном вейвлет-сжатии информации (Filter = 0,8, Std = 0,3 м)
Рис.3. Изображения - исходное (восстановленное по всем вейвлет-элементам) (а) и сжатое в 2,5 раза (б) с отфильтрованными значениями (менее 4 единиц) составляющих вейвлет-разложения по модулю
фильтрации (сжатия) изображений. На рис.3 приведены изображения - исходное (и восстановленное по всем вейвлет-элементам) и сжатое в 2,5 раза с отфильтрованными значениями (менее 4 единиц) составляющих вейвлет-разложения по модулю.
Кружками отмечены несовпадения с оригиналом. При большем сжатии различных несовпадений кружки больше и представляют собой уже шумы, при которых изображение непригодно к использованию. Таким образом, при сжатии изображения рекомендуется также подобрать вначале фильтр в соответствии с пунктами 1-3 в случае рельефа, или при котором отсутствуют видимые несовпадения. И лишь при нем осуществлять вейвлет-сжатие.
Выводы
1. Вейвлет-фильтрация и сжатие приводят к потере части информации. Однако такая фильтрация и сжатие весьма эффективны при известном допустимом пороге (величине СКО) понижения точности сжатой информации. При этом объем вычислений ничтожно мал по сравнению с фильтрацией и сжатием в рядах Фурье и косинусно-синусных преобразованиях JPEG-сжатия.
2. При фильтрации и сжатии изображений рекомендуется подобрать вначале фильтр в соответствии с пунктами 1-3 в случае рельефа. При этом необходимо задавать СКО значения пиксела, при котором допустимо сжатие и отсутствуют видимые несовпадения.
ЛИТЕРАТУРА
1. БагровА.А. Разложение сферических функций по вейвлетам Хаара / А.А.Багров, А.С.Багрова // Известия вузов. Геодезия и аэрофотосъемка. 2008. № 4. С. 6-8.
2. Большаков В.Д. Теория математической обработки геодезических измерений / В.Д.Большаков, П.А.Гайдаев. М.: Недра, 1977. 367 с.
3. Воробьев В.И. Теория и практика вейвлет-преобразования / В.И.Воробьев, В.Г.Грибунин. СПб: ВУС. 1999. 204 с.
4. Гонжа Е.А. О вейвлет-фильтрации цифрового изображения земной поверхности // Известия вузов. Геодезия и аэрофотосъемка. 2017. № 3. С. 105-110.
5. ДьяконовВ.П. Вейвлеты. От теории к практике. М.: Солон-Р, 2002. 448 с.
6. Журкин И.Г. Автоматизированная обработка данных дистанционного зондирования / И.Г.Журкин, Н.К.Шавенько. М.: Диона, 2013. 456 с.
7. Красильников Н.П. Цифровая обработка 2D- и 3D-изображений. СПб: БХВ-Петербург, 2011. 608 с.
8. Лапшин А.Ю. Разработка и исследование методики вычисления гравиметрической высоты квазигеоида и составляющих уклонения отвеса на основе вейвлет-преобразования: Автореф. дис. ... канд. техн. наук / Московский университет геодезии и картографии. М., 2011. 207 с.
9. Мазурова Е.М. Алгоритмы быстрого преобразования Фурье // Известия вузов. Геодезия и аэрофотосъемка. 2004. № 3. С. 18-35.
10. Мазурова Е.М. Двумерное и матричное представление быстрого преобразования Фурье // Известия вузов. Геодезия и аэрофотосъемка. 2004. № 4. С. 3-12.
11. Мазурова Е.М. К вопросу о вычислении аномалии высоты на основе вейвлет-преобразования и быстрого преобразования Фурье в плоской аппроксимации / Е.М.Мазурова, А.С.Багрова // Известия вузов. Геодезия и аэрофотосъемка.
2008. № 4. С. 6-8.
12. Малинников В.А. Анализ методов формирования мультифрактальной меры, основанных на вейвлет-обработке экспериментальных данных / В.А.Малинников, Д.В.Учаев // Известия вузов. Геодезия и аэрофотосъемка. 2007. № 6. С. 57-61.
13.Мала С. Вейвлеты в обработке сигналов: Пер с англ. М.: Мир, 2005. 671 с.
14. Нейман Ю.М. Об адаптации глобальной модели геопотенциала к региональным особенностям. Ч.1 / Ю.М.Нейман, Л.С.Сугаипова // Известия вузов. Геодезия и аэрофотосъемка. 2014. № 3. С. 3-12.
15. Нейман Ю.М. Основы разномасштабной аппроксимации геопотенциала / Ю.М.Нейман, Л.С.Сугаипова. М.: Изд-во МИИГАиК, 2016. 218 с.
16. ПискуновН.С. Дифференциальные и интегральные исчисления: В 2-х т. М.: Наука. 1978. Т.2. 575 с.
17. Уэлстид С. Фракталы и вейвлеты для сжатия изображений в действии. М.: Изд-во Триумф, 2003. 320 с.
18. Фрейзер М. Введение в вейвлеты в свете линейной алгебры: Пер. с англ. М.: Бином. Лаборатория знаний. 2008. 487 с.
19. Чуи К Введение в вейвлеты: Пер. с англ. М.: Мир, 2001. 412 с.
20. Шпарр Д. Точное сельское хозяйство (Precision agriculture) / Д.Шпаар, А.В.Захаренко, В.П.Якушев /Федеральное министерство продовольствия, сельского хозяйства и защиты прав потребителя Германии. Санкт- Петербург - Пушкин.
2009. 398 с.
21. Шовенгердт Р.А. Дистанционное зондирование, модели и методы обработки изображений. М.: Техносфера, 2010.
560 с.
22. ЯковлевА.Н. Введение в вейвлет-преобразование. Новосибирск: Изд-во НГТУ, 2003. 104 с.
23. ЯрмоленкоА.С. Фильтрация геоинформации в рядах Фурье / А.С.Ярмоленко, О.В.Скобенко // Известия вузов. Геодезия и аэрофотосъемка. 2016. № 1. С. 107-113.
24. ЯрмоленкоА.С. Использование вейвлетов в аналитическом представлении дискретных функций графической информации // Известия вузов. Геодезия и аэрофотосъемка. 2008. № 3. С. 20-30.
25. Ярмоленко А.С. Вейвлет-преобразования в кодировании графической информации // Известия вузов. Геодезия и аэрофотосъемка. 2010. № 4. С. 18-25
Авторы: АС.Ярмоленко, д-р техн. наук, профессор, yarmolenko_alex@mail.ru (Новгородский государственный университет им. Ярослава Мудрого, Великий Новгород, Россия), О.В.Скобенко, ассистент, olga-skobenko@mail.ru (Белорусская государственная сельскохозяйственная академия, г. Горки, Могилевская область, Республика Беларусь). Статья поступила в редакцию 19.04.2018. Статья принята к публикации 21.09.2018.
