Биорадиолокация в Engee

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

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

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

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

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

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

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

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

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

Поэтому хочу с вами поделиться относительно простой, но в то же время полезной моделью биорадара. Модель разбита на две части. Первая часть посвящена моделированию перемещения грудной клетки человека. Вторая часть модели – про разработку 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 ноября в Москве. Я и мои коллеги подготовили для вас много интересных практических примеров.

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