Chương trình phân tích một số dị thường trọng lực bằng mạng perceptron trong môi trường matlab

Journal of Science – 2015, Vol. 8 (4), 22 – 31 Part D: Natural Sciences, Technology and Environment 22 CHƯƠNG TRÌNH PHÂN TÍCH MỘT SỐ DỊ THƯỜNG TRỌNG LỰC BẰNG MẠNG PERCEPTRON TRONG MÔI TRƯỜNG MATLAB Nguyễn Hồng Hải1 1ThS. Trường Đại học An Giang Thông tin chung: Ngày nhận bài: 22/09/14 Ngày nhận kết quả bình duyệt: 15/01/15 Ngày chấp nhận đăng: 12/15 Title: The analytic program of gravity anomalies by perceptron network in matlab Từ khóa: Dị thường trọng lực, dị thường Bouguer, mạng perceptron, thăm dò từ và trọng lực, Matlab Keywords: Gravity anomalies, Bouguer anomalies, perceptron network, magnetic and gravity survey, Matlab ABSTRACT The program was developed in order to define the crystal basement beneath the deposits by the values of gravity anomalies measured on the ground without the density sign of the deposits. The analytic program with friendly display can calculate the data of the survey line’s average density sign, the squared average error and the calculation time of program. Besides, it can display graphs of observed gravity anomalies, calculated gravity anomalies and graphs of depth, shape and size of a foreign body. The result of the model program and the analytic results of 5 negative Bouguer anomalies in the Mekong Delta including Chau Doc, Tam Nong, Thap Muoi, Bac Lieu and Ca Mau were consistent with the experimental measurements. TÓM TẮT Chương trình được xây dựng nhằm xác định mặt móng kết tinh nằm bên dưới các lớp trầm tích, từ các giá trị dị thường trọng lực đo trên mặt đất khi không biết hiệu mật độ của các lớp trầm tích. Chương trình có giao diện thân thiện với người dùng, tính toán được hiệu mật độ trung bình của tuyến khảo sát, sai số trung bình bình phương và thời gian tính toán. Bên cạnh đó, chương trình cho phép hiển thị đồ thị dị thường trọng lực quan sát, dị thường trọng lực tính được; đồ thị về độ sâu, hình dạng và kích thước của dị vật. Kết quả trên mô hình và năm dị thường Bouguer âm ở vùng Đồng bằng sông Cửu Long: Châu Đốc, Tam Nông, Tháp Mười, Bạc Liêu, Cà Mau phù hợp với thực nghiệm. 1. MỞ ĐẦU Thăm dò trọng lực là phương pháp Địa Vật lý nghiên cứu vỏ trái đất và thăm dò các khoáng sản có ích, dựa trên việc nghiên cứu sự phân bố của trường trọng lực trên mặt đất. Thăm dò trọng lực có hiệu suất quan sát thực địa tương đối cao. Trong một số trường hợp, việc giải thích kết quả thăm dò trọng lực có thể nhận được kết quả định lượng với độ chính xác tin cậy. Trong thăm dò trọng lực, có 2 loại bài toán: bài toán thuận và bài toán ngược. Nguyên tắc của bài toán ngược là trường được đem phân chia, giải đoán phải đúng là trường của chính dị vật đang quan tâm gây ra. Do đó, cần phải cố gắng loại trừ tối đa các nhiễu. Tuy nhiên không bao giờ đạt được kết quả tuyệt đối. Nếu như bài toán thuận có nghiệm duy nhất, thì bài toán ngược nói chung là đa trị. Do đó, việc giải bài toán ngược thăm dò trọng lực phụ thuộc vào sự phân loại và nhận dạng các yếu tố của bài toán với một lượng lớn các điều kiện phi tuyến; điều đó đòi hỏi cần phải được huấn luyện cách tiếp cận đã tích hợp các điều kiện phi tuyến để giải quyết vấn đề. Journal of Science – 2015, Vol. 8 (4), 22 – 31 Part D: Natural Sciences, Technology and Environment 23 Trong khi đó, nhờ vào khả năng học, lưu lại và khái quát hoá từ các mẫu được huấn luyện hoặc dữ liệu thu thập, mạng neural nhân tạo (Artificial Neural Network: ANNs) trở thành một phát minh mới đầy hứa hẹn của hệ thống xử lý thông tin. Các phép tính toán của neural cho phép giải quyết tốt những bài toán đặc trưng bởi một số tính chất như: sử dụng không gian nhiều chiều, các tương tác phức tạp, chưa biết hoặc không thể theo dõi về mặt toán học giữa các biến. Ngoài ra, phương pháp ANNs còn cho phép tìm ra nghiệm đơn trị của bài toán. Với tầm quan trọng của thăm dò trọng lực trong nghiên cứu thăm dò Địa Vật lý và các tiện ích của mạng perceptron – một loại mạng neural đơn, tôi quyết định nghiên cứu xây dựng “Chương trình phân tích một số dị thường trọng lực bằng mạng perceptron trong môi trường Matlab”, chương trình được xây dựng nhằm xác định mặt móng kết tinh nằm bên dưới các lớp trầm tích, từ các giá trị dị thường trọng lực đo trên mặt đất khi không biết hiệu mật độ của các lớp trầm tích. 2. PHƯƠNG PHÁP 2.1 Phương pháp tiến Nguyên tắc của phương pháp tiến là dựa trên các hiểu biết về địa chất và địa vật lý của vùng khảo sát, mô hình ban đầu của nguồn gây ra dị thường được xây dựng qua trực giác của người phân tích. Tính giá trị dị thường lý thuyết gây ra bởi mô hình này và đem so sánh với dị thường quan sát để hiệu chỉnh các tham số của mô hình sao cho hai dị thường tương thích nhau. Quá trình tính toán gồm ba bước (Blakely, 1995): - Khởi tạo mô hình và tính dị thường của mô hình. - So sánh dị thường của mô hình và dị thường quan sát qua sai số trung bình bình phương để hiệu chỉnh mô hình. - Hiệu chỉnh mô hình nhiều lần cho tới khi sai số trung bình bình phương đạt giá trị mong muốn (đạt giới hạn cho phép) Sơ đồ khối của phương pháp tiến được trình bày trong Hình 1. Hình 1. Ba bước của phương pháp tiến giải đoán tài liệu trọng lực. Trong đó: A là dị thường quan sát; A0 là dị thường tính được; p1, p2,... là các thuộc tính của nguồn gây ra dị thường như độ sâu, độ dày, mật độ,... Journal of Science – 2015, Vol. 8 (4), 22 – 31 Part D: Natural Sciences, Technology and Environment 24 2.2 Mô hình của phương pháp Compact Phương pháp này được đưa ra bởi Last và Kubik (1983) để giải bài toán ngược trọng lực bằng phương pháp bình phương tối thiểu, đây là phương pháp thử sai nhưng cho kết quả nhanh và chính xác cao. Mô hình được sử dụng là mô hình hai chiều (2D) bao gồm các ô hình chữ nhật cố định chứa các giá trị hiệu mật độ như trong Hình 2. Hình 2. Mô hình 2D gồm các khối hình chữ nhật chứa các giá trị hiệu mật độ Hình 3. Dị thường trọng lực gây ra bởi một hình chữ nhật thứ j lên một điểm quan sát thứ i Khi đó dị thường trọng lực do tất cả các khối hình chữ nhật gây ra tại điểm quan sát thứ i được tính bởi công thức: M i ij j j 1 g a .    , i = 1,...N (1) Ở đây j là hiệu mật độ tại khối thứ j; aij là ma trận phần tử đại diện cho ảnh hưởng của các khối thứ j lên trọng lực tại điểm i. Công thức để tính giá trị aij cho bởi:         2 3 4 i j 1 4 3 ij hd j 4 2 j 3 1 r r rdx x log d log 2 r r r a 2G h hz z 2 2                            (2) trong công thức trên, r1, r2, r3, r4 là khoảng cách từ điểm i đến các đỉnh của hình chữ nhật thứ j. Journal of Science – 2015, Vol. 8 (4), 22 – 31 Part D: Natural Sciences, Technology and Environment 25 2 2 2 1 2 2 2 2 2 2 2 3 2 2 2 4 ( / 2) ( / 2) ( / 2) ( / 2) ( / 2) ( / 2) ( / 2) ( / 2) j i j j i j j i j j i j r z h x x d r z h x x d r z h x x d r z h x x d                         (3) và θ1, θ3, θ3, θ4 là góc hợp bởi các cạnh r1, r2, r3, r4 với mặt đất: 1 2 3 4 arctan( / 2) / ( / 2) arctan( / 2) / ( / 2) arctan( / 2) / ( / 2) arctan( / 2) / ( / 2) i j j i j j i j j i j j x x d z h x x d z h x x d z h x x d z h                        (4) Trong các công thức (2), Ghd là hằng số hấp dẫn; h, d là bề rộng và bề ngang của mỗi hình chữ nhật; zj là độ sâu của ô thứ j; xi là toạ độ của điểm quan sát P, xj là toạ độ chiều ngang của ô thứ j. Như vậy, mỗi một ô chữ nhật sẽ chứa một giá trị hiệu mật độ vj và ứng với mỗi điểm quan sát Pi ta sẽ thu được một giá trị dị thường trọng lực do tất cả các ô này gây ra. 3. XÂY DỰNG CHƯƠNG TRÌNH 3.1 Bài toán Bài toán ở đây chỉ giới hạn ở bài toán 2D. Đó là bài toán xác định độ sâu đến mặt móng nằm bên dưới các lớp trầm tích của một bồn trầm tích khi có giá đo của dị thường Bouguer trên một tuyến đo. Để giải quyết vấn đề này, tôi giả sử, mật độ của các lớp trầm tích là đồng nhất và mô hình mặt cắt của bồn trầm tích được xấp xỉ bằng một tập hợp gồm N các tấm chữ nhật (vô hạn theo phương y thẳng góc với tuyến đo) có các cạnh lần lượt song song với trục x (tuyến đo) và trục z (độ sâu); các điểm đo được đặt tại trung điểm của cạnh trên của mỗi tấm chữ nhật. Khi đó, dữ liệu là N giá trị dị thường trọng lực Bouguer (miligal) quan sát trên một tuyến đo có chiều dài L (km), khoảng cách các điểm đo là d (km). 3.2 Tóm tắt thuật giải Để áp dụng thuật giải này vào việc giải bài toán ngược trọng lực tôi sử dụng ba bước tính của phương pháp tiến với mô hình Compact; trong đó, bước tính thứ hai và thứ ba sử dụng thuật toán perceptron (Nguyễn Hồng Hải, 2013). Sơ đồ của thuật giải mô tả ở Hình 4. Tương ứng với M đầu vào sẽ có M trọng số wj cho mạng perceptron. Để áp dụng mạng perceptron vào việc giải bài toán ngược trọng lực, các giá trị trọng số này được lấy là hiệu mật độ j của từng ô hình chữ nhật, chúng được khởi tạo một cách ngẫu nhiên hoặc chọn bằng 1 cho tất cả các ô. Journal of Science – 2015, Vol. 8 (4), 22 – 31 Part D: Natural Sciences, Technology and Environment 26 Hình 4. Sơ đồ khối của phương pháp mạng perceptron trong giải bài toán ngược trọng lực Hình 5. Cấu tạo của mạng neural để giải bài toán ngược trọng lực Thuật toán perceptron áp dụng vào bài toán ngược trọng lực được trình bày trong Hình 5 và mục đích của phương pháp là tìm hình dạng dị vật bằng cách điều chỉnh mật độ của các hình chữ nhật trong mô hình. Mạng perceptron sẽ học theo từng mẫu một. Điều này giúp cho quá trình học và cập nhật trọng số của mạng sẽ nhanh hơn (Nguyen Ngoc Thanh Son, Nguyen Hong Hai, Luong Phuoc Toan & Dang Van Liet, 2012). Trước tiên, giá trị đầu ra của mạng sẽ được tính thông qua quá trình lan truyền tiến. Cột thứ nhất của ma trận đầu vào sẽ nhân với ma trận trọng số Journal of Science – 2015, Vol. 8 (4), 22 – 31 Part D: Natural Sciences, Technology and Environment 27 (là giá trị hiệu mật độ). Giá trị của đầu ra là tổng các giá trị này. Tiếp theo sẽ là quá trình lan truyền ngược. Các giá trị trọng số sẽ được cập nhật bằng cách cộng thêm đạo hàm bậc nhất của bình phương sai số (hàm lỗi). Quá trình này sẽ được lặp lại cho cột thứ hai của ma trận đầu vào. Sau khi đã học xong N mẫu, nếu sai số trung bình bình phương vẫn chưa đạt được giá trị mong muốn thì mạng sẽ quay trở lại mẫu (cột) đầu tiên và tiếp tục học cho đến khi sai số đúng với giá trị yêu cầu. 3.3 Chương trình thuật giải Trên cơ sở phương pháp của thuật giải, tôi tiến hành xây dựng các hàm và chương trình thuật giải được xây dựng cho một dị thường cụ thể (ví dụ: Dị thường Bouguer âm Châu Đốc) bao gồm tính toán số liệu và vẽ đồ thị được lập trình như sau: clc clear all; % xuat file du lieu di thuong o Chau Doc load gAngiang.m t=gAngiang; % xac lap chieu dai cua o chu nhat dh=0.1; % xac lap chieu rong (sau) cua o chu nhat dd=1; % xac lap he so hoc LR=0.01; % xac lap sai so cho phep ErTh=0.001; % xac lap nguong luong tu hoa hieu mat do nguong=-0.32; sodiemdo=length(t); % xac lap kiem tra trong khoang 3km dosau=3; % tinh gia tri Ar cua moi diem g(i) [Ar,k,x,z]=TinhAr(dh,dd,sodiemdo,dosau); [rows,numdatasets]=size(Ar); % moi numdatasets la mot tap hop gia tri Ar de tinh g(i) % tao ma tran trong so (hieu mat do) weights=NN2weights(rows); TDIN=Ar; TDOUT=t; ketqua=ketqua(weights,TDIN,TDOUT,LR,ErTh,nguong); giatrig=thu(ketqua.inputs,Ar); Journal of Science – 2015, Vol. 8 (4), 22 – 31 Part D: Natural Sciences, Technology and Environment 28 figure(1) plot(t,'r') hold on plot(giatrig,'b*') trongso=reshape(ketqua.inputs,k,sodiemdo); figure(2) pcolor(x,-z,trongso); 4. KẾT QUẢ VÀ THẢO LUẬN Dựa trên chương trình đã xây dựng bằng phần mềm Matlab 7.10.0.499 (R2010a), chạy trên nền hệ điều hành Windows Vistavới Processor: Intel(R) Core(TM) 2 Duo2.6GHz; RAM: 1GB; Hard Disk: 250GB., tôi tiến hành phân tích 5 dị thường trọng lực Bouguer ở Đồng bằng sông Cửu Long Hình 6. Hình 6. Bản đồ dị thường trọng lực Bouguer vùng Đồng bằng sông Cửu Long 1. Dị thường âm Châu Đốc (An Giang) 2. Dị thường âm Tam Nông (Đồng Tháp) 3. Dị thường âm Tháp Mười (Đồng Tháp) 4. Dị thường âm Bạc Liêu 5. Dị thường âm Cà Mau Kết quả tính toán dị thường Bouguer Tam Nông (Đồng Tháp) (Hình 7) cho thấy: Hiệu mật độ tìm được Δρ = -0,40691 g/cm3, sai số trung bình bình phương E = 2,8597.10-6, thời gian thực hiện chương trình t = 6,412 s. Đồ thị cho thấy mặt móng kết tinh kéo dài từ phía Tây ở độ sâu 0,3 km và sâu dần đến độ sâu cực đại 1 km và dốc ngược về phía Đông để đạt độ sâu khoảng 0,5 km, độ sâu trung bình của toàn vùng là 0,7 km. Journal of Science – 2015, Vol. 8 (4), 22 – 31 Part D: Natural Sciences, Technology and Environment 29 Hình 7. Kết quả tính toán của dị thường Bouguer Tam Nông (Đồng Tháp) a) độ sâu b) dị thường trọng lực Hình 8. Kết quả phân tích dị thường trọng lực Tam Nông tính bằng thuật giải di truyền (Bùi Thị Nhanh, 2014) Với kết quả tính bằng thuật giải di truyền (Hình 8) (Bùi Thị Nhanh, 2014) cho kết quả gần giống với kết quả trên; tuy nhiên, thời gian thực hiện chương trình lâu hơn t = 23 s với máy tính (chương trình được chạy bằng phần mềm Matlab 7.10 (R2010a), chạy trên nền hệ điều hành Windows 7 ultimate với Processor: Intel(R) Core(TM) i3-3227U CPU @ 1.90GHz; RAM 4.00GB) và sai số trung bình bình phương tối thiểu lớn hơn E = 0,02. So sánh kết quả phân tích năm vùng dị thường thuộc Đồng bằng sông Cửu Long bằng mạng perceptron và thuật giải di truyền được trình bày ở Bảng 1. Kết quả tính toán mật độ ở các vùng khác nhau so với giá trị mật độ trung bình tính từ ba giếng khoan sâu trong vùng là 455 kg/m3 không sai biệt nhiều. Riêng vùng Cà Mau, hiệu mật độ tính được khá nhỏ so với hiệu mật độ trung bình của vùng, đây có thể dị thường Cà Mau nằm trên vùng than bùn nên có mật độ nhỏ. 0 5 10 15 20 25 30 35 40 45 -1 -0.8 -0.6 -0.4 -0.2 0 x (Km) § é s© u (K m ) 0 5 10 15 20 25 30 35 40 45 -16 -14 -12 -10 -8 -6 -4 x (Km) de lt a g (m ga l) Journal of Science – 2015, Vol. 8 (4), 22 – 31 Part D: Natural Sciences, Technology and Environment 30 Bảng 1. Bảng so sánh kết quả phân tích giữa mạng perceptron và thuật giải di truyền Dị thường Châu Đốc Tam Nông Tháp Mười Bạc Liêu Cà Mau Mạng perceptron Δρ (kg/m3) -434,08 -406,91 -414,39 -407,04 -397,18 E 2,2951.10-6 2,8597.10-6 1,7033.10-5 1,4003.10-5 1,8856.10-5 t (giây) 5,913 6,412 6,38 8,33 4,04 Thuật giải di truyền Δρ (kg/m3) -428 -407 -409 -406 -402 E 0,03 0,02 0,03 0,05 0,05 t (giây) 31 23 11 29 13 Bảng 2. Độ sâu cực đại của mặt móng ở năm vùng dị thường thuộc Đồng bằng sông Cửu Long Tên vùng Châu Đốc Tam Nông Tháp Mười Bạc Liêu Cà Mau Δgmax (mGal) -26 -15 -17,7 -20 -10.1 zmax (km) -1,5 -1,0 -1,1 -1,4 -0,7 ztrung bình (km) - 1,0 -0,7 0,7 -0,8 -0,4 Bảng 2 là độ sâu của mặt móng tính được tại một số vùng của Đồng bằng sông Cửu Long. So sánh kết quả trên với độ sâu của các giếng khoan sâu đến mặt móng trong vùng, trừ giếng khoan Cửu Long ở ven biển có độ sâu tới 2120 m, hai giếng khoan còn lại nằm ở phần đất liền của vùng Đồng bằng sông Cửu Long là Cần Thơ (797 m) và Long Xuyên (350 m), kết quả cho thấy độ sâu trung bình từ 0,4 km đến 1 km cho thấy kết quả tính được và độ sâu đo thực nghiệm tương đối phù hợp. 5. KẾT LUẬN Đề tài đã xây dựng phương pháp mới, ứng dụng phù hợp thuật giải mạng perceptron vào phương pháp tiến, có kết hợp mô hình Compact để giải bài toán ngược trọng lực. Điểm đặc sắc: giải bài toán trong trường hợp không thể biết hiệu mật độ; tốc độ giải nhanh, sai số trung bình bình phương thấp. Chương trình có giao diện thân thiện với người sử dụng; chương trình có khả năng luyện mạng bằng cách so sánh dị thường quan sát và dị thường tính toán, tiến tới tính toán được số liệu về: hiệu mật độ trung bình của tuyến khảo sát, sai số trung bình bình phương, thời gian tính toán của chương trình; bên cạnh đó, cho phép hiển thị đồ thị dị thường trọng lực quan sát, dị thường trọng lực tính được; đồ thị về độ sâu, hình dạng và kích thước của dị vật. Qua kiểm tra kết quả phân tích năm dị thường Bouguer âm ở vùng Đồng bằng sông Cửu Long: Châu Đốc, Tam Nông, Tháp Mười, Bạc Liêu, Cà Mau cho kết quả phù hợp với đo đạc thực nghiệm, có khả năng khai thác, triển khai ứng dụng trong các nghiên cứu về thăm dò từ và trọng lực. TÀI LIỆU THAM KHẢO Blakely, R. J. (1995). Potential theory in gravity and magnetic applications. USA: Cambridge University Press. Bùi Thị Nhanh. (2014). Phân tích một số dị thường trọng lực ở vùng Đồng bằng sông Cửu Long bằng thuật toán di truyền. (Luận văn Journal of Science – 2015, Vol. 8 (4), 22 – 31 Part D: Natural Sciences, Technology and Environment 31 thạc sỹ không xuất bản). Trường Đại học Cần Thơ, Cần Thơ, Việt Nam. Last, B. J. & Kubik, K. (1983). Compact gravity inversion. Geophysics,Vol.48, pp 713-721. Nguyễn Hồng Hải. (2013). Ứng dụng phương pháp Compact trong giải bài toán ngược trọng lực. Tạp chí Khoa học Trường Đại học Phú Yên, 4 (09/2013), tr.25-33, ISSN0866 – 7780. Nguyen Ngoc Thanh Son., Nguyen Hong Hai., Luong Phuoc Toan., & Dang Van Liet. (2012). Determination of the crytal basement of some gravity anomalies in Mekong delta area using the forced neural network. Proceedings of the International Scientific conference “Geophysics – Cooperation and Sustainable development”, pp. 371-379. Publishing house for Science and Technology, ISBN: 978-604- 913-091-5.