One–dimensional magnetotelluric inversion by modified differential evolution algorithm
- Ho Chi Minh City University of Technology and Education, Vietnam
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ẫn2). 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ả Holland3. 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 19975. 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ệp6, 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ư trong8. 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 trong5, 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ả Patella11 để tính toán bài toán thuận.
Mô hình 1
Mô hình cấu trúc điện trở suất gồm 3 lớp có dạng K (K-type) trong thăm dò điện: lớp thứ nhất có điện trở suất 100 Ω.m với bề dày 2 km; lớp thứ 2 có điện trở suất 5000 Ω.m với bề dày 8 km và lớp thứ 3 là nửa không gian bên dưới có điện trở suất 50 Ω.m (Bảng 1). Mô hình này được xây dựng dựa theo cấu trúc điện trở suất của khu vực Củ Chi theo kết quả đo MT12, 13. Mô hình cùng với số liệu điện trở suất và pha tính tại 50 chu kỳ từ 10 (s) đến 10 (s) phân bố đều trên thang logarith được trình bày trong Figure 1.
Mô hình cấu trúc 3 lớp dạng K (phần bên trái) và các đường cong điện trở suất và pha tính từ mô hình.
Thuật toán MDE được áp dụng với quần thể gồm 50 cá thể, số lần lặp là 500, sai số căn quân phương (root mean square-rms) để dừng chạy chương trình là 10. Giá trị điện trở suất của các lớp được cho tìm kiếm là 1‒10000 Ω.m và bề dày các lớp được cho tìm kiếm là 10‒10 km. Chương trình được cho chạy tự động với số lớp được đặt bằng 3, 4 và 5. Kết quả điện trở suất và bề dày của các lớp được trình bày trong Bảng 1. Các kết quả này cùng với các giá trị điện suất và pha tính toán từ mô hình được biểu diễn trên Figure 2. Kết quả cho thấy trong cả 3 trường hợp (số lớp bằng 3, 4 và 5), các giá trị điện trở suất và pha tính lại từ kết quả nghịch đảo đều gần như trùng với giá trị quan sát.
Các thông số của mô hình lý thuyết dạng K và kết quả thu được bằng cách sử dụng mô hình ban đầu gồm 3 lớp, 4 lớp và 5 lớp
|
Mô hình 1 |
Kết quả nghịch đảo | ||||||
|
3 lớp |
4 lớp |
5 lớp | |||||
|
Điện trở suất (Ω.m) |
Bề dày (m) |
Điện trở suất (Ω.m) |
Bề dày (m) |
Điện trở suất (Ω.m) |
Bề dày (m) |
Điện trở suất (Ω.m) |
Bề dày (m) |
|
100 |
2000 |
100,00 |
2000,56 |
100,00 |
2001,12 |
100,00 |
1999,15 |
|
5000 |
8000 |
5040,05 |
7998,91 |
5104,35 |
5522,33 |
4810,74 |
4694,75 |
|
50 |
50,00 |
4790,92 |
2477,00 |
8521,54 |
851,58 | ||
|
50,00 |
4684,40 |
2453,87 | |||||
|
50,00 | |||||||
Mô hình cấu trúc 3 lớp dạng K và các kết quả thu được bằng cách sử dụng 3 lớp, 4 lớp và 5 lớp để nghịch đảo (phần bên trái) và các đường cong điện trở suất và pha tính từ mô hình và tính từ kết quả nghịch đảo.
Hàm mục tiêu (hay sai số rms giữa các giá trị điện trở suất và pha tính lại từ mô hình nghịch đảo và giá trị quan sát) theo số lần lặp được trình bày trên Figure 3. Kết quả cho thấy đối với trường hợp mô hình nghịch đảo 3 lớp (số lớp chính xác như mô hình lý thuyết), kết quả nghịch đảo hội tụ đều đặn và đến bước lặp khoảng 170 thì sai số rms gần như không đổi và đến bước lặp 248 thì quá trình lặp kết thúc. Lớp trên cùng và nửa không gian bên dưới cho kết quả chính xác như giá trị mô hình, còn lớp giữa kết quả thu được là 5040,05 m so với lý thuyết là 5000 m vì vậy sai số tương đối chỉ 0,8%. Như vậy trong trường hợp này, khi số lớp đưa vào nghịch đảo bằng với số lớp mô hình, thì kết quả nghịch đảo có thể xem là chính xác.
Trường hợp nghịch đảo sử dụng mô hình ban đầu gồm 4 lớp và 5 lớp, kết quả cho thấy lớp trên cùng và nửa không gian bên dưới hoàn toàn trùng khớp với giá trị của mô hình lý thuyết. Đối với trường hợp nghịch đảo 4 lớp, thì 2 lớp ở giữa thu được có điện trở suất và bề dày gần tương đương với lớp ở giữa của mô hình, và khi biểu diễn lên hình vẽ (Figure 2, phần bên trái) kết quả cho thấy chúng gần như trùng khít với mô hình lý thuyết. Trường hợp sử dụng mô hình 5 lớp để nghịch đảo, khi biểu diễn kết quả trên Figure 2, nhận thấy hai lớp trên cùng và hai lớp dưới cùng gần như trùng với mô hình lý thuyết còn lớp ở giữa có giá trị lớn hơn 2 lớp sát nó và bị nhô ra. Quan sát các đường cong sai số trên hình 3, có thể thấy rằng đối với các trường hợp 4 lớp và 5 lớp thì quá trình hội tụ chậm hơn.
Hàm mục tiêu (sai số) của các thể tốt nhất trong mỗi lần lặp khi sử dụng 3 lớp, 4 lớp và 5 lớp để nghịch đảo
Kết quả cho thấy đối với mô hình lý thuyết, khi số lớp đưa vào nghịch đảo đúng bằng số lớp của mô hình thì kết quả thu được gần như chính xác, còn khi số lớp đưa vào nghịch đảo càng khác so với số lớp thực tế thì kết quả bị sai lệch nhiều.
Mô hình 2
Mô hình này là cấu trúc điện trở suất có dạng H (H-type) trong thăm dò điện, gồm 3 lớp: lớp thứ nhất có điện trở suất 300 Ω.m với bề dày 200 m, lớp thứ 2 có điện trở suất 100 Ω.m với bề dày 300 m và lớp thứ 3 là nửa không gian bên dưới có điện trở suất 500 Ω.m ( ảng 2). Mô hình này được xây dựng dựa theo cấu trúc điện trở suất khu vực Củ Chi thu được bằng phương pháp MT đo ở dải tần số cao của nhóm tác giả Phạm Văn Ngọc và nnk14. Mô hình cùng với số liệu điện trở suất và pha tính tại 50 chu kỳ từ 10 (s) đến 10 (s) phân bố đều trên thang logarith được trình bày trên Figure 4.
Mô hình cấu trúc 3 lớp dạng H (phần bên trái) và các đường cong điện trở suất và pha tính từ mô hình.
Thuật toán MDE được áp dụng để nghịch đảo số liệu quan sát với quần thể gồm 50 cá thể, số lần lặp là 500, sai số căn quân phương để dừng chạy chương trình là 10 tương tự như đã áp dụng đối với Mô hình 1. Giá trị điện trở suất của các lớp bây giờ được cho tìm kiếm là 1‒1000 Ω.m, và bề dày các lớp được cho tìm kiếm là 10‒1000 m. Chương trình đã được chạy tự động với số lớp được đặt lần lượt bằng 3, 4 và 5. Các kết quả điện trở suất và bề dày các lớp được trình bày ở ảng 2. Các cấu trúc thu được từ nghịch đảo cùng với các giá trị điện trở suất và pha tính lại từ chúng được biểu diễn trên Figure 5. Kết quả cho thấy trong cả 3 trường hợp (số lớp bằng 3, 4 và 5), các đường cong điện trở suất và pha tính lại từ kết quả nghịch đảo gần như trùng với giá trị quan sát từ mô hình. Kết quả trong Bảng 2 cho thấy trong cả 3 trường hợp (mô hình nghịch đảo 3, 4, và 5 lớp), điện trở suất và bề dày của lớp trên cùng cũng như điện trở suất của nửa không gian bên dưới được xác định một cách khá chính xác. Đối với mô hình 3 lớp thì các thông số của lớp giữa được xác định chính xác, trong khi đối với các mô hình có số lớp lớn hơn thì các lớp ở giữa có những sai lệch về điện trở suất và độ sâu, tuy nhiên về tổng thể (Hình 5) thì những sai lệch đó cũng không lớn mà về cơ bản vẫn là cấu trúc 3 lớp.
Sai số căn quân phương giữa các giá trị điện trở suất và pha tính lại từ mô hình nghịch đảo và giá trị quan sát được trình bày trên Figure 6. Đối với mô hình nghịch đảo 3 lớp, 4 lớp và 5 lớp, quá trình nghịch đảo kết thúc tương ứng ở bước lặp 216, 386 và 500, cho thấy khi mô hình ban đầu càng phức tạp thì quá trình lặp càng lâu hội tụ.
Hình vẽ biểu diễn kết quả thu được (phần bên trái Figure 5) cho thấy các mô hình nghịch đảo thu được có dạng gần giống với mô hình lý thuyết. Khi số lớp đưa vào nghịch đảo đúng bằng số lớp của mô hình thì kết quả thu được gần như chính xác. Khi số lớp nghịch đảo bằng 4 hoặc 5 thì kết quả thu được có sự sai lệch với mô hình lý thuyết. Tuy nhiên những khác biệt đó không lớn và có thể chấp nhận được trong lĩnh vực địa vật lý thăm dò.
Như vậy, việc thử nghiệm trên các mô hình lý thuyết có dạng đơn giản và tiêu biểu của thăm dò điện trở suất cho thấy thuật toán MDE cho kết quả tốt, có thể áp dụng để nghịch đảo số liệu thực tế.
Các thông số của mô hình lý thuyết dạng H và kết quả thu được bằng cách sử dụng mô hình ban đầu gồm 3 lớp, 4 lớp và 5 lớp.
|
Mô hình 2 |
Kết quả nghịch đảo | ||||||
|
3 lớp |
4 lớp |
5 lớp | |||||
|
Điện trở suất (Ω.m) |
Bề dày (m) |
Điện trở suất (Ω.m) |
Bề dày (m) |
Điện trở suất (Ω.m) |
Bề dày (m) |
Điện trở suất (Ω.m) |
Bề dày (m) |
|
300 |
200 |
302,30 |
197,72 |
302,28 |
197,76 |
301,71 |
198,34 |
|
100 |
300 |
100,78 |
304,11 |
100,77 |
302,75 |
100,55 |
300,45 |
|
500 |
500,10 |
331,17 |
10,00 |
450,74 |
157,82 | ||
|
500,10 |
585,68 |
50,38 | |||||
|
499,97 | |||||||
Mô hình cấu trúc 3 lớp dạng H và các kết quả thu được bằng cách sử dụng 3 lớp, 4 lớp và 5 lớp để nghịch đảo (phần bên trái) và các đường cong điện trở suất và pha tính từ mô hình và tính từ kết quả nghịch đảo.
Hàm mục tiêu (sai số) của các thể tốt nhất trong mỗi lần lặp khi sử dụng 3 lớp, 4 lớp và 5 lớp để nghịch đảo
Á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 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 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.
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
Các khoảng tìm kiếm của các thông số mô hình 3 lớp ban đầu và kết quả nghịch đảo thu được
|
Lớp |
Khoảng giá trị của mô hình3 lớp ban đầu |
Kết quả nghịch đảo 3 lớp | ||
|
Điện trở suất (Ω.m) |
Bề dày (km) |
Điện trở suất (Ω.m) |
Bề dày (km) | |
|
1 |
1‒1000 |
0,1‒10 |
119,04 |
0,760 |
|
2 |
1‒100000 |
0,1‒20 |
4352,99 |
9,999 |
|
3 |
1‒1000 |
- |
41,30 |
- |
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 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 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 robot8, 15, tối ưu dầm trong xây dựng10, đ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.
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.
Các khoảng tìm kiếm của các thông số mô hình 4 lớp ban đầu và kết quả mô hình thu được từ nghịch đảo
|
Lớp |
Khoảng giá trị của mô hình 4 lớp ban đầu |
Kết quả nghịch đảo 4 lớp | ||
|
Điện trở suất (Ω.m) |
Bề dày (km) |
Điện trở suất (Ω.m) |
Bề dày (km) | |
|
1 |
1‒1000 |
0,1‒10 |
118,97 |
0,759 |
|
2 |
1‒100000 |
0,1‒20 |
4343,90 |
10,000 |
|
3 |
1‒1000 |
0,1‒10 |
41,29 |
0,794 |
|
4 |
1‒1000 |
- |
41,28 |
- |
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.