Привет, Хабр! По своей профессиональной деятельности я занимаюсь моделированием и разработкой цифровых алгоритмов в области радиолокации. Однако универ закончил по специальности «Биотехнические и медицинские аппараты и системы», поэтому всегда хотел совместить эти два направления. И для этого как нельзя лучше подходит область биорадиолокации.

Биорадиолокация – набирающий популярность бесконтактный метод измерения жизненно важных показателей (ЖВП) человека, таких, как сердцебиение и дыхание. В отличие от контактных систем радары не нуждаются в прикреплении каких-либо датчиков на поверхность тела пациента. 

Биорадары работают на основе модуляции радиосигнала, которая возникает за счет отражения радиосигнала от тела человека. Эта модуляция обусловлена смещением грудной клетки человека, а сам отраженный сигнал содержит как дыхательные, так и сердечные составляющие, в том числе и помехи окружающей среды, а также собственные шумы приемного тракта. Шумы и помехи удаляются в приемнике при обработке сигнала, чтобы более качественно выделить показатели жизнедеятельности. Радары для обнаружения ЖВП не требуют большой мощности излучения, так как работают на небольших расстояниях (до нескольких метров), и благодаря этому они безопасны для организма. 

Этот способ диагностики не такой распространенный, как хорошо знакомые нам рентген, УЗИ, КТ, МРТ. Из способов диагностики нарушений дыхания наиболее известный и доступный метод – это проверка функции внешнего дыхания, спирометрия, когда человек делает выдох в трубочку, и замеряется объем вышедшего из его легких воздуха. 

Зачем тогда нужна биорадиолокация и где она применяется? Биорадиолокация незаменима во всех случаях, когда нельзя использовать контактные датчики:

• для измерения дыхания и сердцебиения у новорожденных,

• для измерения состояния дыхания у людей с ожогами,

• для измерения параметров дыхания при нарушении сна.

Перспективы у таких радаров есть не только в медицине, но и в охране, а также при поисковых работах МЧС – обнаружение людей под завалами.

Поэтому хочу с вами поделиться относительно простой, но в то же время полезной моделью биорадара. Модель разбита на две части. Первая часть посвящена моделированию перемещения грудной клетки человека. Вторая часть модели – про разработку FMCW-радара с последующим анализом его эффективности и применимости для обнаружения ЖВП. Итак, начнем…

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

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

В контексте радиолокации движение грудной клетки – это цель, которая совершает движения по некоему закону. Вот этот закон-то мы и будем моделировать. Особенность в том, что цель флуктуирует, то есть постоянно возвращается в исходное положение. Это перемещение, то есть расстояние, на которое движется цель, варьируется в зависимости от антропометрических показателей человека. В среднем размер грудной клетки составляет 20 см, а во время фазы вдоха-выдоха происходит изменение примерно на 5%. Таким образом, отслеживаемое перемещение составляет порядка одного сантиметра.  

Построение модели

Дыхательную систему сложно перевести в модель, особенно если учитывать все разнообразные физические процессы, происходящие в организме человека. Поэтому, как часто бывает, лучше подобрать упрощенную модель-эквивалент, которая учитывает не все параметры, а лишь то, что необходимо для решения задачи. Из курса по физическому моделированию известно, что часто такие модели представлены как системы со сосредоточенными параметрами, включающие элементы сопротивления, индуктивности и емкости, которые соответствуют различным физиологическим компонентам системы дыхания. Покопавшись в статьях и книгах (например, Physiological control systems Michael Khoo), я наткнулся на следующую схему:

Физическую интерпретацию этой схемы можно проиллюстрировать следующим рисунком:

Как видно, дыхательные пути делятся на две составляющие: большие дыхательные пути, которые моделируется с помощью сопротивления Rc, и малые Rp. Поток воздуха вызывает расширение полости грудной клетки, что моделируется последовательным соединением емкостей Cl и Cw. Также добавлена емкость Cs, которая учитывает некоторый объем воздуха, отводимый из легких. Увеличение этого объема свидетельствует о заболевании. 

Мы рассмотрели все компоненты электрической схемы рисунка 1. Далее сконцентрируемся на характерных точках схемы, которые показывают давление, возникающее на различных участках модели легких.

• Pa – давление в отверстии дыхательных путей. 

• Paw – давление в центральных дыхательных путях.  

• PA – давление в альвеолах (что-то типа воздушного мешка или воздушного пространства).

• Ppl – давление в плевральной полости (между паренхимой легких и грудной стенкой).  

Для лучшего восприятия можно рассмотреть следующий рисунок:

Эти давления соотносятся с P0, давлением окружающей среды, которое мы можем установить равным нулю. Предположим, что объемный расход воздуха, поступающего в дыхательную систему, равен Q. Тогда целью здесь является выведение математической зависимости между  Pa и Q.

Если провести аналогию между механикой и электрическим эквивалентом (рисунок 1), то поток воздуха  Q – это ток i, а давление P – напряжение U, с учетом этого из схемы рисунка 1, применяя законы Кирхгофа, получаем следующую систему уравнений:

Конечно, можно и не решать эту систему, а перейти к расчету коэффициента передачи схемы рисунка 1, но мы же не ищем легких путей. Кроме того, система уравнений наглядно показывает физику процесса. Итак, реализуем модель легких в Engee. Результат на рисунке 4.

Первое, что мы сделали, – реализовали генератор сигнала или источник, который по сути подключается между точками Pa и P0. 

Поскольку дыхание – это периодический процесс, то нам подойдет генератор синуса. Однако в ходе дыхательного процесса  возникают паузы, поэтому я взял в качестве генератора последовательность прямоугольных импульсов (это блок Pulse Generator) с добавление смещения. Период равен 4 секундам, а длительность импульса – 2 секунды, что соответствует 15 циклам вдоха-выдоха в минуту. 

Далее переходим к реализации уравнений. Двойной контур модели, который содержит 1/Cl, 1/Cw, Rp, Cs и представляет собой первое уравнение, где в качестве интеграла выступает блок Integrator, второе – это контур с Cs и производной по времени (от интеграла перешли к производной) в виде блока Derivative.

Как видно из общей модели, она параметризована, что позволяет имитировать различные патологии. Для нас важно фиксировать изменение объема легких. Результаты моделирования представлены на рисунках 5 и 6. На рис. 5 у нас получилась модель нормального дыхания, где средний объем воздуха составляет 0,5 литра. Ниже рисунок с возможной патологией. Общее время дыхания, которое мы моделируем, составляет 15 секунд.

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

Почему именно FMCW-радар?

FMCW-радар имеет несколько преимуществ при измерении перемещения грудной клетки:

• он определяет как скорость, так и дальность;

• имеет довольно высокое разрешение (вплоть до наноцелей на малом расстоянии); 

• недорогой, что позволяет приобрести уже готовое решение.  

Вкратце рассмотрим принципы работы FMCW-радар (полезные демки можно посмотреть тут). FMCW-радар – Frequency-Modulated Continuous Wave, радар непрерывного излучения с линейной частотной модуляцией (ЛЧМ), схема его работы представлена на рисунке 7.

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

Для выделения частоты дыхания данные отраженных сигналов обрабатываются с помощью быстрого преобразования Фурье. Фазовый сдвиг, вызванный дыханием, будет иметь низкочастотный компонент, обычно в диапазоне от 0.1 до 0.5 Гц (от 6 до 30 дыхательных циклов в минуту). После выделения этих частотных компонентов из спектра сигнала их можно проанализировать для определения частоты дыхания. Это достигается путем фильтрации спектра сигнала на соответствующих частотах и подсчета пиков, которые соответствуют дыхательным циклам.

Теперь приступим к моделированию…

Реализация программного кода

Для наглядности процесса разработки FMCW-радар буду строить в виде скрипта. Так как рабочим языком Engee является Julia (синтаксис можно тут посмотреть), сам код будет реализован на Julia. Также будем применять встроенные в Engee системные объекты, аналог тулбокса фазированных решеток MATLAB.

Первый этап – это инициализация входных параметров.

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

T = 15; # время моделирования
interval = 0.01 # время между измерениями.
t = Vector(0:interval:T-interval) # шаг по времени
N_ot = length(t)/T

fc = 77e9;  # задаем несущую частоту 
c = 3e8;   # скорость света
lambda = c/fc; # длина волны

range_max = 5; # максимальное расстояние до человека
tm = 5.5*range2time(range_max, c); # время развертки функция

range_res = 0.1; # разрешение
bw = rangeres2bw(range_res,c); # полоса резвертки ЛЧМ
sweep_slope = bw/tm; # крутизна пилы

fr_max = range2beat(range_max,sweep_slope,c); # верхняя частота 

maxbreathingRate = 80; # максимальное количество дыхательных циклов в минуту
forwordExpantion = 1; # изменение движения грудной клетки в см (у взрослых 4-12 мм)
maxDisInMinute = 2*maxbreathingRate*forwordExpantion*0.01; # максимаьное перемещение

v_max = maxDisInMinute/60; # максимальная скорость грудной клетки в м/с

fd_max = speed2dop(2*v_max,lambda); # максимальный Доплер
fb_max = fr_max + fd_max;

fs = max(2*fb_max, bw); # расчет частоты дискретизации

Следующими строками создаем зондирующий сигнал. Для этого применяем встроенный в Engee системный объект – EngeePhased.FMCWWaveform() и строим картинку сигнала (рисунок 8) сигнала.

# Формируем генератор ЛЧМ-сигнала
waveform = EngeePhased.FMCWWaveform(SweepTime=tm,SweepBandwidth=bw,SampleRate=fs)
sig = waveform()

# Построение спектрограммы
plotting_spectrogram(sig, fs, tm)

Определяем сценарий моделирования, где цель задается с помощью системного объекта EngeePhased.Target(), а ее движение – с помощью EngeePhased.Platform(). Среда распространения между антенной и грудной клеткой в виде свободного пространства задается через EngeePhased.FreeSpace().

# Задаем сценарий
sweeptime =waveform.SweepTime

senarioBreathingRate = 15; # количество дыхательных циклов в минуту
DisInMinute = 2*senarioBreathingRate*forwordExpantion*0.01;
chest_speed = DisInMinute/60; # скорость по сценарию 
chest_rcs = 1.4; # ЭПР (из справочника)
chest_dist = 0.5;# расстояние от излучателя до объекта

koeff = 60/senarioBreathingRate;
N_ot = (N_ot)*koeff/2; # количество отсчетов за полцикла

# Цель
chestTarget = EngeePhased.RadarTarget(MeanRCS=chest_rcs,PropagationSpeed=c,
    OperatingFrequency=fc);
 
# Движение цели
 chestmotion = EngeePhased.Platform(InitialPosition=[chest_dist;0;0],
    Velocity=[chest_speed;0;0])
 
 # Пространство
channel = EngeePhased.FreeSpace(PropagationSpeed=c,
    OperatingFrequency=fc,SampleRate=fs,TwoWayPropagation=true);

Определяем передатчик и приемник следующим образом: EngeePhased.Transmitter() и EngeePhased.ReceiverPreamp()

# Примерные параметры приемо-передатчика 
 ant_aperture = 6.06e-4; # в м^2
 ant_gain = aperture2gain1(ant_aperture,lambda); # в дБ
 
 tx_ppower = db2pow(5)*1e-3; # в Вт
 tx_gain = 9+ant_gain; # в Вт
 
 rx_gain = 15+ant_gain; # в Вт
 rx_nf = 4.5; # в Вт

# Задаем передатчик и приемник 
 transmitter = EngeePhased.Transmitter(PeakPower=tx_ppower,Gain=tx_gain)
 receiver = EngeePhased.ReceiverPreamp(Gain=rx_gain,NoiseFigure=rx_nf,
    SampleRate=fs,NoiseMethod="Noise power",NoisePower=0);

# Скорость приемо-передатчика
 radar_speed = 0;
 radarmotion = EngeePhased.Platform(InitialPosition=[0;0;0],Velocity=[radar_speed;0;0])

rng = MersenneTwister(2024)

Nsweep = 64;
xr = zeros(Complex,Int64(waveform.SampleRate*waveform.SweepTime),Nsweep); # инициализируем приемный сигнал 
l2r = 2; # коэффициент перевода из литров в расстояние
time_array = collect(V_Time).time
resp_array = l2r.*collect(V_Time).value;

Далее идет самое интересное – это процесс измерения:

# Цикл измерения
chest_speed = -chest_speed;
rangeEstimations = zeros(2,length(t)); # массив для оценок дальности
resp_array = l2r.*(collect(V_Time).value)*0.01;
counter = 0

for i = 1:length(t) 
    println("Шаг №$i")
    counter += 1      
    #  изменение направления дыхания
    if counter == ceil(Int64,N_ot)
        global chest_speed = -1 * chest_speed;
        global counter = 0;
    end   

    # обновление параметров цели
    release!(chestmotion)

    (i > 1) &&  (chestmotion.InitialPosition = [resp_array[i]+chest_dist, 0, 0])
    chestmotion.Velocity = [chest_speed;0;0];
    radar_pos, radar_vel = radarmotion(interval);
    global tgt_pos, tgt_vel = chestmotion(interval);
    for m = 1:1:Nsweep

        # Обновление положения радара и цели
        radar_pos, radar_vel = radarmotion(sweeptime);
        tgt_pos, tgt_vel = chestmotion(sweeptime);

        # Передача ЗС
        sig = waveform();
        txsig = transmitter(sig);
        
        # Канал распространения сигнала
        txsig = channel(txsig, radar_pos, tgt_pos, radar_vel, tgt_vel);  
        
        # Отражение сигнала от цели
        rxsig = chestTarget(txsig);
       
        # propagate the signal
        rxsig = rxchannel(rxsig,radar_pos,tgt_pos,radar_vel,tgt_vel);    
        
        # Прием и демодуляция отраженного сигнала
        rxsig = receiver(rxsig);
        xd = dechirp(rxsig,sig);
        xr[:,m] = xd;  
    end
    
    Dn = floor(Int64,fs/(2*fb_max));
    xr_d = zeros(ComplexF64,ceil(Int64,size(xr,1)/Dn),Nsweep);
    for m = size(xr,2):-1:1
        xr_d[:,m] = decimate(xr[:,m],Dn,"FIR");
    end
    fs_d = fs/Dn;

    # Расчет оценки частоты дыхания
    global fb_rng,_ = rootmusic(pulsint(xr_d,"coherent"),1,fs_d);

    # Перевол из частоты в дальность
    rng_est = beat2range(fb_rng,sweep_slope,c); 
    rangeEstimations[:,i].= [rng_est;i*interval];
end

Для оценки спектра целесообразно удалить постоянную составляющую:

samplingFreq = 1/interval;
TT = interval; 
L = length(t)      
k = (0:TT:L-1); 
rangeEstimations[1,:].= rangeEstimations[1,:].-mean(rangeEstimations[1,:]);

Строим спектр сигнала, результат на рисунке 9:

Y=fft(rangeEstimations[1,:]);
P2 = abs.(Y./L);
P1 = P2[1:Int64(L/2)+1];
P1[2:end-1] = 2*P1[2:end-1];

f = samplingFreq*(0:Int64(L/2))/L;
plot(f,P1,label="", line=(2,:blue,:sticks),xlims=(0,10)) 
scatter!(f,P1,label="",marker=(:circle,3,:yellow)) 
title!("Спектр отраженного сигнала (без патологии)")
xlabel!("Частота, Гц")
ylabel!("Модуль мощности, Вт")

Заключительный этап – нахождение максимальной частоты в спектре для оценки числа цикла вдоха-выдоха в минуту:

num_index = (f.==1).*Vector(1:length(f))
index = num_index[num_index.!=0][1]
peak = maximum(P1);
peakIndex = argmax(P1)

s2 = abs(peakIndex-index); 
breathingRate = f[peakIndex]*60/1;

println("Количество циклов дыхания в минуту = $(breathingRate)");

Оценка дыхания показала:

Как можно увидеть, разница есть. Изначально мы задавали его равным 15, да и на картинке 5 равно 15. Ошибка связана с точностью сетки по частоте при обработке сигнала, и при желании ее можно увеличить…

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

Заключение

Итак, мы рассмотрели две модели: модель легких, которая состоит из блоков среды Engee, и модель FMCW-радара, построенного на системных функция среды Engee. Соединили эти две модели и проверили корректность их работы. Как видите на данном примере, биорадиолокация дает очень наглядное и доступное представление о том, правильно ли функционируют легкие. Всем СПАСИБО за уделенное этой статье время. Надеюсь, материал был интересным и полезным!  

Готов ответить на ваши вопросы в комментариях!

Также хочу пригласить вас на наш семинар про моделирование радиолокационных систем с помощью российской среды Engee. Он состоится 27 ноября в Москве. Я и мои коллеги подготовили для вас много интересных практических примеров.

Регистрируйтесь по ссылке!