Ứng dụng biến đổi Hilbert-Huang để chẩn đoán hư hỏng trong kết cấu dưới của cầu

Nội dung text: Ứng dụng biến đổi Hilbert-Huang để chẩn đoán hư hỏng trong kết cấu dưới của cầu

1. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG ỨNG DỤNG BIẾN ĐỔI HILBERT-HUANG ĐỂ CHẨN ĐOÁN HƯ HỎNG TRONG KẾT CẤU DƯỚI CỦA CẦU PGS. TS. BÙI ĐỨC CHÍNH Trường Đại học Giao thông vận tải Tóm tắt: Bài báo này giới thiệu những kết quả nghiên cứu ban đầu về chẩn đoán hư hỏng của kết cấu dưới công trình cầu dựa trên biến đổi Hilbert-Huang. Bài báo bao gồm: Nội dung của biến đổi Hilbert; việc phân tích một tín hiệu thành các hàm dạng/mode bản chất theo phân tích dạng kinh nghiệm và phân tích dạng kinh nghiệm quần thể; một số kết quả ban đầu trong áp dụng phổ giới hạn Hilbert trong chẩn đoán hư hỏng của kết cấu dưới của công trình cầu. Một số từ viết tắt FFT Fast Fourier Transform (Biến đổi Fourier nhanh). WT Wavelet Transform (Biến đổi Wavelet). HT Hilbert Transform (Biến đổi Hilbert). HHT Hilbert-Huang Transform (Biến đổi Hilbert-Huang). IMS Intrinsic Mode Function (Hàm dạng/mode bản chất). EMD Empirical Mode Decomposition (Phân tích dạng kinh nghiệm). EEMD Ensemble Empirical Mode Decomposition (Phân tích dạng kinh nghiệm quần thể) NHS Nominal Hilbert Spectrum (Phổ Hilbert danh định) MHS Marginal Hilbert Spectrum (Phổ Hilbert giới hạn) 1. Đặt vấn đề Biến đổi Hilbert thích hợp để xử lý các tín hiệu không dừng và giải hẹp. Biến đổi Hilbert được Khi trong công trình/kết cấu xuất hiện các hư định nghĩa như sau [3] : hỏng khuyết tật sẽ dẫn tới thay đổi các đặc trưng động học như: Giảm độ cứng, thay đổi tần số dao 1 x  y t = H x t = P d (1) động tự do; thay đổi dạng/mode dao động của π - t - công trình/kết cấu Kỹ thuật chẩn đoán phát hiện trong đó: hư hỏng khuyết tật dựa vào sự thay đổi các đặc Toán tử H[.] là biến đổi Hilbert; trưng động học đang được nghiên cứu áp dụng. Tuy nhiên trong ứng dụng thực tế, kỹ thuật này P là giá trị chính Cauchy. cũng gặp khá nhiều hạn chế như: Việc xác định Từ (1) bất cứ một tín hiệu z(t) đều có thể biểu chính xác các thông số dao động như các tần số diễn thành tổng của phần thực x(t) và phần ảo riêng và các dạng dao động riêng trên kết cấu y(t) của nó thực khá khó khăn, các ảnh hưởng của các đại iθ t lượng cần đo lớn Một trong các vấn đề dẫn tới zt =xt+iyt =ate (2) các hạn chế trên đó là quá trình xử lý các tín hiệu trong đó: a(t) là biên độ tức thời và (t) là pha tức dao động ghi nhận được ở hiện trường trên công thời, chúng được tính như sau: 2 2 trình/kết cấu thực còn gặp nhiều khó khăn [1, 4]. a t = x t + y t Sau đây xin giới thiệu các nét cơ bản về y t (3) θ t = artan EMD, EEMD, NHS, MHS trong biến đổi Hilbert- x t Huang và áp dụng chúng trong phân tích các dữ Tần số tức thời trong phép biến đổi Hilbert liệu dao động thu được để chẩn đoán hư hỏng được tính như sau : của kết cấu dưới của công trình cầu. dθt ytxt-ytxt  ω t = 2πf t = = (4) 2 2 2. Biến đổi Hilbert-Huang dt x t + y t 2.1 Biến đổi Hilbert Như vậy, phần thực của tín hiệu x(t) có thể được Tạp chí KHCN Xây dựng – số 1/2016 25
2. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG biểu diễn qua biên độ và tần số tức thời như một của chuỗi thời gian. Giá trị trung bình của các hàm phụ thuộc thời gian : đường bao trên và đường bao dưới, m1(t) được trừ đi từ tín hiệu nguyên thủy để nhận được i ω t dt xt=Rzt =Rate (5) thành phần thứ nhất h1(t) của quá trình sàng lọc này: trong đó: R(.) ký hiệu phần thực của tín hiệu phân h t = x t - m t tích z(t). 1 1 (6) Để đảm bảo chắc chắn rằng tần số tức thời Nếu h1(t) là một IMF, quá trình sàng lọc sẽ nhận được từ (4) có ý nghĩa vật lý, pha tức thời được dừng lại. Hai điều kiện để kiểm tra xem (t) cần phải là hàm đơn trị tại bất kỳ giá trị nào h1(t) có phải là một IMF là: (i) Số điểm về không của t. cần phải bằng số của cực trị hoặc không khác quá số cực trị là 1; (ii) h1(t) có tính đối xứng giữa Như vậy phép biến đổi Hilbert có thể được sử đường bao trên và đường bao dưới qua trục dụng để nghiên cứu một chuỗi tín hiệu dưới dạng zero. hàm suy rộng thời gian - tần số. Nhưng thật không may mắn, phạm vi áp dụng phép biến đổi Nói cách khác, quá trình sàng lọc cần được Hilbert thường rất ngặt nghèo, nó đòi hỏi đặc tính lặp lại để lọc tín hiệu h1(t) thành một IMF. Tương của x(t) phải là giải hẹp theo thời gian t. Điều kiện tự, h1(t) được sàng lọc để nhận được thành phần này thường không thỏa mãn khi trong thực tế kỹ sàng lọc thứ nhất h11(t). thuật, các chuỗi tín hiệu - thời gian thường không h t = x t - m t 11 11 (7) dừng và giải rộng. Thí dụ giả thiết rằng, có tín trong đó m11(t) là giá trị trung bình của đường hiệu x(t) = cos(1t)+sin(2t), phép biến đổi Hilbert bao trên và đường bao dưới của h (t). Quá trình sẽ tạo ra một tần số tức thời trung bình thay cho 1 này được tiếp tục cho đến khi h1k(t) là một IMF. các tần số 1 và 2 và bản chất tín hiệu đã bị Hàm h (t) được ký hiệu là thành phần thứ nhất thay đổi. 1k c1(t) = h1k(t). Người ta thường sử dụng độ lệch 2.2 EMD và IMF chuẩn để làm tiêu chuẩn kiểm tra quá trình sàng Để khắc phục vấn đề này, Huang và đồng lọc: nghiệp [3] đã đề nghị phương pháp phân tích 2 h t - h t kinh nghiệm EMD để tách các IMF’s từ một tín n 1,k-1 1,k SD = < ε (8) 2 hiệu theo thời gian, mà mỗi IMF chứa chỉ một i=1 h t 1,k-1 dạng dao động đơn giản (một tín hiệu dải hẹp tại trong đó  thường lấy từ 0,2 đến 0,3 [3]. Để thời điểm đang xét). đẩy nhanh quá trình sàng lọc, người ta cũng Một thuật toán EMD đã được đề nghị để tạo thường sử dụng tiêu chuẩn kiểm tra: ra các IMF’s một cách đơn giản, được gọi là quá n 2 trình sàng lọc (Sifting Process). Có 3 giả thiết đối h t - h t  1,k-1 1,k với EMD: (i) Tín hiệu cần có ít nhất hai cực trị - SD =i=1 < ε (9) n 2 một cực tiểu và một cực đại; (ii) khoảng thời gian  h t i=1 1,k-1 giữa các cực trị (time scale) phải xác định được Tiêu chuẩn kiểm tra quá trình sàng lọc được đặc trưng của chuỗi thời gian và (iii) nếu dữ liệu tính toán để đảm bảo nhận được các IMF’s đảm không có cực trị nhưng bao gồm chỉ duy nhất các bảo được ý nghĩa vật lý đầy đủ của nó. Thành điểm uốn, thì có thể lấy đạo hàm để tìm ra cực trị. phần thứ nhất c1(t) chứa những nội dung chính Khi đã xác định được các điểm cực trị, các điểm xác nhất, hoặc thông tin về tần số cao nhất tại cực đại được nối với nhau bởi một đường cong mỗi điểm thời gian. Phần dư sau quá trình sàng bậc ba và sẽ xác định được một đường bao trên; lọc đầu tiên là: tương tự từ các điểm cực tiểu cũng xác định r t = x t - c t (10) được một đường bao dưới. Đường bao trên và 1 1 đường bao dưới sẽ chứa tất cả các điểm dữ liệu Sau đó r1(t) được sử dụng để thay thế tín hiệu 26 Tạp chí KHCN Xây dựng – số 1/2016
3. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG chưa xử lý x(t) và quá trình sàng lọc tiếp tục để phần được cho bởi: tạo ra các IMF khác. Quá trình sàng lọc có thể n r t = x t - c t (11) dừng theo yêu cầu của quá trình vật lý. Tuy 1i=1 i nhiên, có những tiêu chuẩn chung, ví dụ tổng Như vậy mỗi ci(t) trong công thức (11) chính là bình phương của phần dư nhỏ hơn một giá trị một IMF, là kết quả của quá trình sàng lọc hoàn cho trước hoặc phần dư trở thành một hàm đơn toàn thỏa mãn các yêu cầu của phép biến đổi điệu. Phần dư sau quá trình sàng lọc n thành Hilbert. Hình 1. Tín hiệu ECG được phân tích thành các IMF‘s theo EMD 2.3 Sự gián đoạn của tín hiệu dao động EEMD được giới thiệu nhằm nâng cao hiệu (Intermittence) và EEMD quả của EMD để tránh hiện tượng chồng phổ [5]. Tín hiệu dao động luôn luôn bị ảnh hưởng Phương pháp này thực hiện quá trình lựa chọn của nhiễu và phần lớn có tính không ổn định, đây trên một quần thể các tín hiệu chứa nhiễu, được là một nguyên nhân gây ra sự gián đoạn của dữ tạo ra từ tín hiệu ban đầu cộng với chuỗi nhiễu liệu trong các scale thời gian. Sự gián đoạn là trắng (hay nhiễu Gaussian) với các hệ số khác nguyên nhân gây ra hiện tượng trộn lẫn dao động nhau, mỗi IMF của tín hiệu gốc được xác định (mode mixing/scale mixing), nghĩa là trong một chính là trị trung bình của quần thể các lần thử, IMF tồn tại các kiểu dao động khác nhau, hay trong đó có chứa tín hiệu cộng với nhiễu trắng cùng một kiểu dao động mà tồn tại trong hai IMF với biên độ hữu hạn nhưng không phải là vi phân. và hậu quả của nó là hiện tượng méo tín hiệu Với EEMD có thể tách hoàn toàn các scale của (aliaxing) và mất đi trị trung bình vật lý của tín tín hiệu một cách tự nhiên mà không cần bất kỳ hiệu, ảnh hưởng đến kết quả khi áp dụng HT. chuẩn lựa chọn chủ quan trước đó. Quá trình Huang và đồng nghiệp [3] đã giải quyết hiện EEMD gồm các bước: (i) Khởi tạo tín hiệu mới tượng trên bằng cách giới hạn kích thước của bằng cách cộng một chuỗi nhiễu trắng Gaussian; các scale, về thực chất là xem xét khoảng cách (ii) phân tích quần thể dữ liệu đã cộng nhiễu giới hạn của các cực đại liên tiếp trong một IMF thành các IMF’s dựa trên EMD; (iii) lập lại bước để sao cho IMF chỉ kết hợp các sóng mà khoảng (i) và (ii) nhưng với dãy nhiễu trắng khác trong cách của các cực đại liên tiếp nằm trong giới hạn mỗi lần, và sự lặp này phải đủ lớn để các thành cho phép của chúng. Ngoài ra cũng tiến hành tương tự đối với chiều dài của các zero-crossings phần nhiễu có thể loại bỏ lẫn nhau sau khi lấy để hạn chế ảnh hưởng của vấn đề gián đoạn trung bình của các IMF’s tương ứng. Về thực trong tín hiệu. chất EEMD là một bộ lọc tín hiệu [5]. Tạp chí KHCN Xây dựng – số 1/2016 27
4. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Hình 2. Tín hiệu ECG được phân tích thành các IMF’s theo EEMD Các công trình nghiên cứu hiện nay [5] đã njΦ t n j ω t dt x t = RP a t ei =RP  a t e i (17) chứng tỏ, EEMD là một phương pháp thích hợp i=1i i=1 i để xử lý các tín hiệu dao động không dừng và Trong (17) chứa cả biên độ và tần số của mỗi phi tuyến rất thường gặp trong chẩn đoán đánh IMF như các hàm của thời gian t. Mối quan hệ giá công trình, nhất là các kết cấu dưới của công giữa tần số - thời gian của biên độ được gọi là trình cầu. phổ Hilbert danh định (NHS) của tín hiệu: n j ω t dt 2.4 Phổ Hilbert danh định (NHS) và phổ Hilbert H ω,t = RP a t e i (18) i=1 i giới hạn (MHS) Phổ Hilbert giới hạn (MHS) của tín hiệu, được Áp dụng biến đổi Hilbert với mỗi IMF ci(t) định nghĩa như sau: nhận được (11), sẽ nhận được: T h ω,t = H ω,t dt (19) c  1 i 0 H c t = d (12) i π - t - Trong (19), T là tổng thời gian của tín hiệu. Từ (2) và (12), sẽ nhận được một tín hiệu Như vậy về ý nghĩa có thể thấy, NHS đưa ra một dạng giải tích: độ đo về sự đóng góp trong biên độ của mỗi tần jΦ t số theo thời gian, còn MHS là một phổ tần số đưa z t=c t+jHc t =a te i (13) i i i i ra một độ đo về sự đóng góp của mỗi tần số Trong đó: trong tổng biên độ của tín hiệu đang xét. 2 2 a t = c t + H c t (14) i i i Từ (15), xác định được hàm pha tổng của tín hiệu đang xem xét: H c t i Φ t = arctan (15) i H c t c t n n i i Φt=  Φ t=  arctan (20) i c t Từ (15) nhận được tần số tức thời của IMF tại i=1 i=1 i thời điểm t: Φ(t) là tổng số vòng quay của tín hiệu x(t) dΦ t trong mặt phẳng phức. Đối với một kết cấu, nó là i ω t = (16) một công cụ để trích xuất và biểu diễn pha tức i dt thời của các sóng nội tại di chuyển trong kết cấu Sau khi thực hiện biến đổi Hilbert cho mỗi của bất kỳ thông số động lực học nào có thể đo IMF, tín hiệu ban đầu được biểu diễn bởi phần đạc được của kết cấu. Đây là một thông số quan thực của biến đổi Hilbert: trọng để phát hiện hư hỏng trong kết cấu. 28 Tạp chí KHCN Xây dựng – số 1/2016
5. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG 3. Áp dụng HHT để chẩn đoán hư hỏng trong trưng dao động (có thể là các tần số dao kết cấu dưới của công trình cầu động riêng) được xem như một mốc so 3.1 Một số luận cứ và phạm vi áp dụng HHT sánh của kết cấu cần chẩn đoán. trong chẩn đoán hư hỏng kết cấu xây dựng  Sau mùa lũ, tiến hành đo và thu nhận tín Như đã trình bày ở trên, biến đổi Hilbert- hiệu dao động của kết cấu cần chẩn đoán Huang là một kỹ thuật biến đổi phi tuyến có nhiều dưới tác dụng của các loại xe có tải trọng ưu điểm và rất thích hợp với xử lý các tín hiệu phi khác nhau. tuyến và không dừng. Các luận cứ của việc áp  Sử dụng EEMD để phân tích các dao dụng HHT trong chẩn đoán kết cấu xây dựng đã động nhận được khi đo dao động thành được thảo luận trong nhiều nghiên cứu [1, 4, 6, các IMF‘s. 7]. Có thể tóm tắt ngắn gọn như sau: Ứng xử  Tiến hành tính tần số tức thời, pha tức (tĩnh học và động học) của bất kỳ một kết cấu xây thời và các đặc trưng phổ Hilbert (NHS và dựng nào khi chịu tác động của tải trọng thiết kế MHS) của các hàm IMF‘s. sẽ là tuyến tính và đàn hồi. Khi kết cấu bị hư  Đánh giá sự thay đổi của các sóng nội tại hỏng thì độ bền của nó bị giảm xuống, và khi chịu qua các pha tức thời nhận được để tìm tác dụng của các tải trọng, thậm chí còn nhỏ hơn sự thay đổi phi tuyến của dao động qua tải trọng thiết kế, thì ứng xử của kết cấu sẽ xuất đó nhận được dấu hiệu hư hỏng của kết hiện tính phi tuyến. Sự phi tuyến của ứng xử của cấu. kết cấu có thể là hệ quả của các biến dạng không đàn hồi của vật liệu hoặc có thể là do sự thay đổi  Nếu cầu có nhiều trụ, so sánh phổ Hilbert độ cứng của kết cấu do bị hư hỏng. Dao động đo giới hạn (MHS) của các IMF của các trụ được trên kết cấu gây ra bởi một kích động cho khác nhau dưới cùng một tải trọng xe trước cần được tách ra thành những tần số riêng chạy qua để đánh giá độ suy giảm của phần như tần số của lực kích động, tần số dao MHS và tìm ra các trụ có sự suy giảm độ động tự nhiên của toàn bộ kết cấu hoặc của các cứng. kết cấu cục bộ, hoặc cả hai và các tần số  Nếu cầu có một trụ, so sánh MHS của các nhiễu để xem xét trong miền thời gian. Do HHT IMF‘s của trụ dưới các tải trọng xe chạy là công cụ duy nhất (cho đến hiện nay) có thể qua khác nhau để tìm mức độ và sự thay phân tích các tín hiệu phi tuyến và không dừng, đổi phi tuyến của dao động dưới mỗi loại nên sử dụng HHT có thể xác định được những tải trọng qua đó đánh giá độ suy giảm độ thay đổi của các đặc tính dao động đo được trên cứng của trụ. kết cấu theo miền thời gian, qua đó xác định  Tính toán tỷ số tần số trước và sau khi lũ được vị trí và mức độ của hư hỏng. và dựa vào các tiêu chí đánh giá để xác Theo luận cứ trên, một phương pháp khả thi định kết cấu dưới còn an toàn hoặc đã bị để phát hiện vị trí và mức độ hư hỏng của cầu là suy giảm. sử dụng tải trọng tức thời di động qua cầu, ghi  Nếu kết cấu không an toàn, cần tiến hành nhận và sử dụng HHT để xử lý các tín hiệu dao các bài toán nhận dạng hoặc các phương động ghi nhận đươc. Khi tải trọng này tác dụng pháp trình bày trong [1] để chẩn đoán các trực tiếp lên các điểm hư hỏng thì ứng xử động thông số độ cứng của kết cấu. học của cầu sẽ thay đổi mạnh nhất [1, 4]. 3.3 Một ví dụ thực tế 3.2 Đề xuất trình tự tiến hành Khảo sát một cầu có 2 làn xe, dầm BTCT với Dựa trên các nghiên cứu về HHT trong đánh sơ đồ 3 nhịp giản đơn. Sử dụng đầu đo gia tốc giá kết cấu, trong [1] đã đề xuất trình tự tiến hành theo 3 phương để đo dao động trụ P1 và P2 khi đối với kết cấu dưới của công trình cầu trước và có xe tải nặng qua cầu. Kết quả đo thể hiện trên sau khi có lũ, có thể tóm tắt như sau: hình 3. Dữ liệu này không chỉ ra được bất kỳ sự  Trước mùa lũ, tiến hành thu thập các khác biệt lớn nào giữa hai trụ khi chưa xử lý. Qua thông tin, các kết quả thí nghiệm, đo đạc, phân tích HHT và FFT, nhận thấy sự khác nhau thử tải đã tiến hành để xác định các đặc khá lớn giữa HHT và FFT như trên hình 4. Tạp chí KHCN Xây dựng – số 1/2016 29
6. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Hình 3. Gia tốc đo của trụ P1 và P2 Hình 4. Kết quả xử lý FFT và HHT Kết quả của FFT không thể hiện rõ sự khác 4. Kết luận và kiến nghị nhau về dao động của hai trụ, FFT chỉ cho thấy Qua những kết quả trình bày ở trên có thể rút có sự thay đổi nhỏ về biên độ, nhưng không thể ra một số kết luận và kiến nghị như sau: đưa ra được nhận xét gì với sự thay đổi biên độ này. Trong khi đó kết quả của HHT, MHS lại chỉ  Biến đổi Hilbert-Huang là một phương ra được sự thay đổi mạnh về độ cứng của hai trụ. pháp biến đổi có nhiều ưu điểm: Có thuật Điểm nhọn tần số của trụ P2 là 1,2 Hz và của trụ toán đơn giản, là phương pháp thích nghi P1 là 1,6 Hz cho thấy một sự suy giảm rõ ràng dựa trên đặc tính nội tại của dữ liệu/tín hiệu của tần số và là dấu hiệu giảm độ cứng rõ ràng và rất hiệu quả để các xử lý các dữ liệu/tín của kết cấu trụ P2. Vì tần số tỉ lệ với căn bậc hai hiệu phi tuyến, không dừng, dải rộng và bị của độ cứng kết cấu, nên tỉ số 1,6/1,2 chứng tỏ ảnh hưởng bởi nhiễu. rằng độ cứng trụ P2 đã giảm khoảng 50%. Kiểm  Đối với công trình nói chung, công trình tra trực tiếp thấy rằng trụ P2 bị xói khá nhiều sau cầu và kết cấu dưới của nó nói riêng, các tín trận lũ lớn. hiệu dao động ghi nhận về ứng xử động của 30 Tạp chí KHCN Xây dựng – số 1/2016
7. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG chúng thường là các tín hiệu phi tuyến, Proc. R. Soc. Lond , 454, pp. 903-995, 1998. không dừng, dải rộng và bị ảnh hưởng khá [4]. Larry D. and Olson P. E. (2005), Dynamic lớn bởi nhiễu. Việc xử lý các tín hiệu này Bridge Substructure Evaluation and Monitoring, theo các biến đổi cũ (như FFT, WT, HT ) Pub. No. FHWA-RD-03-089, September 2005. không mang lại các kết quả mong muốn. Những kết quả nghiên cứu đã thực hiện ở [5]. Wu Z. and Huang N. E. (2009), “Ensemble Hoa Kỳ, Việt Nam cho thấy HHT là một công empirical mode decomposition: A noise-assisted cụ thích hợp để chẩn đoán các hư hỏng data analysis method”, The Advances in trong kết cấu dưới của công trình cầu và cần Adaptive Data Analysis, Vol. 1, No. 1, pp. 1-41. được tiếp tục nghiên cứu để hoàn thiện. [6]. Yang J. N., Lei Y., Pan S. W., and Huang N.E. TÀI LIỆU THAM KHẢO (2003), “System identication of linear structures based on Hilbert-Huang spectral analysis; Part I: [1]. Bùi Đức Chính (2010), Nghiên cứu áp dụng Normal modes”, Earthquake Engineering and công nghệ và thiết bị mới trong kiểm tra, đánh Structural Dynamics, Vol. 32 (10), pp. 1443-1467. giá và sửa chữa kết cấu dưới của các công trình cầu ở Việt Nam, Đề tài KH và CN trọng điểm cấp [7]. Yang J. N., Lei Y., Pan S. W., and Huang N.E. Bộ 2009-2010, Hà Nội, 2010. (2003), “System Identification of Linear Structures based on Hilbert-Huang Spectral [2]. Cohen L. (1995), Time-Frequency Analysis : Analysis. Part II: Complex Modes”, Earthquake Theory and Applications, Prentice Hall, Inc., Engineering & Structural Dynamics, Vol. 32 (10), Upper saddle River, NJ, 1995. pp. 1533-1554. [3]. Huang N. E. et al (1998), “The Empirical Mode De- Ngày nhận bài: 13/02/2016. composition Method and the Hilbert Spectrum for Non-linear and Non-stationary Time Series Analysis”, Ngày nhận bài sửa lần cuối: 26/02/2016. Tạp chí KHCN Xây dựng – số 1/2016 31