Open Access 
Abstract
Traditional methods of solving geophysical inverse problems often used an initial structural model and an iterative process to find solutions. During each iteration, the model was linearized and updated, then the error between the model's response and the observed data were calculated and minimized. The inverse results solved by this method were often strongly dependent on the initial model choice, and the minimum of the objective function found could be the local value but not the global minimum. Global search methods such as genetic algorithm (GA) overcame these limitations of traditional solutions. This paper presented the application of modified differential evolution (MDE) algorithm which was developed from the genetic algorithm, to one-dimension (1D) magnetotelluric inversion. The algorithm had been used to invert two 3-layer geoelectrical models. The results showed that when using a 3-layer model for inversion, the obtained parameters of the subsurface almost coincided with those of the theoretical models. In the case the number of layers used for inversion had more than three layers (four layers and five layers), the obtained results showed that there was still a good agreement between the inverted structure and the theoretical model. Therefore, the algorithm was applied to one-dimensional (1D) inversion for data of one magnetotelluric station measured at Cu Chi district, Ho Chi Minh city. The results showed that the subsurface structure down to a depth of 16 km consisted of 3 layers: the top layer had a resistivity of 119 Ω.m with a thickness of about 760 m, the middle resistive layer with a resistivity of 4353 Ω.m with about 10 km thickness, and the lower half-space with a relatively low resistivity of 41.4 Ω.m. This result was consistent with the geoelectric structure information obtained from other magnetotelluric studies carried out in the studied area, demonstrating the practical applicability of the improved differential evolution algorithm.
GIỚI THIỆU
Giải bài toán ngược hay nghịch đảo trong địa vật lý là xuất phát từ các số liệu dị thường hay các trường vật lý (thường được đo đạc trên mặt đất), tính toán để thu được các thông số hình học (gồm hình dạng, thế nằm, kích thước, độ sâu) và tính chất vật lý như điện trở suất, khối lượng riêng, vận tốc truyền sóng… của các đối tượng nghiên cứu nằm bên dưới mặt đất. Bài toán ngược thăm dò MT, cũng như phần lớn các bài toán ngược địa vật lý khác, là loại bài toán phi tuyến tức là các thông số của đối tượng nghiên cứu không có mối liên hệ tuyến tính với trường vật lý hay dị thường do chúng gây ra. Cách giải bài toán ngược thông thường là (thí dụ trong 1 ) sử dụng một mô hình ban đầu về đối tượng nghiên cứu, sau đó thông qua quá trình lặp để điều chỉnh các thông số của mô hình sao cho sai số giữa số liệu tính từ mô hình và số liệu số quan sát đạt đến một giá trị ngưỡng xác định hoặc số lần lặp đạt đến một số lần cho trước. Sai số giữa các số liệu tính toán và quan sát được biểu diễn bằng một hàm sai số (thường là hàm phi tuyến) và việc cực tiểu hàm này được thực hiện bằng cách tuyến tính hóa quanh một mô hình ban đầu thông qua quá trình lặp. Việc tuyến tính hóa bao gồm tính đạo hàm bậc nhất theo các tham số mô hình để thu được ma trận Jacobi dùng để cập nhật mô hình sau mỗi lần lặp. Kết quả nghịch đảo theo cách giải truyền thống này bị phụ thuộc mạnh vào việc lựa chọn mô hình xuất phát dùng để nghịch đảo. Tuy nhiên việc chọn được mô hình xuất phát một cách phù hợp không phải là điều dễ dàng. Thêm nữa, việc tìm kiếm giá trị cực tiểu của hàm sai số bằng cách tuyến tính hóa có thể dẫn đến tìm được các cực tiểu địa phương mà không phải cực tiểu toàn cục.
Trong những năm gần đây, lớp các bài toán tiến hóa phỏng theo thuyết tiến hóa sinh học của Darwin đã được nghiên cứu và phát triển cho các vấn đề khoa học và kỹ thuật nói chung và lĩnh vực địa vật lý nói riêng (thí dụ trong luận án tiến sĩ vật lý của Lương Phước Toàn do PGS. Đặng Văn Liệt và PGS. Trần Vĩnh Tuân hướng dẫn 2 ). Thuật toán tiến hóa được giới thiệu lần đầu tiên vào năm 1975 bởi tác giả Holland 3 . Ngày nay nó đã được phát triển theo nhiều hướng, có thể kể đến một số hướng chính như di truyền, bầy đàn (mô phỏng hành vi tìm kiếm thức ăn của đàn kiến, đàn ong…) 4 . Các thuật toán tiến hóa không đòi hỏi phải tuyến tính hóa hàm mô hình và không dùng cách tính đạo hàm trong quá trình tìm mô hình kết quả. Ưu điểm của chúng là tìm kiếm cực trị toàn cục và giải quyết được các bài toán có mức độ phi tuyến cao. Bài báo này giới thiệu việc áp dụng thuật toán tiến hóa vi sai cải tiến để giải bài toán ngược MT trường hợp một chiều đối với mô hình lý thuyết và số liệu thực tế.
PHƯƠNG PHÁP
Thuật toán tiến hóa vi sai cải tiến (Modified Differential Evolution, MDE)
Thuật toán tiến hóa vi sai, gồm 5 bước: khởi tạo, đột biến, lai ghép, chọn lọc và hội tụ. (Differential Evolution, DE) do Storn và Price đề nghị năm 1997 5 . Việc cải tiến thuật toán DE thành thuật toán MDE đã được thực hiện bởi Sơn và đồng nghiệp 6 , 7 . Các tác giả này đã thực hiện cải tiến ở bước đột biến của thuật toán DE bằng cách sử dụng hai kiểu đột biến rand/1 và best/1 để tạo vector đột biến thay vì chỉ sử dụng một kiểu đột biến rand/1 hoặc best/1 như trong thuật toán DE ban đầu. Với kiểu đột biến rand/1, thuật toán mạnh về tìm kiếm toàn cục, trong khi đột biến best/1, thuật toán mạnh về tìm kiếm cục bộ. Sự cải tiến này giúp cho sự tìm kiếm toàn cục và tìm kiếm cục bộ trong thuật toán MDE trở nên cân bằng hơn so với thuật toán DE. Trong thuật toán MDE, việc áp dụng kiểu đột biến rand/1 được áp dụng khi giá trị ngẫu nhiên trong thuật toán (phân bố đều trong phạm vi [0,1]) được tạo ra có giá trị lớn hơn 0,3; còn nếu giá trị ngẫu nhiên trong thuật toán được tạo ra nhỏ hơn hoặc bằng 0,3 thì áp dụng kiểu đột biến best/1.
Nếu tỉ lệ giữa 2 kiểu đột biến này, rand/1 và best/1, bằng 7/3 thì sự tìm kiếm toàn cục và tìm kiếm cục bộ được cân bằng, tương tự như trong 8 . Một sự khác biệt nữa giữa các thuật toán MDE và DE là trong thuật toán MDE hệ số đột biến F và xác suất lai ghép CR được tạo một cách ngẫu nhiên chứ không cố định như trong thuật toán DE. Mã code của các thuật toán DE và MDE được trình bày chi tiết trong 5 , 6 , 9 , 10 .
Thuật toán MDE đã được lập trình bằng ngôn ngữ Matlab bởi Hồ Hữu Vịnh với tên gọi: Matlab code "A fast and efficient differential evolution for truss optimization problems" 9 . Đây là một chương trình đơn giản nhưng rất hiệu quả với tốc độ hội tụ nhanh và cho lời giải chính xác khi áp dụng để giải bài toán ngược 1D phương pháp thăm dò MT với hàm mục tiêu là sai số căn quân phương (root mean square-rms) giữa điện trở suất và pha của số liệu quan sát và tính toán được cho bởi công thức sau:
Trong đó và là các giá trị điện trở suất và pha quan sát, và là các giá trị điện trở suất và pha tính từ mô hình, N là số chu kì (hoặc tần số).
KẾT QUẢ VÀ THẢO LUẬN
Thử nghiệm thuật toán MDE trên mô hình lí thuyết
Mục này trình bày việc sử dụng thuật toán MDE để giải bài toán ngược 1D đối với 2 mô hình lí thuyết tiêu biểu trong thăm dò điện. Trước tiên cần giải bài toán thuận một chiều để thu được các số liệu điện trở suất và pha của các mô hình. Trong phương pháp MT, việc giải bài toán thuận 1D đối với môi trường phân lớp ngang gồm N lớp có thể được thực hiện dễ dàng nhờ có sẵn các công thức tính do các tác giả nước ngoài công bố. Bài báo này sử dụng công thức của tác giả Patella 11 để tính toán bài toán thuận.
Thử nghiệm thuật toán MDE trên mô hình lí thuyết
Mục này trình bày việc sử dụng thuật toán MDE để giải bài toán ngược 1D đối với 2 mô hình lí thuyết tiêu biểu trong thăm dò điện. Trước tiên cần giải bài toán thuận một chiều để thu được các số liệu điện trở suất và pha của các mô hình. Trong phương pháp MT, việc giải bài toán thuận 1D đối với môi trường phân lớp ngang gồm N lớp có thể được thực hiện dễ dàng nhờ có sẵn các công thức tính do các tác giả nước ngoài công bố. Bài báo này sử dụng công thức của tác giả Patella 11 để tính toán bài toán thuận.
Thử nghiệm thuật toán MDE trên mô hình lí thuyết
Mục này trình bày việc sử dụng thuật toán MDE để giải bài toán ngược 1D đối với 2 mô hình lí thuyết tiêu biểu trong thăm dò điện. Trước tiên cần giải bài toán thuận một chiều để thu được các số liệu điện trở suất và pha của các mô hình. Trong phương pháp MT, việc giải bài toán thuận 1D đối với môi trường phân lớp ngang gồm N lớp có thể được thực hiện dễ dàng nhờ có sẵn các công thức tính do các tác giả nước ngoài công bố. Bài báo này sử dụng công thức của tác giả Patella 11 để tính toán bài toán thuận.
Áp dụng thuật toán MDE cho số liệu thực tế
Số liệu sử dụng để nghịch đảo là số liệu điện trở suất trung bình và pha trung bình của điểm đo MT C12, thuộc địa phận huyện Củ Chi, Thành phố Hồ Chí Minh. Tọa độ của điểm đo là 11°5,728’ vĩ độ Bắc và 106°28,961’ kinh độ Đông, ở độ cao 22 m so với mực nước biển. Đây là một điểm đo thuộc tuyến MT Củ Chi‒Bến Cát được thực hiện vào năm 2007. Tuyến đo này kéo dài khoảng 15,5 km, gồm 15 điểm đo, kí hiệu từ C01 đến C15. Kết quả xử lý của tuyến đo này đã được công bố trong một số công trình nghiên cứu trước đây 12 , 13 . Trong các công trình đó, điểm C12 đã được chọn để phân tích 1D nhằm nghiên cứu ảnh hưởng của dòng điện xích đạo đối với số liệu điện từ trong khu vực khảo sát. Kết quả cho thấy trong khu vực khảo sát, các số liệu điện từ có chu kì lớn hơn khoảng 5 s bị ảnh hưởng bởi dòng điện xích đạo. Do đó các số liệu đo đạc ở trong dải chu kì 0,0127‒4,8 s của điểm C12 đã được sử dụng để nghịch đảo nhằm thu được cấu trúc địa điện của khu vực Củ Chi và thu được kết quả là cấu trúc gồm 3 lớp: lớp điện trở suất thấp ở trên cùng, lớp điện trở suất cao ở giữa và dưới cùng là lớp điện trở suất thấp.
Trong bài báo này, số liệu sử dụng để phân tích 1D là các giá trị trung bình nhân của điện trở suất và pha của điểm C12:
trong đó , và , là điện trở suất và pha theo các hướng vuông góc với nhau. Các đường cong điện trở suất và pha trung bình của điểm C12 được biểu diễn trên Figure 7 .
Trước tiên, một mô hình địa điện ban đầu gồm 3 lớp đã được sử dụng để tìm kiếm lời giải tối ưu của bài toán. Phạm vi tìm kiếm đối với điện trở suất của các lớp từ trên xuống lần lượt trong các khoảng (1‒1000 Ω.m), (1‒10000 Ω.m) và (1‒1000 Ω.m). Bề dày của lớp trên cùng và lớp giữa lần lượt được cho tìm kiếm trong khoảng (100‒10000 m), (100‒20000 m), còn dưới cùng là nửa không gian đồng nhất. Chương trình nghịch đảo sử dụng thuật toán MDE đã được áp dụng với quần thể gồm 50 cá thể, sai số rms để dừng chạy chương trình là 10 -6 giống như đã áp dụng đối với các mô hình lí thuyết đã trình bày ở trên. Vì các khoảng tìm kiếm được đặt khá rộng nên việc tìm kiếm tốn nhiều thời gian, do đó số lần lặp tối đa được đặt bằng 10000. Máy tính sử dụng để chạy chương trình là máy tính bàn cài Window 11 với bộ vi xử lí i5-11400F, RAM 16 Gb. Chương trình chạy xong sau 930,81 s ( khoảng 15,5 phút ). Sai số rms theo số lần lặp được biểu diễn bằng đường màu đỏ trên Figure 8 . Có thể nhận thấy xu hướng chung là các sai số giảm dần theo số lần lặp, chứng tỏ sự hội tụ của quá trình lặp. Tại các bước lặp thứ 7140 đến 7379 sai số giảm nhanh theo số lần lặp, sau đó sai số có giá trị bằng 5.10 -5 và không giảm nữa. Bằng cách vẽ số liệu điện trở suất và pha tính từ kết quả nghịch đảo và số liệu quan sát trên cùng hình vẽ để đánh giá mức độ phù hợp của kết quả nghịch đảo, thì thấy rằng từ bước lặp thứ 7290 số liệu tính lại từ kết quả nghịch đảo đã gần như trùng với số liệu quan sát. Từ bước lặp 7379 trở đi các thông số mô hình thay đổi ít, không đáng kể. Các thông số mô hình tại bước lặp thứ 10000 (lần lặp cuối cùng) được trình bày ở Bảng 3, cho thấy lớp trên cùng dày 760 m có điện trở suất bằng 119 Ω.m, lớp giữa dày khoảng 10 km, có điện trở suất bằng 4353 Ω.m, và nửa không gian bên dưới có điện trở suất 41,3 Ω.m. Cấu trúc địa điện cùng với số liệu quan sát và số liệu tính toán từ cấu trúc này được trình bày trên Figure 9 , có thể thấy số liệu quan sát và số liệu tính toán hoàn toàn chồng khít với nhau.
Figure 7 . Các đường cong điện trở suất biểu kiến trung bình và pha trung bình của điểm đo C12
Figure 8 . Hàm mục tiêu (sai số rms) của cá thể tốt nhất trong mỗi lần lặp khi sử dụng các mô hình 3 lớp và 4 lớp để nghịch đảo
Tiếp theo, mô hình ban đầu gồm 4 lớp được sử dụng để nghịch đảo. Từ kết quả nghịch đảo dùng 3 lớp ở trên ta đã biết sơ bộ giá trị điện trở suất của cấu trúc bên dưới mặt đất nên phạm vi tìm kiếm đối với điện trở suất của các lớp lúc này được thu hẹp bớt để giảm thời gian tính toán. Cụ thể, phạm vi tìm kiếm điện trở suất của các lớp từ trên xuống lần lượt là các khoảng (100 - 150 Ω.m), (2000 - 6000 Ω.m), (10 - 100 Ω.m) và (10 - 100 Ω.m). Bề dày của lớp trên cùng và lớp giữa lần lượt được cho tìm kiếm trong các khoảng (500 - 1000 m), (8000 - 12000 m), (6000 - 10000 m) và dưới cùng là nửa không gian đồng nhất. Chương trình nghịch đảo dùng thuật toán MDE đã được áp dụng với quần thể gồm 50 cá thể, sai số rms để dừng chạy chương trình là 10 -6 như đã áp dụng ở trên. Số lần lặp được đặt bằng 6000. Chương trình tính toán kết thúc sau 571,72 s ( khoảng 9,5 phút). Sai số rms theo số lần lặp được biểu diễn bằng đường màu xanh trên Figure 8 . Tại các bước lặp từ 3293 đến 3439 sai số giảm nhanh theo số lần lặp, sau đó sai số có giá trị bằng 3.10 -5 và không giảm nữa. Bằng cách vẽ số liệu điện trở suất và pha tính từ kết quả nghịch đảo và số liệu quan sát trên cùng hình vẽ thì thấy rằng từ bước lặp khoảng 1600 trở đi thì số liệu tính lại từ kết quả nghịch đảo đã gần như trùng với số liệu quan sát, tuy nhiên các thông số mô hình kết quả vẫn còn thay đổi. Từ bước lặp 3461 cho đến lần lặp cuối cùng (bằng 6000) thì các thông số mô hình thay đổi rất ít, không đáng kể. Các thông số mô hình tại bước lặp thứ 6000 được trình bày ở ảng 4. Có thể thấy các thông số của cấu trúc 4 lớp này thực chất là cấu trúc 3 lớp, vì 2 lớp cuối có giá trị điện trở suất gần như bằng nhau, và mô hình 3 lớp này gần như trùng với mô hình 3 lớp đã thu được ở trên, ngoại trừ lớp điện trở suất cao ở giữa của 2 mô hình có chênh lệch nhau chút ít (4343 Ω.m so với 4353 Ω.m, tương đương sai số tương đối 0,2%). Cấu trúc địa điện cùng với số liệu quan sát và số liệu tính toán từ cấu trúc này được trình bày trên Figure 9 . Trên hình vẽ này có thể thấy các mô hình 3 lớp và 4 lớp gần như trùng nhau, còn số liệu tính từ mô hình hoàn toàn khớp với số liệu quan sát.
Kết quả thu được về cấu trúc điện trở suất của điểm C12 tại khu vực Củ Chi thu được bằng thuật toán MDE trong bài báo này phù hợp với các kết quả đã công bố trước đây 1 về số lớp cũng như về điện trở suất và bề dày của từng lớp. Do đó chứng tỏ hiệu quả của thuật toán MDE trong việc nghịch đảo 1D số liệu MT.
KẾT LUẬN
Thuật toán MDE là một thuật toán tối ưu toàn cục đã áp dụng thành công trong các bài toán cơ kỹ thuật như tạo dáng đi cho robot 8 , 15 , tối ưu dầm trong xây dựng 10 , điều khiển robot 6 , 7 . Việc thử nghiệm thuật toán này cho bài toán nghịch đảo một chiều phương pháp MT đối với số liệu mô hình và số liệu thực tế cho thấy kết quả tốt. Do đó hoàn toàn có thể áp dụng thuật toán này cho việc giải bài toán ngược MT trong các trường hợp phức tạp hơn, chẳng hạn trường hợp môi trường bên dưới gồm nhiều lớp, hoặc môi trường có tính chất 2D... Đây sẽ là những chủ đề mà chúng tôi hướng tới trong các công trình nghiên cứu sau.
Figure 9 . Kết quả thu được bằng cách sử dụng 3 lớp và 4 lớp để nghịch đảo (phần bên trái) và các đường cong điện trở suất và pha quan sát và tính lại từ kết quả nghịch đảo của điểm C12.
LỜI CẢM ƠN
Tác giả chân thành cảm ơn TS. Trần Thiện Huân về những thảo luận bổ ích trong quá trình dùng thuật toán vi sai cải tiến để tính toán và ThS. Nguyễn Văn Thuận đã giúp chỉnh sửa hoàn thiện bài báo.
DANH MỤC CÁC TỪ VIẾT TẮT
MT: Magnetotellurics – Thăm dò từ tellua (Trường điện từ)
MDE: Modified Differential Evolution – Thuật toán tiến hóa vi sai cải tiến
GA: Gene Algorithm – Thuật toán di truyền
DE : Differential Evolution – Thuật toán tiến hóa vi sai
1D: One-Dimensional – 1 chiều
2D: Two-Dimensional – 2 chiều
XUNG ĐỘT LỢI ÍCH
Tác giả tuyên bố rằng không có xung đột lợi ích.
ĐÓNG GÓP CỦA TÁC GIẢ
Tác giả Lưu Việt Hùng xây dựng các mô hình lý thuyết, xử lý số liệu thực tế, vẽ hình và viết bản thảo.
References
- Constable SC, Parker RL, Constable CG. Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics. 1987;52(3):289-300. . ;:. Google Scholar
- Lương PT. Giải bài toán ngược trọng lực bằng thuật toán memetic-áp dụng phân tích tài liệu trọng lực vùng Đồng bằng sông Cửu Long. 2018. . ;:. Google Scholar
- Holland JH. Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence. MIT press Cambridge; 1992. . ;:. Google Scholar
- Grandis H, Maulana Y, editors. Particle swarm optimization (PSO) for magnetotelluric (MT) 1D inversion modeling. IOP conference series: Earth and Environmental Science; 2017: IOP Publishing. . ;:. Google Scholar
- Storn R, Price K. Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces. Journal of global optimization. 1997;11:341-59. . ;:. Google Scholar
- Son NN, Van Kien C, Anh HPH. A novel adaptive feed-forward-PID controller of a SCARA parallel robot using pneumatic artificial muscle actuator based on neural network and modified differential evolution algorithm. Robotics and Autonomous Systems. 2017;96:65-80. . ;:. Google Scholar
- Son NN, Anh HPH, Chau TD. Adaptive neural model optimized by modified differential evolution for identifying 5-DOF robot manipulator dynamic system. Soft Computing. 2018;22(3):979-88. . ;:. Google Scholar
- Huân TT, Ánh HPH. Tối ưu hóa dáng đi ổn định cho robot dạng người kích thước nhỏ sử dụng thuật toán tiến hóa vi phân cải tiến. Chuyên san đo lường, Điều khiển và Tự động hóa; 2018. Contract No.: 63-74. . ;:. Google Scholar
- Ho V. A fast and effective differential evolution algorithm MATLAB Central File Exchange: MATLAB Central File Exchange; 2023 [cited 2023. . ;:. Google Scholar
- Ho-Huu V, Vo-Duy T, Luu-Van T, Le-Anh L, Nguyen-Thoi T. Optimal design of truss structures with frequency constraints using improved differential evolution algorithm based on an adaptive mutation scheme. Automation in Construction. 2016;68:81-94. . ;:. Google Scholar
- Patella D. Interpretation of magnetotelluric resistivity and phase soundings over horizontal layers. Geophysics. 1976;41(1):96-105. . ;:. Google Scholar
- Hung LV, Menvielle M, Son VT, Thang NC, Marquis G, Trieu CD. Studying the deep structure in the Saigon river fault area by the magnetotelluric method. Vietnam Journal of Earth Sciences. 2014;36(3):233-40. . ;:. Google Scholar
- Luu HV, Vu TT, Lee SK, Ngo MS, Tran DT, Dinh TQ, et al. Using the magnetotelluric method to study the deep structure of the Sai Gon river fault in the vicinity of Ho Chi Minh City. VNUHCM Journal of Natural Sciences. 2022;6(2):2103-15. . ;:. Google Scholar
- Van Pham N, Boyer D, Thoa Nguyen TK, Van Nguyen G. Deep Ground‐Water Investigation by Combined VES/MTS Methods Near Ho Chi Minh city, Viet Nam. Groundwater. 1994;32(4):675-82. . ;:. Google Scholar
- Huan TT, Van Kien C, Anh HPH, Nam NT. Adaptive gait generation for humanoid robot using evolutionary neural model optimized with modified differential evolution technique. Neurocomputing. 2018;320:112-20. . ;:. Google Scholar
