Радиотехника и электроника, 2019, T. 64, № 12, стр. 1219-1227

1.1. Магнитокардиографические измерения

Магнитокардиография – это один из новых методов диагностики в кардиологии, который предполагает анализ величин параметров магнитного поля сердца человека, бесконтактно зарегистрированных сверхчувствительной магнитометрической аппаратурой в воздухе, над грудной клеткой пациента. Для получения диагностической информации разработаны и применяются магнитометрические системы (магнитокардиографы) на основе СКВИД. Эти системы разрабатываются на основе разных конструкций преобразователей магнитного потока детекторов поля (осесимметричные градиентометры первого или второго порядка; планарные градиентометры и т.д.) и различной конфигурции. Однако во всех случаях точками измерения магнитного поля, как правило, являются узлы прямоугольной, регулярной решетки, расположенной в плоскости, параллельной поверхности грудной клетки пациента. В каждом из этих узлов (в одном, нескольких или всех, в зависимости от используемого измерителя) синхронно регистрируются величины магнитного поля сердца и, если необходимо, опорный сигнал (одно из стандартных отведений электрокардиограммы – ЭКГ). При исследовании МКГ часто точки измерения (общим числом 36) располагают в узлах равномерной сетки с шагом 4 см (шесть рядов и шесть столбцов по взаимно перпендикулярным направлениям). При этом и для одноканальных, и для многоканальных систем характерный размер области измерений составляет ≈20 см, а плоскость располагают максимально близко к поверхности тела.

Для определенности отметим, что в нашем случае все МКГ были зарегистрированы в обычном, специально не экранированном помещении 12-канальной (девять измерительных каналов и три канала – векторный магнитометр) магнитометрической системой, а основу измерительного канала магнитометрической системы составляли СКВИД-градиентометры второго порядка с аксиальными трансформаторами магнитного потока [3]. Алгоритм измерений предполагает четыре последовательных пространственных позиции измерителя над грудной клеткой пациента. Для согласования МКГ в разных позициях синхронно выполняется запись одного отведения ЭКГ (второе стандартное отведение). Продолжительность записи в каждой позиции измерителя составляет 30 с при частоте регистрации сигнала 1000 Гц.

На этапе предварительной обработки результатов измерений производится морфологический анализ опорного ЭКГ-сигнала и определяется количество кардиокомплексов в записи, их тип (нормальный, экстрасистола, артефакт и т.д.) и местоположение на временной оси максимумов зубцов “R”. В соответствии с найденными комплексами ЭКГ, автоматически выделяются кардиокомплексы в сигналах МКГ, а затем производится их цифровая обработка, анализ и усреднение. В качестве примера на рис. 1а показан магнитокардиосигнал (после цифровой фильтрации), регистрируемый магнитокардиографом [3] с отметками найденных типов кардиокомплексов и сигнал ЭКГ. Один из режимов работы программного обеспечения магнитокардиографа позволяет сформировать временные последовательности магнитокардиосигнала в виде N = 1000 равноотстоящих моментов времени t₁, t₂, …., t_N длительностью 1 мс таким образом, что фиксированный момент времени (нашем случае, t = t₃₃₃) соответствует максимуму зубца “R” усредненного кардиокомплекса. Усредненные кардиокомплексы МКГ одного из пациентов, полученные в соответствующих точках над грудной клеткой, показаны на рис. 1б в виде графиков, а на рис. 1в – в виде “совмещенной” МКГ всех 36 точек плоскости измерений.

Полученный таким образом набор усредненных значений сигнала МКГ Ψ(1, 1, t₁), Ψ(1, 1, t₂), …, Ψ(1, 1, t_N), …, Ψ(i, j, t_k), …, Ψ(6, 6, t_N) синхронизирован как по времени, так и по пространству (в точках решетки и границах области измерений) и является исходными данными для дальнейшего анализа (i, j = 1, 6; k = 1, N).

1.2. Алгоритм преобразования данных МКГ: метод главных компонент

Итак, после регистрации и предварительной обработки магнитокардиосигнала (МКГ-сигнала) известны N × 36 значений величин выходного сигнала СКВИД-градиентометра второго порядка в каждом узле решетки измерений Ψ(1, 1, t₁), …, Ψ(1, 1, t_N),…, Ψ(i, j, t_k), …, Ψ(6, 6, t_N). Диапазон изменения значений параметра Ψ определяется как разрядностью аналогово-цифрового преобразователя, который используется для введения регистрируемого МКГ-сигнала в компьютер, так и расстоянием между датчиком магнитного поля и источником сигнала, которое зависит в том числе от особенностей анатомического строения грудной клетки пациента. Поэтому на первом этапе обработки данных производится нормировка результатов измерений: вычисляются нормированные 36 × N значений МКГ-сигнала Ω (I, t_k) таким образом, чтобы выполнялось условие:

где с, С – константы, I – номер точки наблюдения с координатами (i, j) (I = 1, …, 36).

После нормирования МКГ-сигнала сформируем матрицы ξ_а и ξ_b. Размерности матриц равны 36 × N и N × 36 соответственно. При этом элементы строки матрицы ξ_а с номером “a” совпадают со значениями нормированного МКГ-сигнала для этого момента времени, т.е. Ω(1, t_а), …, Ω(36, t_а). А элементы строки матрицы ξ_b с номером “b” – со значениями нормированного МКГ сигнала для соответствующей точки наблюдения, т.е. Ω(i, t₁), …, Ω(i, t_N). После этого введем квадратичный тензор пространственно-временных характеристик магнитного поля с размерностью 36 × 36 и вычислим значения величин его элементов с помощью свертки матриц ξ_а и ξ_b

Определенный таким образом тензор ξ(a, b) симметричен, его след равен нулю, а его элементы вещественны. Отсюда следует, что с помощью стандартных методов, например метода последовательных вращений, тензор ξ(a, b) может быть приведен к диагональному виду. Таким образом, на заключительном этапе алгоритма преобразования данных определяется вектор собственных значений λ₁, …, λ₃₆ размерностью 36 (значения диагональных элементов матрицы).

Отметим, что оператор, который преобразует тензор ξ(a, b) к диагональному виду, можно применить и для преоб6разования единичной матрицы размерностью 3 × 36. В результате получаем значения элементов α(1, 1), α(1, 2), …., α(I, J), …, α(36, 36) матрицы собственных векторов МКГ-сигнала. Затем вычислим 36 × N значения МКГ-сигнала в новой “собственной” системе координат с помощью следующих соотношений:

и автоматически сформируем последовательность β главных компонент в порядке убывания значений величин энергии МКГ-сигнала (пропорциональна сумме квадратов МКГ-сигнала для всех N моментов времени исследуемого интервала).

Временные ряды МКГ-сигнала в заданной (лабораторной) системе координат и “собственной” системе координат, базисные векторы которой совпадают с собственными векторами, связаны с помощью соотношений, аналогичных (3):

С помощью соотношений (4) можно вычислить значение МКГ-сигнала ${\delta }(i,j,k,Q)$ в узле решетки с номером i, j в момент времени k в предположении, что используются только Q первых элементов сформированной последовательности β главных компонент. Тогда точность описания данных измерений на основе соотношений (1)–(4) можно оценить с помощью значения-среднеквадратической погрешности

которое равно нулю, если параметр Q имеет максимальное значение (равно количеству точек наблюдения, в нашем случае Q = 36). Таким образом, можно утверждать, что для произвольной МКГ-записи существует предельное значение параметра Q (<36), при котором значение среднеквадратической погрешности не превышает заданный уровень (например, 0.1) для произвольного момента времени.

1.3. Численное моделирование и варианты преобразования магнитокардиосигналов

Исходные данные, полученные после цифровой обработки зарегистрированных магнитных сигналов сердца взрослого человека в виде совмещенных МКГ показаны на рис. 2а, и те же зависимости после применения алгоритма ГК-метода (1)–(4) – на рис. 2б. Как видно, исходные кривые МКГ (см. рис. 2а) содержат магнитный шум и артефакты, которые не устранены в процессе цифровой фильтрации и усреднения. Результат применения ГК-метода на рис. 2б отображает существенное повышение отношения “сигнал/шум”, что является важным для получения достоверного диагностического заключения. Порядок ослабления шумового сигнала (фильтрация) после применения ГК-метода (1)–(4) в явном виде показан на примере анализа “зашумленной” МКГ-записи в одной из точек измерения (рис. 3).

Как в первом примере (см. рис. 2), так и во втором (см. рис. 3) преобразование сигналов выполнено при Q = 5. Таким образом, для восстановления МКГ-записи можно хранить не всю исходную выборку, а только следующий набор данных:

1) тензор пространственно-временных характеристик МКГ-сигнала (размерность 36 × 36 = = 1296);

2) временные ряды МКГ-сигнала в пяти точках “собственной” системы координат (размерность 5 × N = 5000).

Другими словами, чтобы восстановить измеренные значения МКГ-сигнала в любой заданной точке измерения для произвольного момента времени необходимо хранить в памяти компьютера 6296 значений, что практически в шесть раз меньше чем исходный объем измерительной информации (6 × 6 × N = 36 000).

1.4. Магнитные карты для собственных значений

Построение и анализ динамической последовательности пространственного распределения величин параметров магнитного поля сердца пациента в границах плоскости измерений (карт магнитного поля) для выбранных моментов времени усредненного кардиоцикла – один из основных и наглядных этапов анализа МКГ. На рис. 4 показаны магнитные карты для 32 моментов времени для той же МКГ записи, что и на рис. 1. При этом временной интервал был сформирован в ручном режиме: шаг (дискретность) по времени составляет 14 мс, а границы интервала совпадают с положением узловых точек “Q” и “Te” (“Te” – точка окончания зубца “Т”) усредненного кардиокомплекса. Магнитные карты показаны в режиме заполнения соответствующим цветом промежутков между последовательностью изолиний равного уровня [4]. Отметим, что каждая карта представляет собой распределение измеренных величин магнитного поля $\Psi (\vec {r},k)$не только в заданных точках плоскости измерений (36 точек), но и в точках, которые не являются узлами измерений, а величины сигнала в этих точках получены путем применения методов интерполяции [5].

Как было показано, временные ряды МКГ-сигнала в заданной (лабораторной) системе координат и “собственной” связаны с помощью соотношений (4). Предположим теперь, что в узле решетки в “собственной” системе координат, где энергия МКГ-сигнала имеет максимальное значение, этот сигнал для всех моментов времени равен 1. Тогда с помощью соотношений (4) получим распределение значений МКГ-сигнала в узлах решетки плоскости измерений. Обозначим это распределение (в “собственной” системе координат) как G₁($\vec {r}$). Аналогично можно вычислить распределения магнитного поля G₂($\vec {r}$), G₃($\vec {r}$), …, G₃₆($\vec {r}$). Тогда схему преобразования (4) можно записать в виде

Другими словами, магнитное поле для произвольного момента времени, усредненного кардиокомплекса можно представить в виде линейной комбинации Q магнитных карт, а их количество может быть равно 1, 2, 3, … или равно количеству точек измерения (в нашем случае 36). Распределения магнитного поля G($\vec {r}$) в “собственной” системе координат можно отобразить также в виде магнитных карт, а затем с помощью (4), (5) пересчитать эти данные в заданные точки исходной плоскости измерений. На рис. 5 в графическом виде показана принципиальная схема работы алгоритма преобразования магнитометрической информации при анализе МКГ.

При реализации алгоритма ГК-метода в программном обеспечении магнитокардиографа магнитное поле для заданного момента времени кардиоцикла в границах области измерений (магнитная карта) представлено на экране компьютера в виде “графической формулы”.

Результаты преобразования данных показаны на рис. 6 для того же пациента, что и на рис. 1. Для магнитных карт в “собственной” системе координат в текстовом виде указаны значения весовых коэффициентов (рис. 6а). Погрешности между измеренным и вычисленным распределениями магнитного поля составили 0.38, 0.12 и 0.09 соответственно (в текстовом виде также указаны на рис. 6б).

1.5. Карты распределения вектора плотности токов для собственных значений

Представление результатов измерений магнитного поля в виде суммы соответствующего набора магнитных карт в математической интерпретации (6) можно рассматривать также и в соответствии с принципом аддитивности: в произвольной точке пространства вектор магнитной индукции является векторной суммой магнитных полей всех источников.

Будем предполагать, что для каждой магнитной карты рис. 6 получено решение обратной задачи магнитостатики и найдено распределение вектора плотности токов в плоскости, параллельной плоскости измерений. Метод решения обратной задачи приведен нами в [6] и основан на применении двойного интегрального преобразования Фурье. При этом система токов распределена в плоскости, которая расположена от плоскости измерений на расстоянии, равном соответствующей координате магнитного диполя (определяется по результатам измерений МКГ-сигнала на предыдущем этапе обработки данных) [7].

Исходя из этого в рамках соотношений (3) введем 36 независимых источников магнитного поля в виде систем плоских токов J₁($\vec {r}$), J₂($\vec {r}$), J₃($\vec {r}$), …, J₃₆($\vec {r}$) в собственной системе координат. При этом магнитное поле, созданное системой токов J_i($\vec {r}$), i = 1, N, в точках наблюдения совпадает с соответствующей магнитной картой G_i($\vec {r}$). Тогда для произвольного момента времени получаем следующую схему преобразования:

Иначе говоря, распределение вектора плотности токов в заданной системе координат для произвольного момента времени усредненного кардиокомплекса можно представить в виде линейной комбинации Q систем токов J_i($\vec {r}$), найденных в системе координат, связанной с собственными векторами. При этом количество карт распределения токов Q может быть равно 1, 2, 3, … или равно количеству точек измерения.

На рис. 7 для того же пациента, что и МКГ на рис. 1, показаны пространственное распределение вектора плотности токов, найденное после решения обратной задачи магнитостатики (рис. 7а) [6] и “графическая формула” для его вычисления по токам J_i($\vec {r}$) (рис. 7б).