УДК 531, 53.072, 617-089

Россия, МГТУ им. Н.Э. Баумана

gouskov_am@mail.ru

kuzn-maria@yandex.ru

Введение

Применение роботизированных хирургических систем имеет много преимуществ. Роботизированные системы предлагают хирургам трехмерное изображение операционного поля и ощущение глубины, а также предоставляют в их распоряжение инструменты, которые в отличие от человеческих рук позволяют проводить длительные операции [1]. Если говорить о послеоперационных обострениях, то риски осложнений, связанных с искусственной вентиляцией легких, инфицированием раны и операционным травматическим стрессом, определяются фактором времени выполнения каждого этапа оперативного вмешательства [2]. Так как по сравнению с человеческой рукой инструменты роботизированных систем имеют управляемые с высокой точностью позиционирования степени свободы и меньшие размеры, то достигается бóльшая эффективность операций. Результаты операций менее зависят от квалификации проводящего их специалиста, тем самым снижается так называемый «человеческий фактор» при проведении длительных и требующих максимальной точности хирургических вмешательств.

В данной работе сделана попытка смоделировать процесс стернотомии (от др.-греч.στέρνον - грудина и τομή - разрез, рассечение) -хирургического маневра, заключающегося в рассечении грудины [2]. Стернотомия выполняется для обеспечения доступа к органам и патологическим образованиям переднего средостения: сердцу, отходящим от него крупным кровеносным сосудам и др. Задача численного моделирования стернотомии решена в программном пакете MATLAB^®  RoboticsToolbox + Simulink^®.

Нужно отметить, что при проведении робот-ассистированных операций управление роботом осуществляется исходя из усилий, действующих на инструмент. Данные усилия, возникающие при разрезании ткани, могут вызывать перемещения инструмента в направлениях, не совпадающих с направлением приложения рабочей силы. Это явление принято называть координатной связью [3].

Если инструмент в процессе операции отклоняется от прогнозируемой траектории, то в оперируемой ткани могут возникнуть серьезные нарушения. Поэтому оценка и учет координатной связи является одной из важнейших задач численного моделирования и проектирования робот-ассистированных операций [4]. В работе рассмотрена рабочая конфигурация робота, для которой произведена оценка перемещения конечной точки робота-манипулятора под действием нагрузки и вычислена матрица податливости. Решение получено методом конечных элементов в программном пакете ABAQUSCAE^®.

1.         Постановка задачи моделирования операции стернотомии

1.1.          Общие сведения

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

Рисунок 1.1 – а) Грудина; б) разрез с помощью электростернотома

Моделируемый в данной работе робот-ассистент должен провести операцию стернотомии в программном режиме с учетом координатной связи.

1.2.          Численное моделирование в MATLAB^® RoboticsToolbox

За простейшую основу моделируемого робота был взят прототип робота-манипулятора Puma 560 (рис. 1.2), рассмотренный в [3].

Рисунок1.2 – Робот-манипуляторPuma 560 (Westinghouse Limited)

После добавления третьего связующего звена к приведенному выше рисунку создана следующая кинематическая схема (рис. 1.3, а). Схема представляет собой систему из семи вращательных шарниров и трех связующих звеньев (). При моделировании численные значения длин связующих звеньев равны:  м,  м, м. Согласно примерам, рассмотренным в [3], был смоделирован робот-манипулятор (рис. 1.3, б). Листинг программы приведен в приложении 1.

а)                                                              б)

Рисунок 1.3 – а) Кинематическая схема робота; б) результат моделирования в MATLAB^® RoboticsToolbox

Далее рассматривается поведение робота при стернотомии. При этом жесткий Г-образный инструмент нагружается продольной силой  и моментом  (рис. 1.4)

Рисунок 1.4 – Упрощенная модель стернотома

Так как толщина кости по длине грудины меняется незначительно ([5, 6]), то усилие, действующее на стернотом, принимается постоянным и равным 2 Н.

Траектория движения задается следующим образом. В ходе операции стернотом проходит расстояние из точки А с координатами (0.2, 0, 0.55) до точки В с координатами (0.3, 0, 0.55) за 2 секунды. Т.е. инструмент совершает перемещение на расстояние 0.1 м вдоль оси X. После моделирования робота и задания движения можно вывести на экран траекторию движения стернотома в процессе операции (рис. 1.5).

Рисунок 1.5 – Траектория движения инструмента в проекции на плоскость

2.         Численное моделирование методом конечных элементов

Для оценки поведения точки приложения сил под действием нагрузки производится расчет перемещения инструмента и определение матрицы податливости для выбранной конфигурации робота.

На рис. 2.1 изображен общий вид рассматриваемого робота-ассистента.

1 – Зажим инструмента с установленным инструментом; 2 – блок приводов ориентации;    3 – предплечье; 4 – плечо; 5 –привод плеча; 6 – опорно-поворотное устройство
Рисунок 2.1 – Общий вид робота

Как было отмечено в пункте 1.2, робот имеет семь степеней свободы и состоит из шести звеньев. Первое звено соответствует установленному инструменту и смоделировано в данной работе упрощенно. Исследуемая конфигурация робота приведена на рис. 2.2.

Рисунок 2.2 – Рассматриваемая конфигурация робота

2.1.          Определение поля перемещений

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

Таблица 1 – Вес звеньев

При определении перемещений робота под действием только сил тяжести, силы, действующие на инструмент, не рассматривались.

В качестве типа конечного элемента выбран четырехузловой линейный элемент C3D4. Разбиение на конечные элементы представлено на рис. 2.3.

Рисунок 2.3 – Конечноэлементная модель робота. Количество элементов при данной разбивке равно 489110, число степеней свободы 2742846

Поле перемещений робота под действием сил тяжести всех звеньев представлено на рис. 2.4.

Рисунок 2.4 – Поле перемещений точек робота под действием сил тяжести

Компоненты вектора перемещений точки приложения рабочих сил инструмента по осям глобальной системы координат  представлены в таблице 2.

Таблица 2 –Перемещения на конце робота для рассматриваемой конфигурации

2.2.          Расчет матрицы податливости

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

где  - вектор перемещений,  - матрица податливости,  - вектор обобщенных сил. Матрица  симметричная и определенно положительная. Матричные вычисления проведены в глобальном декартовом базисе .

Вычисление матрицы податливости производится по следующей схеме:

1.     К выбранной конфигурации робота прикладываются по очереди три единичные силы в точке А в направлении соответствующих осей (рис. 2.5).

Рисунок 2.5 – Положение системы координат конечной точки инструмента

Точка А – точка приложения силы при расчете матрицы податливости, – оси локальной системы координат, расположенной в точке А.

2.     Вычисляются три компоненты вектора перемещений  для каждой силы.

3.     Заполняется соответствующая колонка матрицы податливости

Правильность вычислений доказывается симметричностью полученной матрицы податливости.

В соответствии с проведенными расчетами, ниже посчитана матрица податливости для рассматриваемой конфигураций робота.

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

Пусть единичный вектор сосредоточенной силы, приложенной к инструменту равен

Тогда согласно уравнению , единичный вектор смещения точки приложения силы оказывается равным

Относительное расположение векторов в пространстве представлено на рис. 2.6.

Рисунок 2.6 – Неколлинеарность векторов силы и перемещения

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

2.3.          Моделирование движения робота в SolidWorks®

Для моделирования операции стернотомии рассмотрена произвольная рабочая конфигурация робота, представленная на рис. 2.7, а. Траектория движения инструмента показана на рис. 2.7, б. В соответствии с параграфом 1.2, длина разреза выбрана равной 0,1 м.

а)                                                       б)

Рисунок 2.7 – а) рабочая конфигурация робота; б) траектория движения инструмента

Для определения матрицы податливости в процессе движения робота траектория движения инструмента в оперируемом объекте была разбита на 11 последовательных положений робота, включая начальную и конечную точки разреза. Каждое последующее положение взято на расстоянии 0,01 м от предыдущего. Данные положения представлены на рис. 2.8, где точка А является начальной точкой разреза, а точка В – конечной.

Рисунок 2.8 – Исследуемые положения инструмента

Алгоритм определения матриц податливости для каждого положения описан в параграфе 2.2. Матричные вычисления проведены в локальной системе координат  (см. рис. 2.5).

В соответствии с проведенными расчетами, ниже посчитаны матрицы податливости для 11 указанных положений робота. Индекс указывает на номер положения.

В пределах допустимой погрешности вычислений матрицы податливости являются симметричными.

Для указанных положений робота построен график отклонения реального перемещения (с учетом координатной связи) от заданного. Отклонение от заданной траектории представлено на рис. 2.9: красная линия и красные точки соответствуют программному положению инструмента (длина пути, равная 0.1 м, отображена в единицу), синие отрезки показывают отклонение конца инструмента от программного в одиннадцати равномерно расположенных узлах (за единицу принято отклонение в конечном положении, равное 4.4 мкм).

Рисунок 2.9 – Отклонение траектории инструмента от заданного движения для 11 положений робота

Заключение

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

Благодарность

Работа выполнена при поддержке Министерства образования и науки РФ в рамках федеральной целевой научно-технической программы «Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России» на 2007 – 2012 годы, номер контракта 16.523.11.3011.

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

1. Богданова Ю.В., Гуськов А.М., Жидких И.В., Нарайкин О.С., Шурхай В.А. Проблемы разработки стереотаксических систем наведения хирургического инструмента на примере биопсии головного мозга // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2012. № 12. DOI: 10.7463/0113.0486770
2. Информация о позвоночнике. Режим доступа: http://www.pozvonochnik.info/text/2/content/3/ (дата обращения 20.02.2013).
3. Corke P. Robotics, Vision and Control: Fundamental Algorithms in MATLAB. Springer, 2011. 558 c. DOI: 10.1007/978-3-642-20144-8
4. Банин Е.П., Барышева О.П., Батанов А.Ф., Богданова Ю.В., Гуськов А.М., Козубняк С.А., Нарайкин О.С., Пуценко А.А. Координатная связь медицинского робота-манипулятора // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2012. № 12.DOI: 10.7463/0113.0520630
5. Tillier Y., Paccini A., Durand-Reville M., Chenot J. Finite element modeling for soft tissue surgery based on linear and nonlinear elasticity behavior // Computational Aided Surgery. 2006. Vol. 11, no. 2. P. 63-68.
6. Shields T.W., Ponn R.B., Rusch V.W. General Thoracic Surgery. Lippincott Williams & Wilkins, 2004. 3040 p.
7. Bathe K.J. Finite Element Procedures. Prentice Hall, Englewood Cliffs, New Jersey, 1996. 1037 p.

Приложение 1

Программа численного моделирования условий работы хирургических инструментов с биотканями. Программный комплекс MATLAB Toolbox.

Содержание программы:

% Numerical Simulation of the three-link robotic arm for sternotomy incision

startup_rvc;         %function from matlab sourse

L1 = .2;             % length of link 1, [m]

L2 = .6;             % length of link 2, [m]

L3 = .4;             % length of link 3, [m]

tool_len = .2;       % tool length, [m]

z = .55;             % depth of laparotomy incision

x_startOut = .2;     % start coordinate for laparotomy incision, [m]

x_finalOut = .3;    % end coordinate for laparotomy incision, [m]

% length of the incision = 5 sm

y = 0;               % laparotomy incision will be in (XZ) plane

t_stepOut = .05;     % time step for outside motion, [s]

t_finalOut = 2;

s = 'Rz(q1) Tz(L1) Ry(q2) Rz(q3) Tz(L2) Ry(q4) Rz(q5) Tz(L3) Rz(q6) Ry(q7)'; %10 items

dh = DHFactor(s);

dh.display;

cmd = dh.command('three-link robot arm');

robot = eval(cmd);

q0 = [0 0 0 0 0 0 0];                       % initial conditions for 7 items

robot.tool = transl(tool_len,0,0);          % tool creating *rpy2tr(0,pi/2,0);

figure('Name','three-link robot arm');      % robot plot

robot.plot(q0);

robot.teach();

T1_out = transl(x_startOut,y,z);            % robot trajectory

T2_out = transl(x_finalOut,y,z);

q1 = robot.ikine(T1_out);

q2 = robot.ikine(T2_out);

t_out = (0:t_stepOut:t_finalOut);

Ts = ctraj(T1_out,T2_out,length(t_out));

qi = robot.ikine(Ts);                      

T = robot.fkine(qi);

p = transl(T);

figure('Name','Displacement of the cutting tool in ZX plane');

holdon; grid on; box on;

xlabel('X');

ylabel('Z');

plot(p(:,1), p(:,3),'LineWidth',2);

plot(p(1,1), p(1,3),'ko',p(end,1), p(1,3),'*k','LineWidth',2);

set(gca,'ylim',[0.54 0.56]);               %to avoid unstability

figure('Name','Joint Coordinates');

holdon; grid on; box on;

xlabel('time');

ylabel('q');

plot(t_out,qi,'LineWidth',2);

legend('q1','q2','q3','q4','q5','q6','q7');

x_start = x_finalOut;                       % start/final coordinates for the incision

x_final = 0.35;

dist = abs(x_final-x_start);                % distance for tool motion

x_step = .01;                               % step along trajectory

step_num = dist/x_step;                     % number of steps

% start/final/step time for inside motion

t_start = t_finalOut;

t_final = t_finalOut*2;

t_step = (t_final-t_start)/step_num;

% forces

Fx = 2;

Fy = 0;

Mx = 0;

My = Fx*tool_len;

Mz = 0;

T1 = transl(x_start,y,z);

T2 = transl(x_final,y,z);

q1 = robot.ikine(T1);

q2 = robot.ikine(T2);

t = (t_start:t_step:t_final)';              % time step, [s]

Ts = ctraj(T1,T2,length(t));

qqi = robot.ikine(Ts);

TT = robot.fkine(qqi);

pp = transl(TT);

% assign the load

Q = [Fx,Fy,Fz,Mx,My,Mz];

Torques = robot.jacob0(q2)'*Q';

% conversion qqi (cell-type) to matrix for animation

qi = [qi; qqi];

% animation

figure('Name','incision of 5 sm length');

robot.plot(qi);