Научная статья на тему 'Основной алгоритм обработки ЭКГ для повышения достоверности определения макроэлементов в крови'

Основной алгоритм обработки ЭКГ для повышения достоверности определения макроэлементов в крови Текст научной статьи по специальности «Компьютерные и информационные науки»

CC BY
445
68
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
МАКРОЕЛЕМЕНТИ / ЕКГ / ЕМПіРИЧНі МОДИ / КОМПЛЕКСНА ПЛОЩИНА / ЕТАЛОННИЙ ФРАГМЕНТ / ФАЗОВА ТРАєКТОРіЯ / МЕТОД ОБРОБКИ ЕКГ / ЦИФРОВА ОБРОБКА / МАКРОЭЛЕМЕНТЫ / ЭКГ / ЭМПИРИЧЕСКИЕ МОДЫ / КОМПЛЕКСНАЯ ПЛОСКОСТЬ / ЭТАЛОННЫЙ ФРАГМЕНТ / ФАЗОВАЯ ТРАЕКТОРИЯ / МЕТОД ОБРАБОТКИ ЭКГ / ЦИФРОВАЯ ОБРАБОТКА / MACRONUTRIENTS / EKG / EMPIRIC FASHIONS / COMPLEX PLANE / STANDARD FRAGMENT / PHASE TRAJECTORY / METHOD OF TREATMENT EKG / DIGITAL TREATMENT

Аннотация научной статьи по компьютерным и информационным наукам, автор научной работы — Яковенко И. А., Клочко Т. Р., Леус Е. А.

В статье приведены новые подходы совершенствования алгоритма обработки ЭКГ для повышения достоверности определения макроэлементов в крови. В работе рассматривается обработка данных ЭКГ с помощью метода разложения по эмпирическим модам Хуанга, что позволяет анализировать и обрабатывать нелинейные и нестационарные системы. Его отличительная черта состоит в том, что базисные функции разложения определяются адаптивно на основе обрабатываемых данных. Также в данной статье приведены подходы восстановления эталонных траекторий ЭКГ с помощью преобразования сигнала в комплексную плоскость, что значительно повышает качество съема и обработки электрокардиограмм. Авторами предложено в дальнейшем, для оценки эталонного фрагмента ЭКГ использовать методы обработки изображений, в частности применять контурную согласованную фильтрацию.

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Basic algorithm of treatment EKG for the increase of authenticity of determination of macronutrients in blood

In the article new approaches over of perfection of algorithm of treatment ЭКГ are brought for the increase of authenticity of determination of macronutrients in blood. Treatment of these EKG is in-process examined by means of method of decomposition on the empiric fashions of Huang, that allows to analyses and process the nonlinear and non-stationary systems. His distinguishing feature consists of that, that the base functions of decomposition are determined adaptive on the basis of the processed data. Also in this article approaches over of renewal of standard trajectories EKG are brought by means of the signal shaping in a complex plane, that promotes quality of output and treatment of electrocardiograms considerably. It is offered authors in future, for the estimation of standard fragment of EKG to use the methods of processing of images, in particular to apply the contour concerted filtration.

Текст научной работы на тему «Основной алгоритм обработки ЭКГ для повышения достоверности определения макроэлементов в крови»

Радіоелектроніка біомедичних технологій

РАДІОЕЛЕКТРОНІКА БІОМЕДИЧНИХ ТЕХНОЛОГІЙ

УДК 621.372.542:615.849.19

ОСНОВНОЙ АЛГОРИТМ ОБРАБОТКИ ЭКГ ДЛЯ ПОВЫШЕНИЯ ДОСТОВЕРНОСТИ ОПРЕДЕЛЕНИЯ МАКРОЭЛЕМЕНТОВ В КРОВИ

Яковенко И. А.1, аспирант,

Клочко Т. Р.1, к.т.н., с. н.с.,

Леус Е. А.2, к.м.н.

1 Национальный технический университет Украины

«Киевский политехнический институт», г.Киев, Украина

2

Украинская детская специализированная больница «ОХМАТДЕТ», г.Киев, Украина

Введение

Для оценки состояния организма с позиций его функциональной, гемодинамической сбалансированности и водно-электролитного обмена необходим контроль его жизненно важных элементов, к которым в череде первых следует отнести K+, Ca++, Na+ . При дисбалансе макроэлементов в организме необходимо постоянное отслеживание содержание этих элементов [1]. В клинической практике оно анализируется с помощью биохимического анализа крови, что, в свою очередь, является травматической манипуляцией. Особенно в педиатрической практике такой метод приводит к дополнительной физической и моральной нагрузке пациента.

В авторских работах [2,3] было предложено создание неинвазивного метода для исследования in vivo содержания макроэлементов K+, Ca++, Na+ в крови с помощью обработки данных электрокардиограмм (ЭКГ), так как электрокардиограф является одним из главных средств диагностики сердечнососудистой системы. Известно, что состояние сердечнососудистой системы напрямую связана с балансом макроэлементов в крови [4].

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

Диагностическими признаками являются амплитуды зубцов P, Q, R, S, T, комплекса QRS и временных интервалов P - Q, S - T, R - R, Q - T, T, U, что измеряются следующим образом: разделение ЭКГ на отдельные цик-

176

Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2011.-№46

Радіоелектроніка біомедичних технологій

лы (R - R интервалы); выделение наиболее представительного цикла; оценка эталона ЭКГ на основе усреднения отдельных циклов; измерение амплитудно-временных параметров, характеризующих форму объекта. Кроме того, в результате измерения биоэлектрических процессов в миокарде с помощью электрокардиографа и записи сигнала в цифровой вид, мы столкнулись с рядом проблем, таких как низкочастотный дрейф, шум и т.д., что описано в работе [2,3]. Алгоритмы и методы обычной цифровой обработки в основном предназначены для линейных стационарных систем. Последние десятилетия активно развиваются новые методы анализа нелинейных стационарных и линейных нестационарных систем [5,6].

Таким образом, целью работы является повышение достоверности неинвазивного анализа макроэлементов в крови с применением разложения по эмпирическим модам обработки ЭКГ и восстановление эталонного фрагмента ЭКГ помощью преобразования в комплексную плоскость.

Постановка задачи

В первую очередь восстановим эталонный фрагмент ЭКГ. Представим, что идеальный (эталонный) сигнал ЭКГ ye (t) является периодической

функцией времени. Наблюдаемый сигнал y (t) есть искаженная форма эталонного сигнала y(t ) = f (ye ,$). В работе будем рассматривать следующую модель наблюдаемого сигнала: y(t) = (1 + 8k)ye (t - т) + £ + h(t)

В дискретной форме

y(n)=(1 + 8k)ye (n - т)+ £ + h(n) (1)

где £ - внешнее возмущение (шум), h(n) - низкочастотный дрейф, 8k - коэффициент растяжения-сжатия, т - временная задержка.

При этом будем полагать что £ , т , 8k являются последовательностями независимых случайных величин с нулевым математическим ожиданием. Требуется по наблюдениям сигнала Y = {у(1) у(2)...y(N» восстановить эталонный сигнал y e (n).

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

Вісник Національного технічного університету України "КПІ" 177 Серія — Радіотехніка. Радіоапаратобудування.-2011.-№46

Радіоелектроніка біомедичних технологій

единицу; б) в любой точке функции среднее значение огибающих, определенных локальными максимумами и локальными минимумами, должно быть нулевым. В свою очередь СМФ представляет собой колебательный режим, но вместо постоянной амплитуды и частоты может иметь переменную амплитуду и частоту, как функции времени.

Изложение метода решения задачи Прежде всего, необходимо избавиться от низкочастотного дрейфа (частота дрейфа <0.5 Гц) и шума. Для этого можно применять как обычные КИХ или БИХ фильтры, так и использовать оптимальный фильтр Виннера [7]. Применим к отфильтрованному сигналу дискретное преобразование Гильберта. Для этого промоделируем КИХ фильтр с линейной ФЧХ, имеющий сдвиг фазы на ^, и согласующую линию задержки (Рис.1). Таким

образом, на выходе преобразователя Гильберта два сигнала, сдвинутых друг относительно друга на 90 градусов.

После преобразователя Гильберта мы получим два сигнала ys (n), yc (n), сдвинутых друг относительно друга на ^ . При этом сигнал ys (n)

определяется формулой (1). Представим сигналы ys (n) и yc(n) в виде од-

ного комплексного сигнала:

z (n )= ys(n)+ іУс(n )= ЛІфУ* (2)

Представим эталонный сигнал ye(t) в виде разложения в ряд Фурье:

да да

Ує (t) = Z Ak cos(to)+ Z Bk sin (hat)

k=0 k=0

Обозначим:

да да

yy (t) = Ує (t) = Z Ak c0s(k^t) + Z Bk sin(ka,t) (3)

k=0 k=0

Применяя к (3) преобразование Гильберта, получим:

178

Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2011.-№46

Радіоелектроніка біомедичних технологій

(jj (jj

У* (t) = -Z Ak sin(Ш) + Z Bk cos(k^t)

(4)

k=0

Функции ys (n) и yc (n) связаны с функциями yes (n) и yec (n) формулой (1). Применим к (3) и (4) преобразование (1), учитывая после фильтрации в (1) уйдут слагаемые % и h(n). Имеем:

Ґ я я \

ys (t) = (1 + 5k) Z Ak cos (кш) + Z Bk sin (k®t)

Vk=0 k=0

f я я Л

-x(1+5k) zAksin(k^t)-ZBk cos (ko>t) 1 + 4

Vk=0 k=0 J

f я я

k=0

(5)

+

(6)

yc (t) = (1 + 5k) -Z Ak sin (kwt) + Z Bk cos (kwt)

V k=0 k=0 J

яя

+t(1 + 5k) Z Ak cos(kwt) + ZBk sin(kwt) +4

Vk=0 k=0 J

Учитывая (4), выражение (5) можно представить в виде:

у* (n)=(1+5k) ye (n)- т(1+5k) yec (n)+4 Ус (n) =

= (1 + 5k)ye (n) + T(1 + 5k)y* (n ) + 4

Амплитуда комплексного сигнала (2) с учетом (6) будет определяться выражением:

A = A* (1 + Sk У1 + т2; A* =jyf+yf . (7)

Таким образом, последовательности Y = {у(1), №...y(N)} соответствует траектория в двумерном пространстве с координатами ys и Ус. За счет возмущений, последовательность Y = {у(1), у(2),...y(N)} будет уклоняться от эталонной траектории, определяемой последовательностью Y = {Уе (1), Уе (2),.. •Уе (N)}. При этом, уклонение реальной траектории от идеальной определяется выражением:

<2

5A2 = A2e ((1 + 5k)V 1 + т2 -1)

(8)

Поскольку величины Sk и т имеют нулевое математическое ожидание и независимы, то математическое ожидание отклонения SA также стремиться к нулю. На рис. 2 представлены графики наблюдаемой последовательности одного из ЭКГ во времени и в комплексной плоскости.

При усреднении во временной области необходимо очень точно выделить отдельные циклы во времени и синхронизировать их, что существен-

Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2011.-№46

179

Радіоелектроніка біомедичних технологій

но усложняет алгоритм обработки. А при обработке в комплексной (фазовой) плоскости достаточно просто провести разбиение исходной последовательности (рис. 2 а) на N интервалов, таких чтобы в них гарантировано были все информационные параметры (все зубцы Р, Q, R, S, У), причем интервалы могут не совпадать по размерности друг с другом. Далее разобьем последовательность измеренного сигнала на N интервалов для восстановления эталонного фрагмента ЭКГ. При этом, как указывалось выше, разбиение проводим таким образом, чтобы в каждом интервале находился полный фрагмент ЭКГ (необходимо выбирать длину интервала не больше расстояния между интервалом R - R).

Далее выбираем опорную траекторию в качестве оценки первого приближения. Для этого нужно определить расстояние между всеми N траекториями и выбрать траекторию с наименьшим расстоянием до остальных траекторий.

Пусть была определена опорная траектория (интервал) M0.

Далее для каждой точки m0l є M0 задаем окрестность

Лу , причем m, є Ay . Ay.

б

Рис.2. Последовательность ЭКГ а) изменение сигнала ЭКГ во времени; б) фазовая траектория в комплексной плоскости

причем m0i є лу можно определить например как окружность некоторого

радиуса ri с центром в точке m0i. После этого считаем уточненное значе-

/V

ние m0i по формуле:

m

0i

m0i =

+1

j+i

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

m

n +1

(9)

e

180

Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2011.-№46

Радіоелектроніка біомедичних технологій

Здесь m.j є Asf - все точки (любых траекторий), которые входят в

окрестность Ast, n - количество точек, входящих в окрестность dSi.

На рис.3 показаны результаты восстановления эталонного фрагмента с помощью обычного усредненного сигнала.

а

б

Рис.3. Графики восстановления эталонного фрагмента: а) фазовая и усредненная траектория; б) восстановленный цикл; в) исходный цикл ЭКГ

в

Любой

произвольный

нал 3'^) можно разделить на мейство функций внутрен-мод, придерживаясь изложенной ниже методики. Для наглядности

ем ее на примере разложе-массива данных цифрового

нала который пред-

ставляет собой сумму двух стационарных по амплитуде моник.

Рис.4. Выделение экстремумов функции

сиг-

се-

них

ру-

ния

сиг-

не-

гар-

1. Определяем все экстремумы функции (рис. 4). Через точки максимумов и минимумов проводим две огибающие, используя интерполяцию, например на основе кубических сплайнов (рис.5).

Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2011.-№46

181

Радіоелектроніка біомедичних технологій

значении верхней и среднее значение тЛ^), полу-

'ЛІЛ

(10)

2. Вычисляем текущее локальное среднее mi нижней огибающих.

3. Вычитая из исходного сигнала У чают первое приближение СМФ:

чМ = уС

4. Повторяем операции 1 и 2 и 3, принимая вместо >,(А) функцию

и находим второе приближение к первой функции моды СМФ -функцию:

h-(k) = hl(k)-m-(k) (П)

здесь - среднее значений

огибающих функции ki{k) определенных ее локальными максимумами и локальными минимумами.

Аналогично находим третье и последующие приближения к первой функции моды СМФ. С увеличением количества итераций функция ьп(к) стабилизируется, то есть стремится к неизменяемой форме. В качестве одного из критерия остановки процесса можно воспользоваться величиной

(^l) - Ьї

Рис. 5. Интерполяция экстремумов

6 =

4-і/

w

8,

(12)

где °а пороговое значение по достижению которого процесс вычисления приближений функции моды СМФ останавливается.

5. По достижению ^ порогового значения и остановки итерационного процесса последнее значение функции кп(к) может быть принято за

первую и самую высокочастотную моду СМФ: С*(Ю -

После получения первой моды вычитаем ее из исходной функции

пункты 1-5 и получаем вторую моду

6. Повторяем для остатка М с-(к). Вычитаем ее из остатка Т\-0

г:

Повторяем для остатка П(

и для всех последующих остатков эту процедуру, получим представление исходного сигнала в виде разложения по эмпирическим модам:

п

У{

ci{

ТІ

(13)

где гп - остаток на последнем этапе декомпозиции сигнала.

182

Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2011.-№46

Радіоелектроніка біомедичних технологій

Остановка декомпозиции сигнала должна происходить при максимальном «выпрямлении» остатка, т.е. превращения его в тренд сигнала по заданному интервалу. В частности, если остаток гп является монотонной функцией, то из него не может быть извлечено ни одной функции СМФ. Также процесс декомпозиции можно остановить, если остаток является несущественным по своим значениям в сравнении с исходным сигналом.

На основе разложения (13) можно осуществлять фильтрацию сигнала, исключая СМФ с соответствующими номерами и остаток, который представляет собой медленный тренд либо константу. Результатом будет сигнал вида

jp = Z CL (k), (14)

L

где L - массив индексов СМФ, которые формируют отфильтрованный сигнал.

С увеличением номера СМФ ее частота уменьшается. Первая СМФ представляет собой шум, во второй мы уже видим проявление R-зубца, а последние СМФ и остаток - низкочастотный дрейф. Получив разложение сигнала по функциям СМФ, можно убрать высокочастотный шум и низкочастотный дрейфы. Для этого достаточно просуммировать СМФ функции с индексом от 2 до 9. Результат приведен ниже на рис.6 и рис.7.

В результате, при вычитании из исходного сигнала последние СМФ функции ликвидируется низкочастотный дрейф. Высокочастотный шум при вычитании первой СМФ уменьшился в 2..3 раза. Для более сильного подавления можно вычесть вторую и третью функции СМФ (рис 8).

Этот метод позволяет без каких-либо предположений о харак-

Рис.7. Сигнал после фильтрации

Рис.8. Высокочастотный шум при вычитании второй и третей функции СМФ

Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2011.-№46

183

Радіоелектроніка біомедичних технологій

тере сигнала провести разложение последнего по ортогональному базису. При этом базис является адаптивным и определяется непосредственно по анализируемым данным.

Выводы

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

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

Литература

1. Комаров Ф.И. Биохимическое исследование в клинике/ Ф.И.Комаров, Б.Ф.Коровкин, В.В.Меньшиков. - Элиста, Джангар, 1998.-250 с.

2. Яковенко І.О. Неінвазивний моніторинг складу макроелементів в крові /Яковенко І.О., Клочко Т.Р., Леус О.О.,// ХІ Міжнародна науково-практична конференція «Людина і космос» 2009р.—Дніпропетровськ; Дніпропетровськ; НЦАОМУ, 2009. -

С.265.

3. Яковенко І.О. Визначення макроелементів K+, Ca++, Na+ У крові дитини з подальшим оцінюванням гомеостазу організму / Яковенко І.О, Клочко Т.Р., Леус О.О.// Вісник НТУУ “КПІ” серія приладобудування. - 2009. - №38.

4. Патент 38905 України А61В5/02, А61В5/0295, 2001 р.

5. Huang N.E. Hilbert-Huang transform and its applications / Norden E. Huang, Samuel

S.P. Shen. - World Scientific Publishing Co. Pt. Ltd. 5 Toh Tuck Link, Singapore,1995.

6. Huang N.E. et al. The empirical mode decomposition and Hilbert spectrum for nonlinear and nonstationary time series analysis // Proc.Roy.Soc.London, A,454.-1998, P.903 - 995.

7. Опенгейм А. В. Цифровая обработка сигналов/ Опенгейм А. В., Шафер Р. В.: Пер. с англ/ Под ред. С. Я. Шаца. - М.: Связь, 1979.-416 с., ил.

Яковенко І.О., Клочко Т.Р., Леус О.О. Основний алгоритм обробки ЕКГ для підвищення достовірності визначення макроелементів в крові. У статті приведені нові підходи вдосконалення алгоритму обробки ЕКГ для підвищення достовірності визначення макроелементів в крові. У роботі розглядається обробка ЕКГ за допомогою методу розкладання по емпіричних модах Хуанга, що дозволяє аналізувати і обробляти нелінійні і нестаціонарні системи. Його відмінна риса полягає в тому, що базисні функції розкладання визначаються адаптивно на основі оброблюваних даних. Також в цій статті приведені підходи відновлення еталонних траєкторій ЕКГ за допомогою перетворення сигналу в комплексну площину що значно підвищує якість

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

184

Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2011.-№46

Радіоелектроніка біомедичних технологій

знімання і обробки електрокардіограм. Авторами запропоновано надалі, для оцінки еталонного фрагмента ЕКГ використовувати методи обробки зображень, зокрема застосовувати контурну погоджену фільтрацію.

Ключові слова: макроелементи, ЕКГ, емпіричні моди, комплексна площина, еталонний фрагмент, фазова траєкторія, метод обробки ЕКГ, цифрова обробка.

Яковенко И.А., Клочко Т.Р., Леус Е.А. Основной алгоритм обработки ЭКГ для повышения достоверности определения макроэлементов в крови. В статье приведены новые подходы совершенствования алгоритма обработки ЭКГ для повышения достоверности определения макроэлементов в крови. В работе рассматривается обработка данных ЭКГ с помощью метода разложения по эмпирическим модам Хуанга, что позволяет анализировать и обрабатывать нелинейные и нестационарные системы. Его отличительная черта состоит в том, что базисные функции разложения определяются адаптивно на основе обрабатываемых данных. Также в данной статье приведены подходы восстановления эталонных траекторий ЭКГ с помощью преобразования сигнала в комплексную плоскость, что значительно повышает качество съема и обработки электрокардиограмм. Авторами предложено в дальнейшем, для оценки эталонного фрагмента ЭКГ использовать методы обработки изображений, в частности применять контурную согласованную фильтрацию.

Ключевые слова: макроэлементы, ЭКГ, эмпирические моды, комплексная плоскость, эталонный фрагмент, фазовая траектория, метод обработки ЭКГ, цифровая обработка.

Yakovenko I.A., Klochko T.R., Leus E.A. Basic algorithm of treatment EKG for the increase of authenticity of determination of macronutrients in blood. In the article new approaches over ofperfection of algorithm of treatment ЭКГ are brought for the increase of authenticity of determination of macronutrients in blood. Treatment of these EKG is in-process examined by means of method of decomposition on the empiric fashions of Huang, that allows to analyses and process the nonlinear and non-stationary systems. His distinguishing feature consists of that, that the base functions of decomposition are determined adaptive on the basis of the processed data. Also in this article approaches over of renewal of standard trajectories EKG are brought by means of the signal shaping in a complex plane, that promotes quality of output and treatment of electrocardiograms considerably. It is offered authors in future, for the estimation of standard fragment of EKG to use the methods of processing of images, in particular to apply the contour concerted filtration.

Keywords: macronutrients, EKG, empiric fashions, complex plane, standard fragment, phase trajectory, method of treatment EKG, digital treatment.

Вісник Національного технічного університету України "КПІ" Серія — Радіотехніка. Радіоапаратобудування.-2011.-№46

185

i Надоели баннеры? Вы всегда можете отключить рекламу.