NGHIÊN CỨU PHÁT TRIỂN KỸ THUẬT CT THẾ HỆ THỨ IV ỨNG DỤNG TRONG CÔNG NGHIỆP DẦU KHÍ Việt Nam

Hội nghị khoa học Sau đại học trường Đại học Đà lạt, năm 2014. 1 Hình 2: Các thế hệ của thiết bị CT NGHIÊN CỨU PHÁT TRIỂN KỸ THUẬT CT THẾ HỆ THỨ IV ỨNG DỤNG TRONG CÔNG NGHIỆP DẦU KHÍ VIỆT NAM Phạm Văn Đạo, Đặng Nguyễn Thế Duy, Mai Công Thành, Nguyễn Văn Chuẩn, Bùi Trọng Duy Trung tâm Ứng dụng Kỹ thuật Hạt nhân trong Công nghiệp, Địa chỉ: số 01 đường DT723, P.12, Đà Lạt, Lâm Đồng. E-mail: office@canti.vn; daopv@canti.vn Tóm tắt Cùng với sự phát triển của khoa học và kỹ thuật, việc sử dụng các thiết bị để khảo sát các đối tượng trong nhà máy dầu khí đang được yêu cầu về độ chính xác cao cũng như hiệu quả kinh tế mà nó mang lại. Để đáp ứng nhu cầu kiểm tra, chẩn đoán tình trạng các đối tượng đường ống trong công nghiệp dầu khí ở Việt Nam, nhóm nghiên cứu của Trung tâm Ứng dụng kỹ thuật hạt nhân trong công nghiệp đã tiến hành nghiên cứu, phát triển thiết bị soi cắt lớp thế hệ thứ IV. Đề tài này cũng là kết quả kế thừa và phát triển một cách liên tục qua những đề tài nghiên cứu trước đó tại Trung tâm từ năm 2007 đến nay. Trong bài báo này trình bày các vấn đề về cấu hình thiết bị, thuật toán Tối đa hóa kỳ vọng (EM) và thuật toán Chiếu ngược có lọc (FBP) cho tái tạo ảnh CT thế hệ thứ 4 và một số kết quả khảo sát nhằm đánh giá khả năng vận hành của thiết bị đã được chế tạo. Từ khóa CAT, EM, FBP, chụp cắt lớp điện toán nhanh, tái tạo hình ảnh, soi gamma truyền qua,... I. GIỚI THIỆU Năm 2011, Trung tâm Ứng dụng kỹ thuật hạt nhân trong công nghiệp (CANTI) đã chế tạo thành công thiết bị chụp cắt lớp điện toán đầu tiên theo cấu hình 1 nguồn 1 đầu dò – CT thế hệ thứ I, mang tên GORBIT, đã được ứng dụng rất nhiều trong việc khảo sát khuyết tật các đường ống trong công nghiệp và được xuất khẩu sang các nước trong khu vực theo đơn đặt hàng của IAEA. GORBIT là thiết bị chụp cắt lớp công nghiệp, cho phép chụp ảnh cắt lớp các thiết bị có đường kính tối đa 60cm, Hình 1: Thiết bị chụp cắt lớp GORBIT Hội nghị khoa học Sau đại học trường Đại học Đà lạt, năm 2014. 2 thiết bị này đã được ứng dụng để khảo sát khuyết tật trên các đường ống dẫn khí và một số thiết bị khác trong công nghiệp. Hiện nay trên thế giới, thiết bị chụp cắt lớp điện toán đã phát triển qua nhiều thế hệ, ngày càng nhỏ gọn trong kết cấu cơ khí và rút ngắn thời gian đo. Trong đó thế hệ chụp cắt lớp điện toán thứ IV có cấu hình gồm nguồn và các detector có quỹ đạo nằm trên đường tròn đồng tâm như hình 2 [1], [2], [3]. Ưu điểm của kỹ thuật soi cắt lớp thế hệ thứ IV là có thể rút ngắn thời gian khảo sát và cấu hình được bố trí đơn giản [4]. Nhóm nghiên cứu đã xây dựng và tính toán hình học cấu hình, đồng thời phát triển chương trình dựng ảnh cắt lớp bằng hai thuật toán Tối đa hóa kỳ vọng (EM) và thuật toán Chiếu ngược có lọc (FBP). Sau đó cấu hình vật lý được xây dựng để khảo sát và đánh giá khả năng ứng dụng của phương pháp trong các điều kiện thực tế. II. THUẬT TOÁN TÁI TẠO ẢNH CT Trong cấu hình CT thế hệ thứ 4, nguồn và detector được bố trí theo hình học chùm quạt. Chất lượng hình ảnh tái tạo phụ thuộc vào mật độ phép đo trên một lát cắt qua vật. Tuy nhiên trên cùng số phép đo với số lượng các phép tương đối ít khoảng dưới 64 x 16 (số hình chiếu x số tia), chất lượng hình ảnh được tái tạo bằng thuật toán Tối đa hóa kỳ vọng cho kết quả tương đối tốt hơn thuật toán Chiếu ngược có lọc. Hình 3: Hình học trong CT cấu hình song song – CT thế hệ thứ I (trái) và cấu hình chùm quạt – CT thế hệ thứ III, IV, V (phải)   cos.sin. sin.cos. sty stx         cos.cos. sin.sin. DLy DLx   Hội nghị khoa học Sau đại học trường Đại học Đà lạt, năm 2014. 3 Thuật toán Chiếu ngược có lọc [1], [4], [5] được thực hiện qua 3 bước biến đổi:  Tính biến đổi Fourier của bộ dữ liệu hình chiếu p(t,θ) → P(ξ,θ);  Biến đổi ngược Fourier với bộ lọc cao qua H(ξ);  Chiếu ngược bộ dữ liệu g(t,θ) đã qua biến đổi để thu được hình ảnh μ(x,y). Trong thuật toán Chiếu ngược có lọc, tia tổng được định nghĩa như biểu thức:    tI tI d tp 0ln1),(  (1) Ở đây, I0(t), I(t) là số đếm đo được trong không khí và khi truyền qua vật, d tương đương với khoảng cách giữa nguồn và đầu dò. Ký hiệu P(ξ,θ) là phép Biến đổi Fourier của p(t,θ), biểu thức toán học được viết lại như sau:        dttitpP  2exp),(, (2) Ở đây, g(t,θ) được ký hiệu là phép Biến đổi ngược Fourier của tích chập H(ξ)* P(ξ,θ) theo công thức (3), với H(ξ) là hàm lọc.             dtiPHtg 2exp,, (3) Hình ảnh được tái tạo μ(x, y) thu được từ phép Chiếu ngược của g(t,θ) theo công thức (4):      0 ,sincos),( dyxgyx (4) Hội nghị khoa học Sau đại học trường Đại học Đà lạt, năm 2014. 4 Thuật toán Tối đa hóa kỳ vọng [6], [7] được thực hiện qua 2 bước biến đổi sau:  Bước tính kỳ vọng dựa vào sự suy giảm của cường độ tia gamma theo bề dày hấp thụ của vật liệu (bước E); Giả sử I0j, Ij là số photon ghi được tại đầu dò khi không có vật và có vật, ljk là phần chiều dài đóng góp của pixel thứ k trên hình chiếu thứ j. Sự phát xạ và sự ghi nhận bức xạ là quá trình rời rạc theo thời gian. Số đếm trung bình thu được là: 0 1 .exp jq j jk jk k I l           (5) Hình 4: Lưu đồ thuật toán Chiếu ngược có lọc (FBP) Kết thúc Chiếu ngược g(t,θ) và đưa vào hình ảnh f(x,y) zero padding p(t,θ) → p’(t,θ) Nhân với bộ lọc P(ω,θ)*H(ω)→ G(ω,θ) Biến đổi ngược Fourier G(ω,θ)→ g(t,θ) Hình chiếu cuối cùng Yes No Biến đổi Fourier p’(t,θ) → P(ω,θ) Dữ liệu hình chiếu p(t,θ) tại các góc θ Bắt đầu Chọn hình chiếu mới Hội nghị khoa học Sau đại học trường Đại học Đà lạt, năm 2014. 5 Trong thực tế số photon thu được tại detector là Ij, phân bố Poisson quanh giá trị trung bình trong phương trình (5). Do đó, hàm mật độ xác suất (còn gọi là hàm hợp lý - likelihood) mà các photon Ij thu được tại detector từ hình chiếu thứ j là:   0 1 0 1 .exp exp .exp ! j j j Iq j jk jk q k j j jk jk kj I l p I I l I                            (6) Bởi vì toàn bộ dữ liệu đo được được tạo nên từ các hình chiếu riêng độc lập với nhau, hàm hợp lý của phép đo toàn bộ dữ liệu đơn giản là tích mật độ xác suất của mỗi phép đo hình chiếu riêng lẻ j. Do đó hàm hợp lý của toàn bộ dữ liệu đo được là:    ,j jk j j g I p I  (7)   0 0 1 1 ln , .exp ln ln ! j jq q j jk j jk jk j jk jk j j j j k k g I I l I l I I I                        (8) Để thu được bước E của thuật toán EM, ta cần phải ước đoán dữ liệu đầy đủ. Dữ liệu đầy đủ của trường hợp chụp ảnh cắt lớp là số photon đi vào mỗi pixel dọc theo mỗi hình chiếu. Xét một hình chiếu đơn giản (thứ j) với các pixel được đánh số từ 1 tới (qj -1) từ nguồn đến detector, 1, 1jk q    (hình 5), với pixel thứ qj tương ứng với phép đo của detector: Ij Mấu chốt để viết chính xác hàm xác suất của dữ liệu đầy đủ là cần chú ý đến số photon rời khỏi pixel (phụ thuộc số photon đi vào pixel đó và xác suất Hình 5: Mô hình xác định các giá trị kỳ vọng trên một tia chiếu Hội nghị khoa học Sau đại học trường Đại học Đà lạt, năm 2014. 6 truyền qua). Do đó, xác suất Xj,k + 1 photon rời khỏi pixel k được cho bởi Xj,k photon đi vào pixel đó, tuân theo hàm phân phối nhị thức mà chỉ có hai hệ quả có thể xảy ra: truyền qua và hấp thụ. Hàm mật độ nhị thức tương ứng được cho bởi:        , 1 , , 1,, 1 , , , , , 1 exp 1 exp j k j k j kX X Xj k j k j k j k j k j k j k X p X l l X                        (9) Bởi vì quá trình phát xạ từ nguồn tuân theo phân bố Poisson với giá trị trung bình là 0 jI , nên mật độ xác suất của các photon được phát ra từ pixel đầu tiên là:     ,1 0 0 ,1 ,1 .exp ! jX j j j j I I p X X   (10) Vì mỗi pixel là độc lập với các pixel khác, nên hàm xác suất cho tập hợp toàn bộ dữ liệu đầy đủ chỉ đơn thuần là tích các hàm xác suất riêng lẻ của mỗi pixel dọc theo hình chiếu j :      ,1 , 1 1 , jq j j j k k f X p X p X     (11) Lấy log tự nhiên (ln) hai vế của phương trình (11) ta được:         0 ,1 0 ,1 , , 1 , , , , 1 , , 1 , 1 ln , ln ln ln ln exp ln 1 exp j j j j j j q j k j k j k j k j k j k j k j k k j k f X I X I X X X l X X l X                                  (12) Đối với bước E trong thuật toán EM, chúng ta cần tính toán giá trị kỳ vọng có điều kiện  , | ,j k jE X I  của tập dữ liệu đầy đủ Xj,k dựa trên tập dữ liệu không đầy đủ jI (dữ liệu đo được) và ước đoán giá trị hiện thời của tham số vector μ. Tương ứng với phương trình (5), giá trị kỳ vọng của các photon được truyền ra ngoài pixel thứ 1k là:   1 , , 0 , , 1 exp k j k j k j j t j t t E X I l             (13) Hội nghị khoa học Sau đại học trường Đại học Đà lạt, năm 2014. 7 Tương tự, trung bình phân bố của Ij dựa trên điều kiện Xj,k (tức là giá trị kỳ vọng quan sát Ij được cho bởi Xj,k ) là:   1 0 , , 1 .exp jq j j j t j t t E I I l             (14) Bước ước đoán tính toán kỳ vọng mô hình thống kê của dữ liệu đo Ij (từ đầu dò). Theo K. Lange và R. Carson, biểu thức được viêt như sau [6]:      jkjjjkj IEXEIIXE  ,, | (15) Phương trình (15) là kỳ vọng có điều kiện cần thiết cho thuật toán EM. Tập hợp trên tất cả các hình chiếu, biểu thức đại diện cho kỳ vọng của tập dữ liệu đầy đủ được cho trên tập dữ liệu đo.  Bước tính tối đa hóa cho hàm kỳ vọng (bước M). Trước khi cực đại tham số μ, sử dụng kết quả phương trình (15) để tính toán kỳ vọng của các photon đi vào và rời khỏi pixel k, Nj,k và Mj,k tương ứng. Các giá trị kỳ vọng sau đó được thay vào phương trình (12) và được lấy tổng trên tất cả các hình chiếu để thu được bộ dữ liệu đầy đủ:       , , , , , , , 1 ln exp ln 1 exp jq j k j k j k j k j k j k j k j k M l N M l R               (16) Quá trình cực đại hóa được thực hiện bằng cách tính đạo hàm riêng phần theo μ và cho đạo hàm bằng 0 với mỗi k.     , , , , , , , . 0 exp 1 j k j k j k j k j k j Jk j Jk j k j k l M l N M l          (17) Phương trình (15), (16) và (17) được sử dụng một cách lặp lại lần lượt tạo thành toàn bộ thuật toán EM, điều này dẫn đến sự hội tụ nhanh chóng của thuật toán. Phương trình (17) không thể tìm μk một cách chính xác, ta sử dụng phương pháp xấp xỉ dưới đây:     31 1 1 exp 1 2 12 s O s s s      (18) Hội nghị khoa học Sau đại học trường Đại học Đà lạt, năm 2014. 8 với s = lj,k μj,k (phương trình (18) xuất phát từ khai triển s/(es – 1) trong chuỗi Taylor). Thay thế xấp xỉ này vào phương trình (17) ta có 3 trường hợp [6]: - Trường hợp thứ 1: ta giữ lại 1 số hạng 1 s trong công thức (18) ta được:  , , 1 , , ,. j k j k j Jkn j k j k j k j Jk N M M l        (19) - Trường hợp 2: để lời giải xấp xỉ tốt hơn, ta giữ lại 2 số hạng đầu trong công thức (18) ta được:     , , 1 , , , , 1 . 2 j k j k j Jkn j k j k j k j k j Jk N M M N l         (20) Trong phương trình (19) và (20) thực tế được sắp xếp là ngưỡng trên và ngưỡng dưới tương ứng với nghiệm đúng 1, n j k  theo bất phương trình dưới đây: 1 1 1 1 , 0 2 1s s s e s      (21) - Trường hợp 3: để lời giải xấp xỉ tốt nhất, ta sử dụng cả 3 số hạng trong công thức (18) ta được:       2 , ,2 , , , , , , , ,. . 012 2 j k j k j k j k j k j k j k j k j k j k j Jk j Jk j Jk l l N M N M N M             (22) Với:      jkjjjkjkj IEXEIIXEN  ,,, | (23)      jkjjjkjkj IEXEIIXEM   1,1,, | (24) Từ phương trình (22) ta đặt: 2 , ,. . 0j k j kA B C    (25) Nghiệm của phương trình bậc hai:  21 , 4 2 n j k B B AC A       (26) Hội nghị khoa học Sau đại học trường Đại học Đà lạt, năm 2014. 9 Hình 6: Lưu đồ thuật toán Tối đa hóa kỳ vọng (EM) Maximization )|( )|( 1,, ,, txyxy txyxy IXEM IXEN                 viewj xyxy viewj xy xyxy viewj xy xyxy NMC l NMB l NMA ,, 2 , ,, 2 , ,, 2 12 A ACBBk xy 2 42 ,   - Bộ số liệu đo It , I0 - Hình ảnh ước đoán ban đầu µy,x0 - Điều kiện hội tụ eps    allpixels k xy k xy eps, 1 ,  k xyxy ,,   No Yes Hình ảnh tái tạo           1 1 ,,0, exp.)( x i iyiyxy lIXE            1 ,,, exp).()( n xi iyiyxyt lXEIE  )()()|( ,, txyttxy IEXEIIXE  Expectation Hội nghị khoa học Sau đại học trường Đại học Đà lạt, năm 2014. 10 III. KẾT QUẢ NGHIÊN CỨU 1. Cấu hình thiết bị Với mục tiêu ban đầu xây dựng cấu hình nghiên cứu, khảo sát các thông số, đặc trưng của cấu hình thiết bị để tiến đến chế tạo một thiết bị chụp cắt lớp chuyên dụng cho đối tượng đường ống công nghiệp dầu khí. Dựa trên kinh phí của đề tài, cấu hình ban đầu được lắp đặt như sau: - Nguồn phóng xạ: Cs – 137 hoạt độ 20 – 50 mCi; Se – 75 hoạt độ 1,5 Ci; - Dectector NaI(Tl) kích thước ½ x 1 inch (12 detector); - Kích thước hình ảnh tối đa: 512 x 512 pixel; - Đường kính vật thể tối đa: 600 mm; - Dải mật độ: 0 – 7.8g/cm3; - Thời gian trung bình/lát cắt :1,5 giờ. Phần mềm tái tạo ảnh: Phần mềm tái tạo hình ảnh CT cấu hình thế hệ thứ IV được phát triển trên ngôn ngữ lập trình Visual Basic 6.0 với các thuật toán tái tạo ảnh: Tối đa hóa kỳ vọng (EM), Chiếu ngược có lọc (FBP) và các thuật toán xử lý hình ảnh khác. 2. Kết quả thực nghiệm a. Thí nghiệm 1 Thử nghiệm chụp ảnh cắt lớp trên phatom có hình dạng, kích thước và vật liệu như hình 8. Số phép đo 128 x 240 (Số hình chiếu x số tia). Hình 7: Thiết bị thử nghiệm CT cấu hình chùm quạt Hình 8: Phantom thử nghiệm chụp CT Hội nghị khoa học Sau đại học trường Đại học Đà lạt, năm 2014. 11 Kết quả Từ bộ dữ liệu đo được của cấu hình chùm quạt trên, hình ảnh được tái tạo bằng thuật toán EM và thuật toán FBP như hình 9. b. Thí nghiệm 2 Thí nghiệm được thực hiện trên mẫu đường ống bố trí như hình 10. Với đường kính ngoài của đường ống 275mm, kích thước khối ghỗ 80 x 90mm, đường kính của khối teflon 65mm, bề dày của khối paraffin 25mm, thanh lục giác có bán kính 5mm và khuyết tật bên ngoài đường ống lõm 4mm. Số phép đo 64 x 512 (Số hình chiếu x số tia). Kết quả Hình ảnh được tái tạo bằng thuật toán Chiếu ngược có lọc của phantom như hình 11 (b), (c) với các thang màu Gray và Jet. Hình 9: Phantom (a), hình ảnh tái tạo bằng thuật toán EM (b) và thuật toán FBP (c) (a) (b) (c) Hình 11: Phantom (a), hình ảnh tái tạo thuật toán FBP (b) (c) (a) (b) (c) (c) Hình 10: Bản vẽ và kích thước của phantom (c) Hội nghị khoa học Sau đại học trường Đại học Đà lạt, năm 2014. 12 IV. KẾT LUẬN Từ các kết quả tái tạo và phân tích hình ảnh, nhóm nghiên cứu rút ra một số kết luận như sau: - Việc ứng dụng kỹ thuật chụp cắt lớp điện toán với cấu hình thế hệ thứ IV cho kết quả hình ảnh là có thể chấp nhận được với độ phân giải hiện tại 4 mm; - Tuy độ phân giải hình ảnh không tốt bằng cấu hình CT thế hệ thứ I trước đó (1mm) [8], nhưng với cấu hình mới này thời gian chụp cắt lớp được rút ngắn hơn 8 lần. Qua đó nhóm nghiên cứu cũng đề ra một số hướng tiếp tục phát triển kỹ thuật này và hoàn thiện thiết bị trong tương lai để có thể triển khai dịch vụ ngoài hiện trường: - Nâng cao độ phân giải về vật liệu cũng như độ phân giải không gian bằng việc sử dụng nguồn đa năng lượng và nghiên cứu thiết kế cơ cấu chuẩn trực cho các detector. - Tiếp tục xây dựng và phát triển phần mềm tái tạo hình ảnh, hỗ trợ thêm các thuật toán xử lý ảnh nhằm nâng cao chất lượng hình ảnh tái tạo. Hội nghị khoa học Sau đại học trường Đại học Đà lạt, năm 2014. 13 TÀI LIỆU THAM KHẢO [1] IAEA, Industrial Process Gamma Tomography, Final report of a coordinated research project 2003–2007. [2] Jongbum KIM, Monte Carlo Simulation for the Design of Industrial Gamma – ray Transmission Tomography, Korea Atomic Energy Research Institute, 2011. [3] Jiang Hsieh. Computed Tomography principles, design, artifacts, and recent advances. 2nd Ed, 2009. [4] Emssion Tomography – The Fumdamental of PET and SPECT, Miles N. Wernick, Elsevier, 2004. [5] Edwin L.Dove, Notes on Computerized Tomography – imaging Fundamental, 2003. [6] Kenneth Lange and Richard Carson, EM Reconstruction Algorithms for Emission and Transmission Tomography, 1984. [7] Zeljko V.Kuzeljevic & prof. Muthanna H. Al–Dahhan. Expectation – Maximization (EM) algorithm and its use for CT Imaging. [8] Báo cáo đề tài CS/08/06 – 01: Nghiên cứu chế tạo thử nghiệm thiết bị CT công nghiệp loại một nguồn – một detector quy mô phòng thí nghiệm. Hội nghị khoa học Sau đại học trường Đại học Đà lạt, năm 2014. 14 RESEARCH AND DEVELOP THE FOURTH GENERATION CT TECHNIQUE APPLY FOR PETROLEUM INDUSTRY IN VIETNAM Phạm Văn Đạo, Đặng Nguyễn Thế Duy, Mai Công Thành, Nguyễn Văn Chuẩn, Bùi Trọng Duy Center for Applications of Nuclear Technique in Industry, Address: 01 - DT723 street, 12 ward, Dalat city, Lamdong. E-mail: office@canti.vn; daopv@canti.vn Abstract: Along with the development of science and technology, the use of equipment to examine objects in the oil and gas plants is being asked about high accuracy as well as economic efficiency which it brings. In order to meet the needs test, diagnose the condition of pipelines in the petroleum industry in Vietnam, the research group of the Center for application of nuclear techniques in industry (CANTI) have studied and developed the 4th generation CT device. This project is also result from inheriting and continuously developed through previous research projects at CANTI from 2007 to now. In this paper presents the issues of device configuration, Expectation Maximization algorithm (EM) and Filter Back Projection algorithm (FBP) for 4th generation CT image reconstruction and some of the surveyed results in order to evaluate the operation of the equipment has been manufactured. Key words CAT, EM, FBP, fast CT, CT reconstruction, Transmission Tomography,...