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

ISSN 0868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2008, том 18, № 2, c. 61-65

ИНФОРМАЦИОННО-ИЗМЕРИТЕЛЬНЫЕ СИСТЕМЫ

УДК 519.218.82

© Т. Н. Князева, Л. В. Новиков, Н. И. Орешко

УДАЛЕНИЕ НЕСТАЦИОНАРНОГО ШУМА ИЗ ЭКСПЕРИМЕНТАЛЬНЫХ ДАННЫХ

В статье приводится двухступенчатый метод очистки данных от шума, дисперсия которого изменяется в случайные моменты времени. На первом этапе на основе предложенного метода кластеризации вейвлет-коэффициентов, вычисленных по максимально накладывающемуся дискретному вейвлет-преобразованию, определяются участки шума с кусочно-постоянной дисперсией. На втором этапе производится расчет модифицированных порогов для каждого участка с последующим пороговым удалением шума. Предлагаемый подход по сравнению с традиционными методами очистки от шума, в которых порог рассчитывается по всей выборке данных, позволяет значительно эффективнее удалять локальный шум.

ВВЕДЕНИЕ

Результаты измерений многих физических процессов представляют собой аддитивную смесь квазидетерминированного сигнала и нестационарного шума, дисперсия которого изменяется в случайные моменты времени (гетероскедастический шум). В этом случае наиболее эффективным средством борьбы с шумами является вейвлет-технология, изначально предназначенная для исследования нестационарных процессов [1, 2].

В частности, вейвлетные подходы широко используются в задачах непараметрической регрессии для оценивания функции / на основе наблюдений у (^) в дискретные моменты времени ti:

У (t, ) = f (t,) + £, i = 0,1,2,...,N -1,

(1)

где /(^) — истинное значение сигнала; N — объем выборки; ^ — шумовая составляющая, распределенная по нормальному закону с параметрами £ ~ N(0,1); ai — среднеквадратическое отклонение (СКО) шума. (В дальнейшем для простоты полагаем ti =;). При этом стандартные ге-тероскедастические модели описания дисперсии погрешности измерений (шума) представляют собой достаточно гладкую и, кроме того, известную кривую [3].

При использовании вейвлетов пороговое значение в алгоритмах удаления шума формируется на основании оценки СКО, вычисленной по всей выборке [1]. Однако это приводит к заметному снижению качества очистки на участках, где шум значительно отличается по СКО от глобально вычисленного значения.

В настоящей работе решается задача удаления шума, дисперсия которого изменяется скачком

в случайные моменты времени, т. е. c — кусочно-постоянная функция: c = const при значениях i из некоторого замкнутого интервала, образующего k-й кластер, k = 1,2,...,M, M — число кластеров, M << N . Для этого в работе предлагается двухступенчатый способ очистки от шума, где на первом этапе на основе кластеризации (разбиения отсчетов на группы с близкими свойствами) вейв-лет-коэффициенты (ВК) каждого уровня декомпозиции определяют группы данных с постоянной дисперсией, а затем на втором этапе производится расчет порогов для каждого кластера и пороговое удаление шума. При этом используется предложенное в [4] максимально накладывающееся дискретное вейвлет-преобразование (МНДВП, "maximal overlap discrete wavelet transform").

МАКСИМАЛЬНО НАКЛАДЫВАЮЩЕЕСЯ ДИСКРЕТНОЕ ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЕ

МНДВП является избыточным неортогональным преобразованием. Несмотря на это оно обеспечивает более качественную очистку от шума и обладает рядом преимуществ, в частности инвариантностью относительно сдвига и способностью работать с выборками произвольного объема.

На первом уровне вейвлет-преобразования ВК W1(i) и аппроксимирующие коэффициенты V1(i) рассчитываются по формулам:

W (i) = £g,y (( -1)modN),

l=0

V (i ) = £ hl, y ((i -1 )mod N),

где

l=0

62

Т. Н. КНЯЗЕВА, Л. В. НОВИКОВ, Н. И. ОРЕШКО

Г/ -1, если г -1 > 0, (г -1)modN = <!

[г - 1 + N, в остальных случаях;

, ¡г1 — коэффициенты высокочастотного и низкочастотного фильтров соответственно. Они могут быть получены из фильтров Добеши &1, к1, например, путем нормирования на л/2, т. е.

& = gl /Л

а=у л, таким образом что

ь-1 Ь-1

£ &1 = 0 и £ а = 1, где Ь — длина фильтра. 1=0 1=0

ВК и аппроксимирующие коэффициенты ] -го уровня вычисляются по формулам

^ (г) = £ &У 1 ((г - 2'-11 )N),

1=0 Ь-1

(3)

V (г) = £а1У]-1 ((г - 2'-11 )modN);

1=0

при этом У0 (г )= X (г). Эти два соотношения составляют пирамидальный алгоритм МНДВП. Обратное МНДВП вычисляется с помощью обратного пирамидального алгоритма по следующей формуле:

V-1 (г ) = £ ((г + 2-11 ) N

+

1=0 Ь-1

(4)

+ £011У] (( + 2-11 )modN).

ПРОВЕРКА ГИПОТЕЗЫ О ГЕТЕРОСКЕДАСТИЧНОСТИ ШУМА

Проверка данных на наличие в них шума с кусочно-постоянной дисперсией осуществляется по ВК первого уровня декомпозиции Щ(г) (определенным по формуле (2)), соответствующим высокочастотной шумовой составляющей. Для этого ВК Ж1 (г) разбиваются на р одинаковых интервалов с помощью окна выбранной длины, соблюдая компромисс между точностью оценивания дисперсии и возможностью отслеживания локальных изменений дисперсии. По каждому интервалу определяется выборочная дисперсия Д*, к = 1,...,р .

Затем проверяется гипотеза Н0 об однородности полученной последовательности из выборочных дисперсий Д*. Известны способы проверки этой гипотезы: с использованием критериев Барт-летта [5], Кохрана или на основе вычисления нормализованной суммы квадратов дисперсий ВК первого уровня декомпозиции [6]. В настоящей

работе используется критерий Бартлетта, суть которого заключается в следующем. Для р выборок Х1,Х2,...,Хр вычисляются выборочные дисперсии

Д ,Д,...,Др. При заданном уровне значимости д принимают гипотезу Н0 (о равенстве дисперсий Д = Д =... = Др) в случае выполнения следующего неравенства: А/В <Х% (р -1), где

Х2-д (р -1) — X2 -распределение Пирсона с р -1

степенями свободы. Величины А и В вычисляются по формулам, приведенным в [5].

УДАЛЕНИЕ ШУМА

Традиционная схема удаления шума включает этапы ДВП, пороговую обработку ВК с вычислением порогов по глобальной дисперсии и обратного ДВП. Предлагается модифицированный алгоритм удаления шума, который применяется в случае выполнения гипотезы о его гетероскедастич-ности, отличительной чертой которого является введение процедуры кластеризации ВК каждого уровня декомпозиции и расчет порогов по кластерам.

Прямое ДВП. Сначала выполняется прямое ДВП по формулам (3). Для снижения граничных эффектов на каждом уровне декомпозиции необходимо увеличивать выборку на длину фильтра слева и справа. Для этого можно использовать, например, полиномиальное или симметрическое продолжение в зависимости от вида сигнала [7].

Кластеризация ВК. На разных уровнях декомпозиции вследствие особенностей МНДВП границы кластеров сдвигаются, поэтому для корректной очистки сигнала от шума необходимо проводить кластеризацию на каждом уровне декомпозиции, вычисляя границы кластеров с одинаковой дисперсией. В работе предлагается использовать алгоритм, основанный на сравнении дисперсий с применением F-распределения Фишера [5] и дальнейшим уточнении границ. Преимущества этого алгоритма состоят в том, что, во-первых, он не требует априорного знания числа кластеров, которое автоматически определяется в процессе кластеризации, во-вторых, алгоритм обеспечивает высокую точность разделения выборки на кластеры за счет разработанной процедуры уточнения границ.

Сначала проводится предварительная кластеризация. На каждом уровне j декомпозиции последовательность ВК разбивается на р интервалов заданной длины (например, по 10 коэффициентов). Далее на каждом интервале определяется вы-

борочная дисперсия 5*к (к = 1,...,р), где к —

номер интервала. Определяются все дисперсии, попадающие в доверительный интервал при заданном уровне значимости q :

Рч,2 (г,г Щ < 5),к < ^2 (г,г) 5* ,

где ^ — распределение Фишера с числом степеней свободы г, 5* — медиана из р дисперсий.

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

В результате выполнения алгоритма на разных уровнях декомпозиции может оказаться неодинаковое число pj кластеров (к = 1,2,...,pj), каждый

из которых объединяет ВК с индексами, образующими неперекрывающиеся подмножества К, к .

В этих обозначениях ВК, входящие в к -й кластер , -го уровня декомпозиции имеют индексы

1 <с К ,

У/2

(5)

где — среднее значение ВК, входящих в кластер к ; zk — число ВК в к -м кластере.

Пороговые значения рассчитываются в зависи-

мости от выбранного типа порога. Например, если порог многоуровневый универсальный, то пороговые значения определяются для каждого уровня вейвлет-преобразования и для каждого кластера по формуле

(6)

На следующем этапе, с учетом рассчитанных пороговых значений и выбранного правила усечения (сравнения с порогом) производится корректировка ВК. Например, если правило усечения жесткое [1], то вейвлет-коэффициенты, прошедшие пороговую обработку, Ж^ () вычисляются по формуле

Ж (() =

Ж, (), если

0,

если

Ж ()|> т; Ж (()< т,

(7)

Это означает, в частно-

/ е К,,к, т. е. Ж, (

сти, что каждый кластер может включать в себя несколько участков коэффициентов Ж, () с индексами, например, [1-20], [45-80] и т. п.

Так как начальная кластеризация проводилась по выборочным дисперсиям, рассчитанным на определенных интервалах, то на границах этих интервалов некоторые ВК могли попасть не в свой кластер, поэтому необходимо уточнение границ кластеров. Уточнение границ участков кластеров осуществляется путем минимизации или максимизации локальных дисперсий двух соседних кластеров при перемещении границы, которая их разделяет.

Пороговая обработка. По полученным кластерам на каждом уровне декомпозиции пересчиты-ваются дисперсии и вычисляются пороги. Формула расчета СКО для к -го кластера на , -м уровне декомпозиции имеет вид

где Т = ^г, к при ; е К, к.

На заключительном этапе выполняется обратное МНДВП, обеспечивающее восстановление исходного сигнала, используя измененные ВК всех уровней разложения и аппроксимирующие коэффициенты последнего уровня по формуле (4).

МОДЕЛИРОВАНИЕ

Эффективность применения рассмотренных выше алгоритмов проверена на модели (1), состоящей из шести пиков гауссовой формы / () (см. рис. 1, а), искаженных медленно линейно нарастающим шумом £ и шумом со скачками дисперсии с, в случайные моменты времени и случайной продолжительности (см. рис. 1, б). Проверка шума на гетероскедастичность проводилась по ВК первого уровня декомпозиции, вычисленных по формуле (2), разбитым на р = 103

интервала (рис. 2, а). В качестве фильтров Нп, gn приняты коэффициенты вейвлетов Добеши 4-го порядка. Использовался критерий Бартлетта при уровне значимости q = 0.95. Отношение величин Л/Б оказалось равным 1702.5, что больше квантиля %2 -распределения Пирсона, взятого из таблиц, — 126.5741 [5]. Следовательно, гипотеза Н0 об однородности дисперсии, как и следовало ожидать, в данном случае не выполняется. Этот факт иллюстрируется на рис. 2, б, где показаны величины выборочных дисперсий в каждом интервале. Затем вычислялись ВК по формулам (3) до уровня , = 5; на каждом уровне проводилась кластеризация и уточнение границ кластеров в соответствии

64

Т. Н. КНЯЗЕВА, Л. В. НОВИКОВ, Н. И. ОРЕШКО

Щ (7)

Рис. 1. Удаление нестационарного шума. Кривые: (а) — полезный сигнал; (б) — наблюдаемый сигнал; (в) — результат удаления шума классическим методом; (г) — результат удаления шума предлагаемым методом

А.

Рис. 2. Проверка гипотезы на гетероскедастичность шума.

а — вейвлет-коэффициенты первого уровня декомпозиции; б — выборочные дисперсии Д*

3

2

б 1

-1 -2

100

200

300

400

500

600

700

800

900

1000

к

1_и_Ш_I_1_

100

200

300

400

500

600

700

800

900

1000

Рис. 3. Разбиение на кластеры вейвлет-коэффициентов: первого и третьего (б) уровней декомпозиции после уточнения границ

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

ния с порогами по формуле (7) выполнялось обратное МНДВП по формуле (4) для вычисления

оценки полезного сигнала /(ti). На рис. 1 приведены оценки функции /(ti), полученные двумя

а

способами: путем удаления шума классическим методом (график (в)) с помощью общего порога и предлагаемым методом (г) с помощью порога, вычисленного по кластерам. Из рис. 1 (в) видно, что участки сигнала, соответствующие самым за-шумленным кластерам не очистились, поскольку общий порог оказался для данных кластеров ниже необходимого уровня. График (г) на рис. 1 показывает, что за счет вычисленного порога для каждого кластера удалось значительно подавить шум.

ЗАКЛЮЧЕНИЕ

Предложенный метод очистки от шума путем кластеризации ВК расширяет возможности порогового удаления шума и позволяет значительно эффективнее удалять локальный шум по сравнению с традиционными методами очистки от шума — практически полностью.

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант 06-01-00457-а).

СПИСОК ЛИТЕРАТУРЫ

1. Donoho D., Johnstone I. Adapting to Unknown

Smoothness via Wavelet Shrinkage // J. Amer. Statist. Assoc. 1995. V. 90. P. 1200-1224.

2. Vidakovic B. Statistical Modeling by Wavelets. New York: John Wiley & Sons, 1999. 382 p.

3. Raw lings J., Pantula S., Dickey D. Applied Regression Analysis. Spinger, 1989. 658 p.

4. Percival D., Walden A. Wavelet Methods for Time Series Analysis. London: Cambridge University Press, 2000. 594 p.

5. Кобзарь. А.И. Прикладная математическая статистика. М.: ФИЗМАТЛИТ, 2006. 816 с.

6. Inclan C., Tiao G. Use of Cumulative Sums of Squares for Retrospective Defection of Change of Variance // JASA. 1994. V. 427. P. 913-935.

7. Strang G., Nguyen T. Wavelets and Filter Banks. Boston: Wellesley-Cambridge Press, 1996. 512 p.

ОАО "Научно-инженерный центр Санкт-Петербургского государственного электротехнического университета" (Князева Т.Н., Орешко Н.И.)

Институт аналитического приборостроения РАН, Санкт-Петербург (Новиков Л.В.)

Материал поступил в редакцию 18.04.2008.

REMOVAL OF NON-STATIONARY NOISE FROM EXPERIMENTAL DATA

T. N. Knayzeva1, L. V. Novikov2, N. I. Oreshko1

1 Research & Engineering Center of Saint-Petersburg Electrotechnical University 2Institute for Analytical Instrumentation RAS, Saint-Petersburg

In this paper we present two-step method for denoising. The noise variance is changed in random moments of time. The first step of our method consists of clustering wavelet coefficients which are computed by maximal overlap discrete wavelet transform. Then one defines of noise parts with piecewise constant variance. The second step consists of calculation of the modified thresholds for each part and removing noise on basis of these thresholds. The approach practically completely allows to remove local noise in comparison with traditional methods of denoising in which the threshold calculate on all sample of data.