CC BY
27Программные фильтры данных мониторинга артериального
давления и пульса
Гурьева В.М. (1), Котов Ю.Б. ([email protected]), Эсселевич И.А. (3)
(1)Московский областной НИИ Акушерства и Гинекологии, (2)Институт Прикладной Математики им. М.В.Келдыша РАН, (3)Московский Физико-Технический Институт (Государственный Университет)
Предлагаемая работа является частью исследования, имеющего целью разработку алгоритмов автоматического анализа суточных записей артериального давления и частоты пульса больных акушерской клиники. В этой работе сделана попытка моделирования «врачебных» методов отсеивания ошибочных наблюдений.
Первый шаг сделан на основе прямых ответов врача на вопрос «какие измерения он считает ошибочными?». Ответы врача не были ни категорическими, ни исчерпывающими. По этой причине были испытаны различные варианты алгоритмов, не противоречащие требованиям врача.
Отработка алгоритмов производилась на записях давления и частоты пульса беременных (в том числе, больных гестозом), прошедших обследование в Московском областном научно-исследовательском институте акушерства и гинекологии (МОНИИАГ). Гестозом принято называть патологическое состояние беременных, обычно сопровождаемое повышением артериального давления, наличием белка в моче и отеками. Тяжелый гестоз приводит к гибели плода или заметному отставанию в развитиии. Для женщины гестоз опасен возможной эклампсией (нарушение мозгового кровообращения, приводящее часто к гибели женщины).
Гестоз является одним из часто встречающихся осложнений беременности. В настоящее время он выявляется у 10-25% всех беременных, и эта доля не имеет тенденции к снижению. В структуре причин материнской смертности в настоящее время гестозы занимают одно из первых мест [1].
Точность и своевременность диагностики гестоза оказывается недостижимой при использовании обычных эпизодических измерений артериального давления.
Кровь в сосудах кровеносной системы человека находится под избыточным давлением относительно атмосферы. Величина этого давления различна в разных сосудах и непостоянна во времени. По мере продвижения по большому кругу кровообращения от аортального (выходного) клапана сердца по артериям уменьшающегося калибра и далее по венам до впадения в правое предсердие давление систематически понижается от максимального (порядка 100-150 мм.рт.ст.) до практически нулевого. При этом, в соответствии с сокращениями сердца, давление пульсирует от максимального (систолического)
до минимального (диастолического). Наибольшие величины этих давлений можно наблюдать в аорте, рядом с сердцем. Практически давление принято измерять на крупных артериях, где оно отличается от аортального лишь незначительно.
Для измерения артериального давления обычно применяется неинвазив-ный метод, т.е. датчик измерительного прибора не вводится внутрь кровеносного сосуда, а давление оценивается путем определенных манипуляций вне организма.
Результат измерения зависит от многих факторов, например, от реакции пациентки на саму процедуру измерения. При эпизодическом измерении нередко регистрируются показатели, далекие от характерных для данной больной.
В последние годы в кардиологической практике используется метод суточного мониторирования, заключающийся в том, что в течение 24-28 ч. с интервалом 15-30 минут у пациента осуществляется периодическое измерение артериального давления и частоты пульса при помощи носимого пациентом портативного аппарата (монитора давления и частоты пульса).
В используемом нами варианте метода датчик прибора измеряет давление в манжетке и регистрирует осцилляции давления, возникающие в сосуде. Когда сосуд пережат, потока крови нет, нет и осцилляций. Начинаем снижать давление в манжетке. С какого-то момента кровь начинает впрыскиваться дискретными порциями в те интервалы, когда давление в сосуде превышает давление в манжетке. Появляются осцилляции. Считается, что в этот момент измеряемое давление соответствует систолическому (верхнему) давлению в сосуде. Снижаем давление в манжетке дальше. Пульсации происходят за счет дозаполнения, пока стенка сосуда не до конца растянута. Наконец, когда сосуд заполнен, он становится слишком жестким, чтобы реагировать на пульсовые перепады давления. Осцилляции пропадают. Считается, что в этот момент давление воздуха в манжетке соответствует диастолическому (нижнему).
В исследовании были использованы аппараты фирмы "МЫйеЛ" модели АВРМ-02 и АВРМ-02о. Прибор может делать измерения как в автоматическом режиме, так и в ручном. Для работы в автоматическом режиме его программируют заранее на определенное время исследования. Можно задавать различные частоты измерений в дневное и ночное время. Мы выбирали интервал между измерениями, равный 15-ти минутам днем и 30-ти ночью.
Данные, полученные в ходе мониторирования, накапливаются в памяти прибора. По окончании измерений их переносят в базу данных АВРМВАББ, обслуживающую аппарат. Для дальнейшей обработки их экспортируют в базу формата dBASE.
На рис.1 представлены два примера диаграммы, выданные программой ABPMBASE по результатам суточного мониторирования двух пациенток (Ииш = 107 и Ииш = 160). По оси абсцисс - время суток (часы), по оси орди-
нат - систолическое и диастолическое давление (мм. рт. ст.) и частота пульса (уд/мин). Шкалы общие, поскольку численные диапазоны изменения всех трех величин совпадают. Ломаная пунктирная линия - ЧСС (частота сердечных сокращений), вертикальные столбики - АД (артериальное давление). Нижний конец этого столбика представляет нижнее (диастолическое) давление, а верхний - верхнее (систолическое) давление. Кружками отмечены внеплановые измерения. Приподнятая горизонтальная пунктирная линия соответствует дневному времени, сниженная - ночному. Ночные и дневные измерения маркировались согласно дневнику самонаблюдения беременной.
170 150
100
а) №лп=107
>: время |час.) У: сист. АД, днаст. АД Qmu.pi.ci>, ЧСС (мин)
170 150
100
8 11 14 17 20
Ь) №Ш1=160
170 150
100----И" ^
й б ре рал (час.) V: сист. АД, д и ас т. ЛЦ (ым.рт.ст). ЧСС (Ллнн)
170 150
100
12 15 18 21
Рис.1. Примеры записей
В ночное время величины давлений и частота пульса понижаются. Обычно это благоприятный признак. Кроме того в моменты времени, отмеченные числами «1» и «2», зарегистрированы кратковременные урежения пульса, возможно, вызванные помехами. Отметим также резкое снижение давления вблизи точек «3» и «4». По мнению врача эти «неправильные» значения не характеризуют состояние организма, и их не следует принимать во внимание при автоматическом анализе. Наша цель состоит в том, чтобы, по возможности, исключить «неправильные» точки на предварительном этапе, до начала содержательного анализа наблюдений.
В рамках данного исследования мы ограничены только информацией о давлении и частоте пульса в течение суток. Нам следует, по возможности, избавиться от помех, искажающих физиологические закономерности, которые мы ищем.
Эту процедуру отделения приемлемых измерений от остальных будем называть фильтром.
Вернемся к примеру реальной записи (рис.^). Давление и пульс изменяются согласованно, т.е. большему значению частоты пульса соответствуют большие значения артериального давления (систолического и диастолическо-го). Однако вблизи точек, отмеченных цифрами '1', '2', '3', согласованность нарушается.
а) 1\1ит=107
1 го
100-
80"
60.
10.
го
р_сИа о 1
о 2 о о О о : 8 ОО %со с О о
........О. . .......'9°А; '¿¿Ъ ■ 0 ■ о
о о ° 4
о 3
го 40 60 80 1 00 1 20
Ь) 1\1ипп=160
1001.......
р_сИа
90-........
80
70-
60"
50-
о О
о О
Е0 70 80 90 1 00 ^Г
Рис.2. Диаграммы рассеяния для записей рис.1
На рис. 2 приведены диаграммы рассеяния для записей рис.1. По оси абсцисс отложены значения частоты пульса (Иг), а по оси ординат - диасто-лического давления (р_<Иа). Дальнейшее изложение будет вестись в пространстве этих двух переменных. Но аналогично можно строить фильтры и для других пар переменных, например, частоты пульса и систолического давления. Каждой точке на плоскости Ьт, p_dia соответствует одно измерение из записи рис.1. Диаграмма рассеяния представляет сравнительно плотное облако. Точки, которые вызвали замечания врача ( '1', '2', '3' на рис.^ и рис^), отстоят от центра облака относительно далеко.
Будем строить фильтр, как замкнутый контур на плоскости Ьт, p_dia, окружающий плотное облако «хороших» точек и отделяющий далеко отстоящие от него «плохие» точки. Рассмотрим несколько вариантов аппроксимации этого контура.
Проще всего строить разделяющий контур как совокупность независимых интервалов по выбранным координатам. Определив для каждой переменной допустимые интервалы, мы получаем на плоскости Ьт, p_dia граничные линии для каждой из переменных. Пересечение этих границ дает прямоугольник со сторонами, параллельными координатным осям.
Для практического применения остается решить вопрос о положении границ допустимого интервала по каждой координате.
Рассмотрим распределение точек в проекции на соответствующую координату. На рис. 3 представлена гистограмма распределения по частоте пульса для больной Ыит = 107, по оси абсцисс - частота пульса, по оси ординат -доля измерений от общего числа измерений.
Рис. 3. График плотности распределения Иг
Плотность точек высока в небольшой области, вне которой лежат «хвосты», содержащие малую часть измерений. Для исключения этих «хвостов» обычно используют квантильные границы, отсекающие определенную долю наблюдений. Напомним, что для произвольной функции распределения Р(х) случайной величины X величина Ец есть ц-квантиль [2], если она является решением уравнения
Ц = Р (О . (1)
Тем самым, если выбрать положение границ так, чтобы были отсечены только "хвосты" распределения с заданными квантильными уровнями ц и 1-ц, мы получим простейший однокоординатный фильтр.
Техническая трудность этого метода состоит в том, что именно в хвостах эмпирические распределения содержат мало точек, и их локальная структура там оказывается очень грубой. Поэтому будем аппроксимировать
эмпирическую функцию распределения подходящей стандартной функцией распределения. Приблизительная симметрия эмпиричесого распределения по каждой переменной позволяет взять в качестве аппроксимирующего распределения нормальное.
В качестве оценок среднего и стандартного отклонения будем использовать выборочное среднее Е и выборочное стандартное отклонение 8.
Для распределения Иг у больной Ииш = 107 (рис.3.) находим Е = 81, 8 = 12. Аппроксимирующая нормальная кривая дана на рис.3 в виде гладкой сплошной линии.
а) 120"
100-
8060-
1020
20 АО 00 80 100
Ь)
1003080706050.
4030
60 70 80 30 100 110 [1Г
Рис.4. Пример прямоугольного фильтра
Рассмотрим результаты применения двумерного варианта этого фильт-ра.На рис. 4 показаны диаграммы рассеяния (рис.2) и контур прямоугольника фильтра. По оси абсцисс - частота пульса, по оси ординат - диастолическое давление. Выбранная квантильная доля обозначена буквой р. Параметры аппроксимирующих нормальных распределений: рис. 4a - для Иг (Е=81, 8=12), а дляр_Ша (Е = 73, 8 = 10); рис. 4Ь - для Иг (Е = 89, 8 = 10), а дляр_Ша (Е = 65, 8 = 13).
1\1ит=1 07 А=0 р=0.01 Ьг: Е = 81, г = 12
р_&а: Е = 73, г = 10
Р_с11а □ 1
□ 2 о О О о 8 оо °от о о
. .С» . . .........'¿Л' о ■ 8 °°° 0сь о°9 о о
О О о ■ ° ° 0
□ 3
Нит=160 А=0 р=0.01 Ьг:Е = 89, г = 10
р_с11а: Е = 65, э = 13
р_сйа п5
о ■ о о
о 3 4 ........0 ' ' ' ° о о о °
о ° о о о о 0 о о о ■ ,0. ° ..°...........
ос о ° О 0 о ° о 8 о 8 О : о 0 00 о о :
Р оо . . . О . Р. оо о О ОО г, 8 о ° ■ и О о . о
□ 2 1 о о 9 ■ о°°....... 8
Заметим, что на рис. 4Ь точки '3' и '4', удовлетворяющие условиям пропускания фильтра, «на глаз» кажутся более далекими, чем точки '2' и '5', этим условиям не удовлетворяющие. Это, видимо, связано с тем, что точки «облака» преимущественно лежат вдоль диагонали прямоугольника, т.е. «облако» наклонено.
Попытаемся учесть этот наклон. Для этого воспользуемся координатными осями, повернутыми относительно исходных на некоторый угол ф. Дальнейшие операции те же, что и при построении фильтра-«коридора», описанного выше. Угол ф найдем из условия минимума ^ для одной из переменных в повернутых координатах.
В случае a) минимум ^ достигается при ф=90°, т.е. стороны прямоугольника параллельны осям Ьт и p_dia. А в случае Ь) ^ достигает минимума при ф=150°, т.е. прямоугольник поворачивается на 30°. При этом получаем распределения с параметрами (Е = 12.3, ^ = 8.2) и (Е = 109, ^ = 13). Результат действия такого фильтра приведен на рис.5.
По оси абсцисс - частота пульса, по оси ординат - диастолическое давление.
Ыит = 160 ф = 30° р=0.01
р»_
3 ................ ........... \ \ о о 00 0 \ :
О \ °° 0 0 V о о 8о \ 0 э : \ ...........\........
........... ........V'.........0 ' \ 9 со X 0 1п\ % 0 0 о : У
Рис.5. Наклонный «коридор»
Видно, что результаты оказались лучше, чем при использовании «коридора» без наклона, например, точка '5' (и точка '2') не отброшены, зато оказались отброшенными '3' и '4'. Для случая Ыит = 107, точки, расположенные вблизи углов, попадают в разрешенную область. Т.е. этот фильтр имеет нежелательные «слепые зоны» в углах прямоугольника.
Выберем контур, сглаживающий углы, например, эллипс.
Положение эллипса задается координатами его центра и наклоном большой оси, а размеры - длинами осей. Поместим центр эллипса в «медианную» точку, т.е. точку, координатами которой являются медианы распределений по соответствующим переменным (med(hr), med(p_dia)). При этом новые координаты точек можно записать в виде:
\x. = hr - med(hr) 1 г ^ . (2) = p _ diai - med(p _ dia)
В такой системе координат эллипс задается параметрической формулой:
fx = a cos(t)
. (3)
^У = b sin(t)
Рассмотрим модель «нарезанных» секторов. Практического значения эта модель не имеет, но хорошо иллюстрирует некоторые понятия и приемы, используемые в последующих моделях.
Разделим полный полярный угол 360° на секторы одинаковой угловой ширины. Распределение точек в секторе по радиусу задает функция распределения F(r). Эмпирическая функция распределения F(r) определяется равенством
±h (r)
F (r) = м-
n . (4)
Г0, r < rk где Ik (r) = \- k
l1, r > rk
Здесь n - число элементов в выборке, rk - значение r для k-ого элемента выборки.
Рис.6. Эмпирические функции распределения точек по радиусу для секторов с разными осевыми направлениями
Эмпирические распределения точек по радиусу в секторах могут не совпадать. На рис.6 представлены распределения для секторов с направлениями оси, заданных значениями «осевого» угла ф, и угловой полушириной а = 20°. По оси абсцисс - полярная координата г, по оси ординат - значение эмпирической функции распределения. Угловая полуширина является свободным параметром. Она должна быть достаточно большой, чтобы имело смысл рассматривать распределение точек в секторе, и достаточно маленькой, чтобы не
потерять чувствительность к направлению оси сектора. В нашей задаче приемлемым является значение а = 20°.
Возьмем один из секторов, например, с осевым направлением ф = 20°, за «образцовый». Тогда, предполагая распределения в секторах одинаковыми с точностью до масштаба по радиусу, можно привести их к «образцовому», используя линейное преобразование по радиусу:
К= к, г,, (5)
где ф - осевое направление сектора, в котором осуществляется преобразование. На рис.7 приведены два варианта такого преобразования. Коэффициенты кф вычислялись с помощью непараметрической регрессии (а) и метода наименьших квадратов (Ь).
а)
В преобразованных координатах распределения в разных секторах оказались близки, что позволяет строить общее аппроксимирующее распределение для всех секторов диаграммы.
образованных координатах
На рис.8 представлена эмпирическая плотность распределения, г - полярная координата, / - плотность распределения. Распределение несимметрично, поэтому для аппроксимации можно выбрать модель распределения Вейбулла:
Г(г) = 1 - , Х> 0, а > 0. (6)
Здесь г - полярная координата, Г(г) - функция распределения, а, X - параметры распределения.
Эта модель распределения удобна в прикладных задачах, поскольку позволяет аппроксимировать широкий диапазон эмпирических зависимостей.
Квантильная точка гг уровня г ( 0 < г < 1) распределения задается выражением
1
' 1п(1 - г)л
"X
Для аппроксимации эмпирического распределения можно подобрать параметры распределения Вейбулла, используя регрессию (способ 1) вида
1п(- 1п(1 - Г)) = 1п(Х) + а 1п(г). (8)
Во многих случаях неплохие результаты дает подбор коэффициентов а и X по квантильным точкам эмпирического распределения (способ 2). Для вычисления достаточно знать, как минимум, две квантильные точки. Лучшие результаты дает вычисление по трем точкам (медиана т, нижний д1 и верхний д2 квартили, квартилями называются квантили уровней 0.25 и 0.75). Подставляя в (8) 3 пары значений (г = д1; Г=0.25), (г = т; Г = 0.5), (г = д2; Г = 0.75),
(7)
у
получим 3 уравнения. Из 1-ого и 3-го (для д1 и д2) находим а, из всех трех, суммируя, определяем 1п ( X ):
1.57253
а =
1п(д 2/ д1)
Х = ехр
- 0.428593 -а
1п(д1) + 1п(т) + 1п(д2) 3
(9)
В нашей задаче мы использовали второй способ, как более простой и, тем не менее, дающий удовлетворительные результаты. На рис.9 приведено распределение точек по радиусу в преобразованных координатах.
По оси абсцисс отложено значение радиуса г точки в преобразованных координатах, а по оси ординат - доля измерений со значениями г не больше заданного. Кружками отмечено эмпирическое распределение (плотность которого приведена на рис. 8), сплошной линией - аппроксимирующее его распределение Вейбулла с параметрами а = 1.37 и X = 0.026.
0.8
04-
0.2-
о^г-0'""' --"
.....ж.......
А Г
............. / &
"У.......
10
20
30
Рис.9. Эмпирическое распределение точек по радиусу в преобразованных координатах с аппроксимирующим его распределением Вейбулла
Теперь мы выбираем критический квантильный уровень (1-Ь) и возвращаемся к первоначальным координатам, получая контур фильтра. На рис. 10 представлена диаграмма рассеяния, на которую наложен контур полученного фильтра, соответствующий квантильному уровню 0.95, указанному с помощью параметра Ь = 0.05.
Ыит = 107 Ь=0.05
о / о о О О 8 ОО «О® > '.."Д....
.........Ф0г-°' "о" 8 со^Геч^ о ° °° 0с= °о 1 /
о V. о д ' ° ° О
_ О
Рис.10.Контур фильтра в модели «нарезанных» секторов
У этой модели есть недостаток - она не учитывает относительное угловое положение точек в рассматриваемых секторах, в то время, как секторы приходится брать широкими, чтобы в них попадало достаточное количество точек для построения хотя бы грубых распределений. Поэтому в дальнейшем будем использовать модель «вращающегося» сектора. Она позволяет внести поправку на угловое положение точек в секторах. Будем поворачивать сектор с постоянной угловой полушириной а = 20°, плавно меняя значение осевого угла ф, и для каждого такого сектора рассматривать распределение точек в нем по радиусу. Поскольку точек конечное число, состав набора, покрытого сектором, будет меняться только, когда сектор будет «сползать» с точек диаграммы рассеяния или «заметать» новые. Поэтому рассматривать придется конечное число распределений.
В модели «нарезанных» секторов контур фильтра испытывает разрывы по границам секторов, а мы хотим получить непрерывный контур. Вернемся к построению «эллиптического» фильтра.
В полярных координатах эллипс задается формулой:
г = Г0 + р ^2(ф-ф0), (10)
где ф0 - угол наклона большой оси эллипса, р и г0 - константы.
Для эллипса с постоянной плотностью точек внутри него линия, соединяющая одинаковые квантильные точки Сц(ф) уровня ц распределений по радиусу в разных секторах, т.е. ^ из (1), будет синусоидой с периодом 180° по значениям осевого угла секторов.
Гц=Ц[Г0 +Р С°^(ф — Ф0)]. (11)
Поэтому попытаемся линию, соединяющую одинаковые квантильные точки эмпирических распределений в секторах, аппроксимировать синусоидой.
25-
го.
15
С(ф) ц = 0.5
О ф-оо^ г
........Т ■ у*
юо гоо зоо
Рис.11. Зависимость положения медианной точки от направления ф
На рис.11 приведена зависимость Сц(ф) для д=0.5. По оси абсцисс отложено значение осевого угла сектора ф, по оси ординат - Сц(ф).
На графике заметен большой пик в диапазоне углов от 130 до 170 градусов. На диаграмме рассеяния (рис. 2a) в диапазоне углов от 110° до 190° очень низка плотность точек, и именно сюда попадают точки '1' и '2', отстоящие особенно далеко от центра облака, что и приводит к появлению пика. На рис. 11 отрезками с цифрами указан диапазон секторов в которые попадают точки '1' и '2'.
Аппроксимацию функции Сц(ф) синусоидой с периодом 180° по значениям осевого угла секторов будем производить выделением второй гармоники ряда Фурье.
Для эллипса с центром в начале координат и постоянной плотностью точек все гармоники, кроме нулевой и второй, отсутствуют. Для реальных экспериментальных данных это не обязательно. Оставим для эмпирической функции Сц(ф) только эти гармоники. Тогда получим для нее аппроксимирующее выражение:
Сц(ф) = а0 + а2 cos2ф + Ь2 sin2ф = а0 + сcos2(ф-ф0), (12)
здесь а0, а2, Ь2 - коэффициенты Фурье, с и ф0 определяются из последнего равенства. Коэффициенты Фурье вычислялись численым интегрированием функции Сд(ф).
Эксцентриситет построенного эллипса е определяется из соотношения:
8= а°+С . (13)
а0 - с
На рис.12 приведена функция Сц(ф), та же, что и на рис.11, вместе с аппроксимирующей синусоидой. Аппроксимирующие параметры: а0 = 11, а2 = 2.3, Ь2 = -0.3, с = 2.3, е = 1.53, ф0 = 4°.
Рис.12. Функция Сц(ф), д = 0.5, с аппроксимирующей синусоидой
Вернемся к описанному выше пику. Напомним, что на формирование пика большое влияние оказывают точки '1' и '2' диаграммы рассеяния (рис^), и именно их мы не хотим пускать в дальнейшую работу. В «идеальном» случае точки maxl, mini, max2, min2 будут максимумами и минимумами синусоиды. А построенная синусоида через указанные точки не проходит.
Мы видим, что выбор ц = 0.5 в качестве граничного уровня отбора точек не совсем удачен. Испытаем более высокие квантильные уровни. Таким образом, мы исключаем середину облака, где точки лежат плотно, и используем точки вблизи контура, ограничивающего облако.
Рис.13. Функции Сц(ф), ц = 0.8 и ц = 0.9, с аппроксимирующими их синусоидами
На рис.13 приведены результаты аппроксимации для квантильных уровней 0.8 и 0.9. Для ц = 0.8 получили: a0 = 18, a2 = 4.1, b2 = -2.5, c = 4.9, s = 1.75, ф0 =15.8°, а для ц=0.9 получили: a0 = 23, a2 = 5.7, b2 = -3.4, c= 6.6, s= 1.82, ф0 = 15.3°. Далее будем использовать уровень ц от 0.8 до 0.9.
Из рисунка видно, что результаты стали лучше, вклад пика в положение гармоники не так велик, и синусоида намного ближе к точкам maxl, mini,
тах2, тт2. Т.е. уровень, действительно, лучше брать выше (точки ближе к «хвосту» распределения). Но мы не можем брать порог очень высоким, поскольку в «хвостах» содержится мало точек. Это видно на рис.13 для ц=0.9, где появляется еще один резкий пик (ф ~ 290°), обусловленный точкой '3' (рис.
Итак, аппроксимируя функцию Сц(ф) синусоидой, мы получили искомый эллипс. Произведем линейное преобразование радиусов точек диаграммы, уравнивающее размеры секторов:
Г '(ф) = Г (ф)
а0 + с
а0 + с cos2(ф-ф0)
(14)
ф =ф
В приведенных координатах Г, ф' точки эллипса лежат на окружности, что позволяет построить общее распределение точек облака по радиусу. Это эмпирическое распределение будем аппроксимировать непрерывным. Как и в модели «нарезанных» секторов, в качестве аппроксимирующего мы выбрали распределение Вейбулла. Теперь, выбирая критический квантильный уровень (1-Ь) и возвращаясь к первоначальным координатам, получаем эллиптическую границу.
N11171 = 107 Ф = 15.3° м = 0.9
120
100-
40-
20 40 60 ВО 100
Рис.14. Семейство «эллиптических» фильтров
На рис. 14 представлена диаграмма рассеяния, на которую наложены контуры эллиптических фильтров, соответствующих разным квантильным уровням (1-Ь), указанным с помощью параметра Ь.
Описанный метод выбора направления осей эллипса основан на разной удаленности точек эллипса от центра в разных направлениях и безразличен к количеству точек в каждом из секторов, лишь бы их было достаточно для построения хотя бы грубых распределений. Но если предположить, что точки внутри эллипса расположены равномерно, то сектор постоянной ширины
должен «захватывать» больше точек в направлении большой оси. Описываемые дальше «эллиптические» фильтры будут учитывать это свойство.
Для сектора 5фа с осевым направлением ф и угловой полушириной а введем функцию С^г(ф,а), учитывающую и число лежащих в секторе точек диаграммы рассеяния, и удаленность их от центра облака. Припишем каждой точке диаграммы рассеяния (обозначим ее Л,) единичную массу т,. Тогда функцию СНг(ф,а) можно ввести разными способами. Например:
Скгп(ф,а) = СкгтгП (ф,а) =2 т,г,", п > 0
Л'^Фа . (15)
т1 = 1V/
Функции СНгп(ф,а) (п > 0) в разной степени учитывают удаленность точек диаграммы рассеяния, попавших в сектор, от центра облака. СНг0(ф,а) описывает только число точек диаграммы рассеяния, попавших в сектор, и может рассматриваться как аналог массы сектора. СИг1(ф,а) может рассматриваться как аналог момента сил, СНг2(ф,а) - аналог момента инерции.
Пусть мы построили аппроксимирующий эллипс, задаваемый формулой:
г (ф) = а0 + с cos2(ф-ф0). (16)
где ф0 - угол наклона большой оси эллипса, а0 и с - константы, а = а0 + с -длина большой полуоси эллипса, Ь = а0 - с - длина малой полуоси эллипса. При равномерном распределении точек внутри эллипса СИгп(ф,а) ~ [г(ф)]п+2. Будем считать, что и в нашем случае выполняется такая зависимость. Для направления ф определим функцию еп+2(ф, а):
СНгп (ф, а) + СНгп ((ф +180) mod360, а) ^
Скгп((ф + 90)mod360, а) + Скг ((ф + 270)mod360, а) ~
. (17)
и+2
(ф, а)
и+2 ,
(ф) + гп+2(ф +180)
п+2
(ф)
гп+2 (ф + 90) + гп+2 (ф + 270) гп+2 (ф + 90)
г(ф)
г (ф + 90)
п+2
вп+2(фа,
а)=еп+2, где е = а/Ь - эксцентриситет эллипса. Покажем, что в
точке ф = фо функция еп+2(ф, а) достигает максимума. Действительно:
е п+2(ф, а)
г (ф) г (ф + 90)
а0 + с cos2(ф-ф 0) с cos2(ф-ф 0)
а
max еп + (ф, а) =
п+2
п+2
а0 + с cos2(ф-ф 0) а0 + с cos 2(ф + 90 - ф 0)
1 + 2-
с cos2(ф-ф 0) , - с cos2(ф - ф0)
п+2
п+2
(18)
п+2
1 + 2-
с
а0 с у
= е п+2(ф 0, а) =
а уЬ,
п+2
= е
п+2
ф
Поэтому можно искать направление большой полуоси путем максимизации функционала еп+2(ф,а) по ф.
На рис.15 представлены графики функции еп+2(ф, а), построенные для функционалов СНг0(ф,а) (1) и СНг1(ф,а) (2) (а=20°). По оси абсцисс отложено значение угла ф в градусах, а по оси ординат значение функционала еп+2(ф, а). Для СНг0(ф,а) (1) получили ф0 = 23°, е = 1.7, для С^(ф,а) (2) ф0 = 211°-180° = 31°, е = 2.1.
е2 1) т Ыит = 107 Ф0= 23°
2) тг Ыит = 107 ф = 31°
...1............
А 1
.........„<5 . . 9 . . .....<&..........
/х/ & 6 X УМ а
о1-1-1-1—
О 100 200 300
Рис.15. Графики функции еп+2(ф, а) для разных п
С увеличением п максимумы функции еп+2(ф, а) становятся более узкими и частыми. Это можно объяснить тем, что чем выше степень у г, тем большее влияние оказывают далекие точки (с большим значением г). При очень больших п влияние сохранят только они, т.е. построение эллипса будет производиться по точкам с максимальным значением г по каждому направлению, т.е. по точкам, которые мы не сохраняем для дальнейшей обработки. Поэтому мы остановились на п от 0 до 2.
1)171 N4171=107
р_сНа о
/ - <4 о ^^"Т :8 оо °осо ° ^ А» о о8 о / / \/\ь=о.ш ""X \ь=о.оэ Ь=0.0Е
о
20 40 60 ВО 100 1 20
2) тг Мит=107
р_сйа о
О у
(1 ......'8°=ь й,45' "о': '6' '/Асч 3 ^ ^ 0<ъ о 9 о (/\ь=0.01 С\ь=0.03 О~Х,Ь=0.05 ^\Ь=0Я7
О
Рис.16. «Эллиптические» фильтры. 1) ф0 = 23°, в = 1.7. 2) ф0 = 31°, в = 2.1.
Итак, направление и эксцентриситет аппроксимирующего эллипса определены. Дальнейшие действия те же, что и при построении предыдущего эллиптического фильтра. В итоге получаем эллиптическую границу на плоскости (Иг, р_<Иа).
На рис. 16 представлена диаграмма рассеяния с наложенными на нее контурами семейств «эллиптических» фильтров для случаев (п = 0, 1). В каждом семействе представлены фильтры, отличающихся выбором квантильного уровня (1-&).
Видно, что «слепых» зон практически не осталось. И фильтры отсеивают
те точки, которые врач считает далекими.
12
На рис.17 приведен график в (ф, а) для момента 10-ой степени (ф0 = 153°, в = 2.37). Отметим, что ф = ф0 = 153° почти в точности направление на точки '1' и '2' (рис.2а), что подтверждает вышесказанное. Локальные максимумы '1' и '2' соответствуют одному и тому же наклону большой оси эллипса, поскольку отличаются на 180°.
е12 тг10 Ыит = 107 ф = 153°
2 О
1 О
о со о
Ьсссач&оооо—о—оосс с-чуссо-О— —о^оажжхочоосхххао-со
О юо 200 300 <Р
Рис.17. График функции е12(ф, а)
На рис. 18 представлена диаграмма рассеяния с наложенными на нее контурами семейств «эллиптических» фильтров для случая п = 10.
Видно, что в этом случае эллипс действительно строится по самым дальним точкам, которые мы хотим отбросить. Поэтому в дальнейшем разумно использовать моменты невысоких степеней, например, не выше второй.
тг10 Мит=107
Рис.18. Семейство «эллиптических» фильтров для случая п = 10
Этот метод дает приемлемые результаты. Основной его недостаток в том, что при определении ф0 и е используется точечная оценка (для нас не важно поведение функционала еп+2(ф, а), нам важно лишь положение точки его максимума), которая дает неплохие результаты при большом количестве данных. Но у нас измерений в одной суточной записи не так много. Поэтому хотелось бы найти метод, полнее использующий имеющуюся информацию.
Вспомним, что СИгп(ф,а) ~ [г(ф)]п+2. Тогда
п+2Скгп(ф,а) * г(ф) = г0 +рШБ2(ф-ф0). (19)
Это рассмотренный выше случай, только вместо Сц(ф) имеем [С^гп(ф,а)]1/(п+2). Будем аппроксимировать функцию [С^гп(ф,а)]1/(п+2) синусои-
дой с периодом 180° по значениям осевого угла секторов, оставляя только нулевую и вторую гармоники ее ряда Фурье:
n+2Chrn (ф, а) = a0 + a2 cos 2ф + b2 sin 2ф = a0 + c cos 2(ф - ф0) (20)
Эксцентриситет вычисляем по формуле (13).
1/2
На рис.19 приведена функция [Chr0^,a)] вместе с аппроксимирующей синусоидой. Аппроксимирующие параметры: a0 = 3, a2 = 0.21, b2 = -0.46, c = 0.50,s= 1.40, ф0 = 33°.
[Chrn(<p,a)]1/2
Рис.19. Функция [С%Го(ф,а)]1/2 с аппроксимирующей синусоидой
Видно, что аппроксимирующая синусоида хорошо отражает положение точек maxl, mini, max2, min2. Результаты действия полученного фильтра на поток исходных данных отражены на диаграмме, представленной на рис.20. Результаты фильтрации не хуже полученных предыдущим методом, но этот фильтр полнее использует имеющуюся информацию.
p_dia о
о ^о :s оо
...... /((( ........"Ь" °°° 0Сэ о°9 о V/AA ÍCXXv
^ ■ ■ ■ 'o'ó':' ■ .......
о
20 40 60 80 100 hr
Рис.20. Семейство «эллиптических» фильтров, ф0 = 33°, s = 1.40
Итак, мы рассмотрели фильтры для выделения в суточной записи давления и пульса измерений, пригодных, с точки зрения врача, для последующего анализа и диагностики. Построение и оценка фильтров производились от самых простых конструкций (покоординатные коридоры) к более сложным («эллиптические»). Такой подход позволил рассмотреть достоинства и недостатки этих вариантов и наметить пути их реализации для практического использования. Наилучшим, на наш взгляд, является последний из описанных фильтров. Он вобрал в себя все лучшее, что было в других фильтрах, уменьшив число недостатков, и не обрел при этом новых.
Чтобы не перегружать читателя большим числом различных записей, все иллюстрации делались всего для двух записей, на наш взгляд, наиболее ярко отражающих особенности каждого из построенных фильтров.
Разработанные фильтры в большей степени, чем многие использующиеся в настоящее время на практике, учитывают индивидуальность каждой записи, что особенно актуально в задачах, связанных с медицинскими исследованиями.
Литература
1) Савельева Г.М. Вестник Российской ассоциации акушеров-гинекологов. 1998, 2, 21-26
2) Холлендер М., Вульф Д.А. Непараметрические методы статистики. -М., Финансы и статистика, 1983. - 518с.
3) И.М. Гельфанд, Б.И. Розенфельд, М.А. Шифрин "Очерки о совместной работе математиков и врачей", Москва, "Наука", 1989
4) В.М. Гурьева, Л.С. Логутова, Ю.Б. Котов, В. А. Петрухин "Суточный мониторинг артериального давления и частоты сердечных сокращений при диагностике гестоза", Российский Вестник акушера-гинеколога,1,2003
5) Стентон Гланц "Медико-биологическая статистика", перевод с английского, Москва, Практика, 1999
6) Закирова Н.И. Вестник Российской ассоциации акушеров-гинекологов 1996, 4 53-54
7) Mushambi М.С., Halhgan АЖ, Wilhamson К., ВГ I Anaesth 1996, 76 133148
