1. На главную
  2. Электронные версии
  3. Доклады Российской академии наук. Физика, технические науки
  4. T. 493, № 1, 2020

Доклады Российской академии наук. Физика, технические науки, 2020, T. 493, № 1, стр. 84-89

ВАРИАЦИОННЫЙ МЕТОД РАЗДЕЛЕНИЯ АОРТАЛЬНОЙ И ЛЕГОЧНОЙ КОМПОНЕНТ ЗВУКА ВТОРОГО ТОНА СЕРДЦА

А. И. Корольков 1, *, В. Г. Андреев 1, В. В. Грамович 2, А. М. Алеевская 2, Т. В. Мартынюк 2, академик РАН О. В. Руденко 1

1 Московский государственный университет имени М.В. Ломоносова
Москва, Россия

2 Национальный медицинский исследовательский центр кардиологии Минздрава России
Москва, Россия

* E-mail: korolkov@physics.msu.ru

Поступила в редакцию 29.05.2020
После доработки 29.05.2020
Принята к публикации 02.06.2020

Полный текст (PDF)

Аннотация

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

Ключевые слова: вариационный метод, легочная гипертензия, звук второго тона сердца, аортальная и легочная компоненты, давление в легочной артерии

Звук второго тона сердца S2, выслушиваемый при аускультации, может быть использован для разработки новых методов неинвазивной диагностики легочной гипертензии [1, 2]. Сигнал S2 состоит из аортальной A2 и легочной P2 компонент. В норме S2 содержит несколько колебаний и имеет длительность порядка 60 мс. Легочная гипертензия (ЛГ) в ряде случаев приводит к некоторой задержке времени закрытия легочного клапана и расщеплению S2 на компоненты. У таких пациентов наблюдается корреляция между длительностью задержки между компонентами и величиной давления в легочной артерии (ДЛА) [1–3]. В настоящее время известно несколько методов определения временной задержки между компонентами A2 и P2. В [3–5] развиваются методы частотно-временного анализа, в [6] предлагается использовать статистический метод независимых компонент, в [7] развивается оптимизационный подход. Недостатком частотно-временных методов является уширение спектральных компонент сигнала S2 при его выделении из записанной фонокардиограммы во временном окне конечной длительности. При использовании преобразования Вигнера–Вилла в частотно-временном представлении сигнала появляются интерференционные члены, минимизация которых представляет достаточно трудную задачу [8, 9]. Метод независимых компонент базируется на довольно жестком и по существу недоказанном предположении о статистической независимости легочной и аортальной компонент. Настоящая работа посвящена развитию вариационного метода разделения сигнала звука второго тона сердца на компоненты путем аппроксимации их модельными сигналами. Идея этого метода была предложена в [7]. Каждая компонента сигнала S2 моделируется гармоническим сигналом конечной длительности с экспоненциальной огибающей. Свободными параметрами модели, подлежащими оптимизации, являются амплитуда, показатель экспоненты, длительность, частота и фаза. Параметры оптимизируются в два этапа. На первом этапе вычисляются параметры огибающих, на втором – частота и начальная фаза каждой из компонент. Задержка между компонентами вычисляется также двумя способами: по параметрам огибающих функций и путем определения максимума корреляции между исходным сигналом второго тона и восстановленными компонентами.

В качестве модельного сигнала S2(t) используем сигнал, временной профиль которого был предложен в [3]. Сигнал звука второго тона S2(t) моделируется как сумма аортальной компоненты ${{A}_{2}}(t)$, легочной компоненты ${{P}_{2}}(t)$ и шума n(t):

(1)
${{S}_{2}}(t) = {{A}_{2}}(t) + {{P}_{2}}(t) + n(t).$

Компоненты имеют вид чирпированных импульсов с огибающей E0(t):

(2)
${{A}_{2}}(t) = {{\alpha }_{1}}{{E}_{0}}(t)\sin (2\pi (24.3t + 451.4\sqrt t )),$
(3)
$\begin{gathered} {{P}_{2}}(t) = {{\alpha }_{2}}{{E}_{0}}(t - {{t}_{0}}) \times \\ \, \times \sin (2\pi (21.83(t - {{t}_{0}}) + 356.34\sqrt {t - {{t}_{0}}} )), \\ \end{gathered} $
(4)
${{E}_{0}}(t) = \left[ {1 - \exp \left( { - \frac{t}{8}} \right)} \right]\exp \left( { - \frac{t}{{16}}} \right)\sin \left( {\frac{{\pi t}}{{60}}} \right).$

Здесь t – время в миллисекундах, t0 – время задержки между импульсами, α1, 2 – действительные числа порядка единицы. Функция E0(t) представляет собой огибающую, которая отлична от нуля при 0 < t < 60 мс. Величина времени задержки t0 может варьироваться в пределах 10–70 мс. При t0 < 30 мс компоненты визуально не разделяются, и S2 выглядит как один импульс. Мгновенная частота на фронте импульсов достигает значений 180–220 Гц, затем в течение первых 15–20 мс частота быстро уменьшается до значений 70–80 Гц, после чего медленно спадает до уровня 45–50 Гц. Такое поведение частоты соответствует теоретической модели, где легочный клапан, моделируемый упругой натянутой мембраной в жидкости, возбуждался коротким импульсом, возникающим при резком закрытии клапана [10]. Шум $n(t)$ в модельном сигнале считался полосовым с нулевым средним значением. Он присутствует во всей полосе частот полезного сигнала (40–220 Гц). Отношение сигнала к уровню шума (С/Ш) вычислялось по формуле

${{{{\text{С}} \mathord{\left/ {\vphantom {{\text{С}} {\text{Ш}}}} \right. \kern-0em} {\text{Ш}}} = \int\limits_{ - \infty }^\infty {{{{\left( {{{A}_{2}}(t) + {{P}_{2}}(t)} \right)}}^{2}}} dt} \mathord{\left/ {\vphantom {{{{\text{С}} \mathord{\left/ {\vphantom {{\text{С}} {\text{Ш}}}} \right. \kern-0em} {\text{Ш}}} = \int\limits_{ - \infty }^\infty {{{{\left( {{{A}_{2}}(t) + {{P}_{2}}(t)} \right)}}^{2}}} dt} {\int\limits_{ - \infty }^\infty {{{{\left( {n(t)} \right)}}^{2}}dt} }}} \right. \kern-0em} {\int\limits_{ - \infty }^\infty {{{{\left( {n(t)} \right)}}^{2}}dt} }}.$

Регистрация сигналов фонокардиограмм (ФКГ) и электрокардиограмм (ЭКГ) проводилась врачами ФГБУ “НМИЦ кардиологии” у пациентов с ЛГ различной степени тяжести. Для записи сигналов был разработан и создан комплект аппаратуры, на использование которого в условиях клиники было получено разрешение Этического комитета кардиоцентра. В соответствии с общепринятыми правилами проведения биомедицинских исследований, у всех пациентов было получено письменное информированное согласие на участие в настоящем исследовании.

В качестве приемника акустических сигналов использовался медицинский стетоскоп с мембраной, которая соединялась штатным звукопроводом с микрофоном. Сигнал микрофона усиливался и поступал на один из входов осциллографа. На второй вход осциллографа подавался сигнал ЭКГ. Оцифрованные сигналы ФКГ и ЭКГ с выхода осциллографа поступали на персональный компьютер для записи и последующей обработки. Сигналы ФКГ регистрировались в четырех положениях на грудине на выдохе при задержке дыхания. Длительность записи ФКГ составляла 12–15 с и содержала около 10 полных сердечных циклов, которые в дальнейшем обрабатывались. Для пациентов с ЛГ проводилась процедура катетеризации правых отделов сердца (КПОС) и определялось давление в правом желудочке и легочной артерии. Как правило, измерения ФКГ проводились в тот же день, что и процедура КПОС, что позволяло считать известной величину давления в легочной артерии. По записанной ЭКГ определялась длительность сердечного цикла Tс как интервал между соседними R-зубцами на ЭКГ. Считалось, что сигнал S2 начинается спустя 0.3Tс после R-зубца. Длительность выделяемого сигнала выбиралась равной 150 мс, что немного превышает типичную длительность S2, однако такой интервал позволял учесть физиологические особенности пациентов.

Алгоритм реализуется в несколько этапов.

1. Выделенный из ФКГ сигнал второго тона S2 длительностью 150 мc пропускается через полосовой фильтр с полосой пропускания в диапазоне 40–220 Гц, после чего нормируется на единицу. Получается сигнал, который обозначим $S_{{\text{2}}}^{{\text{m}}}(t)$.

2. Огибающая ${{E}_{{\text{s}}}}(t)$ вычисляется с использованием преобразования Гильберта:

(5)
$\begin{gathered} {{E}_{{\text{s}}}}(t) = {\text{|}}S_{{\text{2}}}^{{\text{m}}}(t) + i \cdot H(t){\text{|}}, \\ H(t) = \frac{1}{\pi }{\text{v}}{\text{.p}}{\text{.}}\int\limits_{ - \infty }^{ + \infty } {\frac{{S_{{\text{2}}}^{{\text{m}}}(\tau )}}{{t - \tau }}d\tau } , \\ \end{gathered} $
где i – комплексная единица.

3. Огибающая аппроксимируется суммой функций

${{E}_{{{\text{A2}}}}}(t) = {{\alpha }_{{{\text{A2}}}}}\left[ {1 - \exp \left( { - \frac{{(t - {{t}_{{{\text{A2}}}}})}}{\gamma }} \right)} \right] \times $
(6)
$\begin{gathered} \, \times \exp \left( { - \frac{{(t - {{t}_{{{\text{A2}}}}})}}{\beta }} \right)\sin \left( {\frac{{\pi (t - {{t}_{{{\text{A2}}}}})}}{{\delta t}}} \right), \\ {{E}_{{{\text{P2}}}}}(t) = {{\alpha }_{{{\text{P2}}}}}\left[ {1 - \exp \left( { - \frac{{(t - {{t}_{{{\text{P2}}}}})}}{\gamma }} \right)} \right] \times \\ \end{gathered} $
$\, \times \exp \left( { - \frac{{(t - {{t}_{{{\text{P2}}}}})}}{\beta }} \right)\sin \left( {\frac{{\pi (t - {{t}_{{{\text{P2}}}}})}}{{\delta t}}} \right).$

Функции ${{E}_{{{\text{A2}}}}}(t)$ и ${{E}_{{{\text{P}}2}}}(t)$ обращаются в ноль при $t < {{t}_{{{\text{A2}},{\text{P}}2}}}$ и $t > {{t}_{{{\text{A2}},{\text{P}}2}}} + \delta t$.

Решается задача оптимизации функционала

(7)
$\begin{gathered} Q({{\alpha }_{{{\text{A2}}}}},{{\alpha }_{{{\text{P}}2}}},\gamma ,\beta ,{{t}_{{{\text{A2}}}}},{{t}_{{{\text{P}}2}}},\delta t) = \\ \, = \left\| {{{E}_{{\text{s}}}}(t) - {{E}_{{{\text{A2}}}}}(t) - {{E}_{{{\text{P2}}}}}(t)} \right\|_{2}^{2} \\ \end{gathered} $
по семи параметрам, которые определяют вид функций EA2(t) и EP2(t). Здесь $||\cdot|{{|}_{2}}$ обозначает евклидову норму. Для уменьшения размерности множества допустимых параметров вводятся следующие ограничения на диапазон изменения параметров:

(8)
$\begin{gathered} 1\;{\text{мc}} < \gamma < 20\;{\text{мc,}}\quad 2\;{\text{мc}} < \beta < 32\;{\text{мc,}} \\ 50\;{\text{мc}} < \delta t < 70\;{\text{мc}}{\text{.}} \\ \end{gathered} $

Выбор диапазона допустимых параметров был сделан на основании значений в модели (1)–(4), а также по результатам обработки сигналов, измеренных у пациентов. Известно, что функционал Q является невыпуклым по переменным ${{t}_{{{\text{A2}}}}}$ и ${{t}_{{{\text{P2}}}}}$ [11]. Иными словами, при оптимизации определяется лишь локальный минимум функционала, и, следовательно, результат значительно зависит от начальных значений ${{t}_{{{\text{A2}}}}}$ и ${{t}_{{{\text{P2}}}}}$. Для преодоления указанной трудности задача минимизации решается несколько раз для различных пар ${{t}_{{{\text{A2}}}}}$ и ${{t}_{{{\text{P2}}}}}$. Результат, соответствующий минимальному значению Q, принимается за абсолютный минимум. Численные эксперименты показали, что приемлемые результаты получаются при минимизации функционала с 45 различными парами начальных значений ${{t}_{{{\text{A2}}}}}$ и ${{t}_{{{\text{P2}}}}}$, покрывающих область $({{t}_{{{\text{A2}}}}},{{t}_{{{\text{P2}}}}})$ в диапазоне от 0 до 100 мс c шагом 10 мс. Для каждой пары значений ${{t}_{{{\text{A2}}}}}$ и ${{t}_{{{\text{P2}}}}}$ оптимизация производится нелинейным методом наименьших квадратов [11].

4. Вычисляется задержка между легочной и аортальной компонентами по огибающей:

(9)
${{t}_{0}}_{{{\text{env}}}} = {{t}_{{{\text{P2}}}}} - {{t}_{{{\text{A2}}}}}.$

Отметим, что полученная таким образом задержка полностью игнорирует фазовый состав сигнала $S_{2}^{{\text{m}}}$.

5. На данном этапе вычисляется задержка между компонентами с использованием фазы аортальной и легочной компонент. Фаза аппроксимируется линейной функцией

(10)
${{\phi }_{{{\text{A2}}}}} = 2\pi {{f}_{{{\text{A2}}}}}t + {{\theta }_{{{\text{A2}}}}},\quad {{\phi }_{{{\text{P2}}}}} = 2\pi {{f}_{{{\text{P2}}}}}{\kern 1pt} t + {{\theta }_{{{\text{P2}}}}},$
т.е. компоненты ищутся в виде

(11)
$\begin{gathered} A_{{\text{2}}}^{{\text{r}}}(t) = {{E}_{{{\text{A2}}}}}(t)\sin (2\pi {{f}_{{{\text{A2}}}}}t + {{\theta }_{{{\text{A2}}}}}), \\ P_{{\text{2}}}^{{\text{r}}}(t) = {{E}_{{{\text{P2}}}}}(t)\sin (2\pi {{f}_{{{\text{P2}}}}}t + {{\theta }_{{{\text{P2}}}}}). \\ \end{gathered} $

Средняя частота ${{f}_{{{\text{A2,P2}}}}}$ определяется по максимальному значению фурье-компоненты сигнала $S_{{\text{2}}}^{{\text{m}}}$, умноженного на огибающую E1, 2(t):

(12)
${{f}_{{{\text{A2,P2}}}}} = \mathop {\arg \max }\limits_f \left| {\int\limits_{ - \infty }^\infty {{{E}_{{{\text{A2,P2}}}}}(t)S_{{\text{2}}}^{{\text{m}}}(t)\exp ( - i \cdot 2\pi ft)dt} } \right|.$

Начальная фаза ${{\theta }_{{{\text{A2,P2}}}}}$ определяется путем решения следующей однопараметрической оптимизационной задачи:

(13)
$\begin{gathered} {{\theta }_{{{\text{A2,P2}}}}} = \\ \, = \mathop {\arg \max }\limits_{\theta \in (0,2\pi )} \left| {\int\limits_{ - \infty }^\infty {{{E}_{{{\text{A2,P2}}}}}(t)S_{{\text{2}}}^{{\text{m}}}(t)\sin (2\pi {{f}_{{{\text{A2,P2}}}}}t + \theta )dt} } \right|. \\ \end{gathered} $

6. Вычисляется задержка между компонентами корреляционным методом. Для этого вычисляются положения максимумов кросскорреляционных функций между исходным сигналом и восстановленными компонентами:

(14)
$\begin{gathered} {{{\bar {t}}}_{{{\text{A2}}}}} = {\text{max}}\left( {\int\limits_{ - \infty }^\infty {(S_{2}^{{\text{m}}}(t) - P_{{\text{2}}}^{{\text{r}}}(t))A_{{\text{2}}}^{{\text{r}}}(t + \tau )d\tau } } \right), \\ {{{\bar {t}}}_{{{\text{P2}}}}} = {\text{max}}\left( {\int\limits_{ - \infty }^\infty {(S_{2}^{{\text{m}}}(t) - A_{{\text{2}}}^{{\text{r}}}(t))P_{{\text{2}}}^{{\text{r}}}(t + \tau )d\tau } } \right). \\ \end{gathered} $

Время задержки по фазе вычисляется по формуле

(15)
${{t}_{{\text{0}}}}_{{{\text{phi}}}} = {{\bar {t}}_{{{\text{P2}}}}} - {{\bar {t}}_{{{\text{A2}}}}}.$

В результате мы имеем две оценки для времени задержки, полученные по формулам (9) и (15). Величина задержки ${{t}_{{{\text{0phi}}}}}$, рассчитанная по фазе, не должна значительно отклоняться от значения ${{t}_{{{\text{0env}}}}}$. Значительное несоответствие между ${{t}_{{{\text{0phi}}}}}$ и ${{t}_{{{\text{0env}}}}}$ является признаком некорректности работы алгоритма.

РЕЗУЛЬТАТЫ

На рис. 1 приведен результат восстановления модельного сигнала при задержке t0 = 30 мс, при которой компоненты визуально не разделяются. Отношение С/Ш в модельном сигнале равнялось 10. Несмотря на значительный уровень шума на частотах полезного сигнала, аппроксимирующий сигнал достаточно хорошо соответствует исходному модельному сигналу. Профили компонент также хорошо восстановились, а рассчитанные задержки ${{t}_{{{\text{0env}}}}} = {{t}_{{{\text{0phi}}}}}$ = = 33.9 мс, что на 3.9 мс (13%) больше заданного значения.

Рис. 1.

Модельный сигнал (сплошная линия) и результат его аппроксимации вариационным методом (точки) (а). Вид восстановленных аортальной (сплошная линия) и легочной (штриховая) компонент (б).

Для оценки погрешности предложенного метода были проведены численные расчеты по определению времени задержки при различных уровнях шума. При типичном для измеренных ФКГ отношении С/Ш = 20 задержки, превышающие 20 мс, определяются с погрешностью порядка 10–12%. Для задержек 10–20 мс погрешность возрастает до 20–25%. С ростом отношения С/Ш погрешность практически линейно уменьшается до уровня 6–8%.

Предложенный метод был использован для разделения на компоненты сигналов S2, записанных у пяти пациентов с ЛГ. На рис. 2 представлен результат аппроксимации измеренного сигнала и восстановленные компоненты для пациента со средним давлением в легочной артерии (ср. ДЛА) 69 мм рт.ст. При таком высоком значении ср. ДЛА две компоненты можно выделить визуально. Измеренное время задержки между компонентами составило 54.1 мс. Значения времени задержки у пяти пациентов возрастало с ростом ср. ДЛА, начиная с 25.8 мс (ср. ДЛА = 24 мм рт. ст.). Задержки, определенные по огибающим компонент и их фазам, практически совпали.

Рис. 2.

Измеренный сигнал (сплошная линия) и результат его аппроксимации вариационным методом (точки) (а) у пациента с ЛГ со ср. ДЛА = 69 мм рт. ст. Вид восстановленных аортальной (сплошная линия) и легочной (штриховая) компонент (б).

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Основным достоинством предложенного вариационного метода разделения сигнала второго тона на компоненты является относительная простота алгоритма и возможность автоматизации измерений. Он достаточно устойчив по отношению к уровню шума и вариациям основных параметров самого сигнала. Существенным ограничением является рост погрешности для задержек менее 15 мс. Это связано с недостаточно хорошей аппроксимацией высоких частот на фронтах компонент в первые 15–20 мс. Это время составляет менее четверти от длительности основного низкочастотного сигнала. Но именно большая продолжительность сигнала с практически одинаковой частотой обеспечивает хорошую аппроксимацию всего сигнала. Дополнительным результатом, который также можно использовать в качестве диагностического признака ЛГ, является средняя частота легочной компоненты сигнала ${{f}_{{{\text{P2}}}}}$. Как показано в [10], повышение давления в ЛА и связанный с этим рост натяжения ее упругих стенок приводит к повышению средней частоты колебаний стенок при их импульсном возбуждении. Этот эффект, в свою очередь, должен приводить к увеличению ${{f}_{{{\text{P2}}}}}$, которая вычисляется в процессе реализации алгоритма. Значения ${{f}_{{{\text{P2}}}}}$, вычисленные у пяти пациентов с ЛГ, не превышают 45 Гц, и они, как правило, ниже средней частоты аортальной компоненты ${{f}_{{{\text{A2}}}}}$. Для пациентов со ср. ДЛА, превышающим 55 мм рт. ст., наблюдается увеличение ${{f}_{{{\text{P2}}}}}$ до 53 Гц, что превышает значение ${{f}_{{{\text{A2}}}}}$.

Список литературы

  1. Алмазов В.А., Салимьянова А.Г., Шляхто Е.В. Аускультация сердца. СПб: Издательство СПбГМУ, 1996. 458 с.

  2. Abrams J. Current Concepts of the Genesis of Heart Sounds I. First and Second Sounds // JAMA. 1978. V. 239. № 26. P. 2787–2789.

  3. Xu J., Durand L.-G., Pibarot P. Extraction of the Aortic and Pulmonary Components of the Second Heart Sound Using a Nonlinear Transient Chirp Signal Mo-del // IEEE Transactions on Biomedical Engineering. 2001. V. 48. № 3. P. 277–283.

  4. Allen J.B., Rabiner L.R. A Unified Approach to Short-Time Fourier Analysis and Synthesis // Proc. IEEE. 1977. V. 65. № 11. P. 1558–1564.

  5. Boudreaux-Bartels G.F., Parks T.W. Time-varying filtering and signal estimation using Wigner distribution synthesis techniques. IEEE Trans. Acoust. Speech Signal Processing, 1986. P. 442.

  6. Nigam V., Priemer R. A Procedure to Extract the Aortic and the Pulmonary Sounds from the Phonocardiogram // Proc. 28th IEEE EMBS An. Int. Conf. 2006. P. 5715–5718.

  7. Sæderupa R.G., Hoanga P., Wintherc S., Bøttcherd M., Struijkb J., Schmidtb S., Østergaard J. Estimation of the Second Heart Sound Split Using Windowed Sinusoidal Models // Biomedical Signal Processing and Control. 2018. V. 44. P. 229–236. https://doi.org/10.1016/j.bspc.2018.04.006

  8. Mecklenbräuker W., Hlawatsch F. The Wigner distribution: theory and applications in signal processing. Elsevier Science, 1997. P. 480.

  9. Djebbari A., Bereksi-Reguig F. Detection of the valvular split within the second heart sound using the reassigned smoothed pseudo Wigner–Ville distribution // BioMedical Engineering OnLine. 2013. V. 12. № 37. http://www.biomedical-engineering-online.com/content/12/1/37

  10. Андреев В.Г., Грамович В.В., Выборов О.Н., Мартынюк Т.В., Родненков О.В., Руденко О.В. Колебания полулунного клапана, моделируемого упругой натянутой мембраной в жидкости // Акуст. журн. 2019. Т. 65. № 6. С. 647–652.

  11. Boyd S., Vandenberghe L. Convex optimization. Cambridge: Cambridge University Press, 2004. P. 742.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Доклады Российской академии наук. Физика, технические науки

Архивы выпусков Информация о журнале Отправить рукопись в журнал