ĐẠI HỌC QUỐC GIA HÀ NỘI
TRƢỜNG ĐẠI HỌC KHOA HỌC TỰ NHIÊN
---------------------
NGUYỄN KIM DŨNG
GIẢI BÀI TOÁN NGƢỢC 3D XÁC ĐỊNH SỰ PHÂN BỐ MẬT ĐỘ CỦA
ĐÁ MÓNG THEO TÀI LIỆU DỊ THƢỜNG TRỌNG LỰC
LUẬN VĂN THẠC SĨ KHOA HỌC
Hà Nội - 2012
ĐẠI HỌC QUỐC GIA HÀ NỘI
TRƢỜNG ĐẠI HỌC KHOA HỌC TỰ NHIÊN
---------------------
NGUYỄN KIM DŨNG
GIẢI BÀI TOÁN NGƢỢC 3D XÁC ĐỊNH SỰ PHÂN BỐ MẬT ĐỘ CỦA
ĐÁ MÓNG THEO TÀI LIỆU DỊ THƢỜNG TRỌNG LỰC
Chuyên ngành: Vật lý địa cầu
Mã số: 60.44.15
LUẬN VĂN THẠC SĨ KHOA HỌC
NGƢỜI HƢỚNG DẪN KHOA HỌC:
PGS.TS. Đỗ Đức Thanh
Hà Nội - 2012
MỞ ĐẦU
Chƣơng 1 – CÁC PHƢƠNG PHÁP XÁC ĐỊNH DỊ THƢỜNG TRỌNG LỰC ĐỐI VỚI CÁC VẬT THỂ CÓ DẠNG HÌNH HỌC ĐỀU ĐẶN
1.1. Những khái niệm cơ bản ..................................................................... 2 1.2. Các biểu thức tích phân tổng quát về đạo hàm của thế trọng lực ...... 3 1.3. Bài toán thuận cho những vật thể có dạng hình học ........................... 6 1.3.1. Hình cẩu hoặc điểm vật chất .................................................... 6 1.3.2. Thanh vật chất nằm ngang, hình trụ tròn nằm ngang ............... 8 1.3.3. Nửa mặt phẳng vật chất nằm ngang .......................................... 9 1.3.4. Hình hộp vuông góc ............................................................... 11 1.3.5. Lăng trụ thẳng đứng ................................................................. 12 1.3.6. Bậc thẳng đứng ....................................................................... 12 1.3.7. Bậc nghiêng ........................................................................... 15
Chƣơng 2 - PHƢƠNG PHÁP XÁC ĐỊNH DỊ THƢỜNG TRỌNG LỰC CỦA
CÁC RANH GIỚI 2D VÀ 3D
2.1. Phƣơng pháp xác định dị thƣờng trọng lực của ranh giới trầm tích 2D trong miền không gian .............................................................. 18 2.1.1. Xác định dị thƣờng trọng lực của ranh giới trầm tích trên cơ sở phân chia nó thành các đa giác có tiết diện ngang bất kỳ ...... 18 2.1.2. Trƣờng hợp mật độ dƣ thay đổi tuyến tính theo độ sâu ........... 19 2.1.3. Trƣờng hợp mật độ dƣ thay đổi theo quy luật hàm mũ theo chiều Sâu ................................................................................................. 20 2.1.4. Mật độ dƣ thay đổi theo dạng hàm hypepol ............................ 21 2.1.5. Mật độ dƣ thay đổi theo dạng hàm parabolic ........................... 23 2.1.6. Xác định dị thƣờng trọng lực của ranh giới trầm tích trên cơ sở phân ........................................................................................ 24
2.2. Phƣơng pháp xác định dị thƣờng trọng lực của ranh giới trầm tích
3D trong miền không gian. ............................................................... 25 2.2.1. Cơ sở lý thuyết .......................................................................... 26 2.2.2. Xây dựng chƣơng trình giải bài toán ngƣợc 3D và tính toán thử nghiệm trên mô hình ................................................................. 29 2.3. Các phƣơng pháp xác định độ sâu của bể trầm tích 3D trong miền tần số ............................................................................................. 31 2.3.1. Nâng cao độ chính xác của việc tính dị thƣờng trọng lực trong
miền tần số số bằng phƣơng pháp "trƣợt mẫu" (Shift-sampling) ... 32
2.3.2. Xác định dị thƣờng trọng lực của ranh giới 3D trọng lực trong miền tần số .................................................................................................... 33
2.3.3. Xây dựng chƣơng trình giải bài toán ngƣợc và tính toán
thử nghiệm trên các mô hình ................................................................ 36
Chƣơng 3 - MÔ HÌNH HÓA VIỆC GIẢI BẢI TOÁN NGƢỢC XÁC ĐỊNH SỰ PHÂN BỐ MẬT ĐỘ MÓNG KẾT TINH 3.1. Các phƣơng pháp giải bài toán ngƣợc ............................................. 38 3.1.1. Xác định sự phân bố mật độ đá móng theo phƣơng pháp lựa chọn .......................................................................................... 38 3.1.2. Xác định sự phân bố mật độ đá móng theo phƣơng pháp trực tiếp .......................................................................................... 41 3.2. Thuật toán và sơ đồ khối .................................................................. 42 3.2.1. Thuật toán ................................................................................. 42 3.2.2. Sơ đồ khối ............................................................................... 43 3.3. Tính toán thử nghiệm trên mô hình .................................................. 44 3.3.1. Mô hình 1 .............................................................................. 45 3.3.2. Mô hình 2 ............................................................................... 50 3.3.3. Mô hình 3 ................................................................................. 54 KẾT LUẬN ........................................................................................... 58 TÀI LIỆU THAM KHẢO .................................................................... 59
Danh mục các hình vẽ
Hình 1.1 Xác định thế các đạo hàm của một chất điểm ................................ 4 Hình 1.2 Xác định thế và đạo hàm của các vật thể hai chiều ..................... 6 Hình 1.3 Xác định thế và các đạo hàm của vật thể hình cầu ...................... 7 Hình 1.4 Trƣờng trọng lực của hình cầu ....................................................... 8 Hình 1.5 Tƣờng trọng lực của hình trụ tròn nằm ngang ............................... 9 Hình 1.6 Xác định thế và các đạo hàm của nửa mặt phẳng vật chất nằm ngang ......................................................................................... 10 Hình 1.7 Trƣờng trọng lực của nửa mặt phẳng vật chất nằm ngang ......... 11 Hình 1.8 Bậc thẳng đứng ........................................................................... 13 Hình 1.9 Trƣờng trọng lực trên bậc thẳng đứng ........................................ 15 Hình 1.10 Bậc nghiêng ................................................................................ 15 Hình 2.1 Vật thể hai chiều có tiết diện ngang bất kỳ ................................ 18 Hình 2.2 Xấp xỉ vật thể có tiết diện ngang ................................................ 20 Hình 2.3 Việc phân chia mỗi cạnh đa giác ................................................. 20 Hình 2.4 Ranh giới phân chia mật độ và xấp xỉ nó bằng các lăng trụ ........ 24 Hình 2.5 Mô hình lăng trụ 3 chiều ............................................................. 26 Hình 2.6 Mô hình khối đa diện .................................................................. 28 Hình 2.7 Kết quả xác định độ sâu bể trầm tích trong miền không gian ..... 31 Hình 2.8 Kết quả xác định độ sâu bể trầm tích trong miền tần số ............. 37
Hình 3.1. Sơ đồ khối giải bài toán ngƣợc 3D xác định sự
phân bố mật độ của đá móng ..................................................... 43 Hình 3.2 Mô hình sự phân bố mật độ dƣ trong đá móng ........................... 45 Hình 3.3 Địa hình các ranh giới và các thành phần trƣờng tƣơng ứng ...... 46 Hình 3.4 Trƣờng quan sát .......................................................................... 46 Hình 3.5 Tƣơng quan giữa trƣờng phông bậc 3 và các mức nâng trƣờng .. 47 Hình 3.6 Trƣờng phông khu vực ............................................................... 47 Hình 3.7 Trƣờng móng dƣ ........................................................................ 47 Hình 3.8 Trƣờng móng dƣ ở lần lặp cuối .................................................. 48 Hình 3.9 Sai số giữa trƣờng móng dƣ và trƣờng móng dƣ ở lần lặp cuối .. 48 Hình 3.10 Tốc độ hội tụ .............................................................................. 48 Hình 3.11 Kết quả tính toán sự phân bố mật độ dƣ trong đá móng ........... 49 Hình 3.12 Sai số giữa sự phân bố mật độ theo mô hình và tính toán ........ 49 Hình 3.13 Địa hình các ranh giới và các thành phần trƣờng tƣơng ứng .... 50 Hình 3.14 Trƣờng quan sát ......................................................................... 50 Hình 3.15 Tƣơng quan giữa trƣờng phông bậc 3 và các mức nâng trƣờng .. 51 Hình 3.16 Trƣờng phông khu vực ............................................................. 51 Hình 3.17 Trƣờng móng dƣ ....................................................................... 51 Hình 3.18 Trƣờng móng dƣ ở lần lặp cuối ................................................ 52
Hình 3.19 Sai số giữa trƣờng móng dƣ và trƣờng móng dƣ ở lần lặp cuối 52 Hình 3.20 Tốc độ hội tụ .............................................................................. 52 Hình 3.21 Kết quả tính toán sự phân bố mật độ dƣ trong đá móng ........... 53 Hình 3.22 Sai số giữa sự phân bố mật độ theo mô hình và tính toán ......... 53 Hình 3.23 Địa hình các ranh giới và các thành phần trƣờng tƣơng ứng .... 54 Hình 3.24 Trƣờng quan sát ........................................................................ 54 Hình 3.25 Tƣơng quan giữa trƣờng phông bậc 3 và các mức nâng trƣờng .. 55 Hình 3.26 Trƣờng phông khu vực ............................................................... 55 Hình 3.27 Trƣờng móng dƣ ...................................................................... 55 Hình 3.28 Trƣờng móng dƣ ở lần lặp cuối ................................................ 56 Hình 3.29 Sai số giữa trƣờng móng dƣ và trƣờng móng dƣ ở lần lặp cuối 56 Hình 3.30 Tốc độ hội tụ .............................................................................. 56 Hình 3.31 Kết quả tính toán sự phân bố mật độ dƣ trong đá móng ............ 57 Hình 3.32 Sai số giữa sự phân bố mật độ theo mô hình và tính toán ........ 57
MỞ ĐẦU
Trƣớc đây, khi nghiên cứu cấu trúc địa chất các nhà địa vật lý chủ yếu tập
trung nghiên cứu trong việc xác định độ sâu, hình dạng của các mặt ranh giới nhƣ:
Moho, conrat, đáy trầm tích Kainozoi với giả thiết sự phân bố mật độ trong trầm
tích và mật độ đá móng là không đổi. Trong những năm gần đây, với việc tìm thấy
dầu trong đá móng, việc nghiên cứu sự bất đồng nhất của mật độ cũng nhƣ sự phân
bố mật độ của đá móng đã hấp dẫn sự quan tâm của nhiều nhà nghiên cứu địa vật lý
trong nƣớc. Tuy nhiên, cho tới nay, sự phân bố mật độ của đá móng mới chỉ dừng
lại trong việc phân tích định tính hoặc cũng chỉ đƣợc xác định bằng phƣơng pháp
tƣơng quan nên độ chính xác vẫn còn nhiều hạn chế. Để góp phần vào việc nghiên
cứu này, trong phạm vi của bản luận văn, chúng tôi tiến hành nghiên cứu kết hợp tổ
hợp phƣơng pháp bóc lớp dị thƣờng với việc giải bài toán ngƣợc 3D theo phƣơng
pháp bình phƣơng tối thiểu Marquart, xây dựng thuật toán và chƣơng trình máy tính
xác định sự phân bố mật độ trong đá móng theo tài liệu dị thƣờng trọng lực. Thuật
toán và chƣơng trình xây dựng đƣợc tính toán thử nghiệm trên các mô hình 3D
nhằm nghiên cứu khả năng áp dụng của phƣơng pháp.
Khóa luận này đƣợc chia làm ba chƣơng sau:
- Chương 1: Các phƣơng pháp xác định dị thƣờng trọng lực đối với các vật
thể có dạng hình học đều đặn.
- Chương 2: Phƣơng pháp xác định dị thƣờng trọng lực của các ranh giới 2D
và 3D.
- Chương 3: Mô hình hóa việc giải bài toán ngƣợc xác định sự phân bố mật
độ móng kết tinh
1
Chƣơng 1
CÁC PHƢƠNG PHÁP XÁC ĐỊNH DỊ THƢỜNG TRỌNG LỰC
ĐỐI VỚI CÁC VẬT THỂ CÓ DẠNG HÌNH HỌC ĐỀU ĐẶN
1.1 Những khái niệm cơ bản.
Sau khi tu chỉnh số liệu đo đạc bằng máy trọng lực, ta thành lập bản đồ hoặc
đồ thị đạo hàm bậc nhất, hoặc bậc hai của thế trọng lực (dị thƣờng trọng lực nói
chung). Giải thích địa chất dị thƣờng trọng lực bao gồm phân tích các quy luật phân
bố của nó trên mặt đất (hoặc gần mặt đất) và mối liên hệ để nó giải quyết các nhiệm
vụ khác.
Giải thích địa chất dị thƣờng trọng lực đƣợc hình thành nhƣ sau: Dựa vào số
liệu trọng lực đo đạc và số liệu địa chất, địa vật lý sẵn có, vào kinh nghiệm giải
thích trọng lực tại các vùng tƣơng đƣơng, ta có thể đƣa ra những kết luận địa chất
về vùng cho trƣớc tƣơng ứng với nhiệm vụ địa chất đề ra.
Nhiệm vụ giải thích dị thƣờng trọng lực đƣợc phân loại dƣới hai hình thức:
Phân tích định tính và định lƣợng.
Khi giải thích định tính cần xác định:
- Các yếu tố địa chất chắc chắn ảnh hƣởng lên trƣờng trọng lực cũng nhƣ các
trƣờng vật lý khác (nếu nhƣ các phƣơng pháp Địa vật lý khác cũng đƣợc áp
dụng).
- Vị trí của yếu tố địa chất hoặc vật quặng.
- Vùng hoặc khu vực cần phải tiến hành nghiên cứu tỉ mỉ hơn.
- Điểm hoặc vùng nhỏ tại đó tại đó có thể đặt đƣợc các lỗ khoan hoặc đào
hầm lò.
- Khả năng và điều kiện để phân tích định lƣợng.
Trong trƣờng hợp tổng quát, có bốn yếu tố địa chất chính gây nên dị thƣờng
trọng lực:
- Cấu tạo các lớp trầm tích.
- Địa hình mặt nền kết tinh.
2
- Cấu tạo bên trong của nền kết tinh.
- Cấu tạo sâu vỏ Trái đất.
Khi minh giải định tính ta tiến hành mô tả một cách hệ thống các vùng dị
thƣờng và dị thƣờng riêng biệt, chỉ rõ bản chất địa chất dị thƣờng đƣợc mô tả với
xác suất lớn nhất, đƣa ra những đề nghị về việc tiến hành những nghiên cứu tiếp
theo.
Công tác phân tích định lƣợng đƣợc tiến hành khi:
- Có lƣợng thông tin đầy đủ hoặc tƣơng đối đầy đủ về địa chất của vùng, do
đó có khả năng hình thành mẫu vật lý của môi trƣờng địa chất dùng để
phân tích dị thƣờng trọng lực.
- Tác dụng của một trong những yếu tố địa chất gây nên dị thƣờng trội hơn.
Để đảm bảo yêu cầu này, trong thực tế ngƣời ta sử dụng các phƣơng pháp
biến đổi trƣờng.
- Các yếu tố địa chất trong vùng tƣơng đối ổn định có thể sử dụng một hoặc
tổ hợp phƣơng pháp phân tích chung.
- Các số liệu đo đạc có độ chi tiết và chính xác cao.
- Cơ sở lý thuyết phân tích tốt.
1.2. Các biểu thức tích phân tổng quát về đạo hàm của thế trọng lực
Để thuận tiện cho việc tính toán sau này ngƣời ta viết lại các biểu thức tích
phân tổng quát của thế hấp dẫn và các đạo hàm của chúng khi giải các bài toán
thuận và nghịch.
Thế V tại điểm với tọa độ x1, y1, z1 đƣợc biểu diễn bằng công thức:
(1.1)
trong đó:
3
O
y
A(x1,y1,z)
r
dm(x,y,z)
x
B
z
Hình 1.1: Xác định thế các đạo hàm của một chất điểm
Từ đó ta có:
(1.2)
(1.3)
(1.4)
(1.5)
Trong đó:
K : là hệ số hấp dẫn
Vxz: Đạo hàm của thế trọng lực theo phƣơng nằm ngang x
Vyz: Đạo hàm của thế trọng lực theo phƣơng nằm ngang y
g: Đạo hàm của thế trọng lực theo phƣơng thẳng đứng z
Đặt gốc tọa độ tại điểm quan sát A, tức là trong công thức đặt x1=y1=z1=0 thì
ta có:
4
(1.6) (1.7)
(1.8) (1.9)
(1.10)
Trong thực tế thƣờng gặp các vật thể có dạng kéo dài một hƣớng. Với độ
chính xác khá đủ có thể xem các vật thể đó là các vật thể hai chiều. Việc giải bài toán
thuận và nghịch đối với bài toán hai chiều đơn giản hơn nhiều so với bài toán ba
chiều. Để chuyển từ bài toán ba chiều về bài toán hai chiều thì trong công thức cho
vật thể ba chiều ở trên cần cho một biến chạy từ - +.
Để làm ví dụ, chúng ta xét trƣờng hợp Vz( g). Từ công thức ba chiều (1.2)
ta có:
Cho biến y chạy từ -+, lúc đó ta có:
(1.11)
Trong biểu thức (1.11), biến số y chạy từ - + còn các biến số (x,y) di
chuyển trong giới hạn tiết diện ngang S của vật thể.
Nếu đƣa vào biến số mới :
Thì: (1.12)
từ đó: (1.13)
Cũng nhƣ trong trƣờng hợp ba chiều, nếu đặt điểm quan sát tại gốc tọa độ
tức cho x1=y1=0 thì:
5
(1.14)
Có thể viết lại (1.14) trong hệ tọa độ cực.
Từ hình (1.2) ta có:
Lúc đó (1.14) có dạng:
(1.15)
Trên cơ sở các bài toán tổng quát trên chúng ta xét các bài toán cụ thể:
Hình 1.2: Xác định thế và đạo hàm của vật thể hai chiều 1.3 . Bài toán thuận cho những vật thể có dạng hình học
1.3.1. Hình cầu hoặc điểm vật chất
Trong thực tế, thƣờng gặp các vật thể địa chất tƣơng đối có dạng đẳng thƣớc,
kích thƣớc ngang của chúng theo các hƣớng cùng một bậc. Khi tính toán tác dụng
trọng lực của các vật thể này, ngƣời ta thƣờng xem chúng có dạng hình cầu hoặc là
điểm vật chất. Các vật thể địa chất này thƣờng rất khác nhau: các vật quặng dạng ổ,
dạng bƣớu, các vòm mối, các lỗ hổng cáctơ…
Khảo sát vật thể hình cầu tâm C nằm trong mặt phẳng xoz với các tọa độ xc=
x, yc=y, zc= h (hình 1.3). Khối lƣợng của toàn bộ hình cầu là M, nằm tại tâm hình
cầu. Vì thế ta không phải tính các tính phân khối trên.
6
Hình 1.3: Xác định thế và các đạo hàm của vật thể hình cầu Tại gốc tọa độ:
(1.16) (1.17)
(1.18) (1.19)
(1.20) (1.21)
Để thuận tiện cho việc tính toán ta biến đổi các công thức này khác đi. Đặt
tâm hình cầu dƣới gốc tọa độ (0,0,h), chỉ cần thay đổi dấu của x mà thôi, tức là:
(1.22) (1.23)
(1.24) (1.25)
(1.26) (1.27)
7
Hình 1.4: Trƣờng trọng lực của hình cầu 1.3.2. Thanh vật chất nằm ngang, hình trụ tròn nằm ngang
Các vật thể địa chất dạng này là các cấu tạo dài (các nếp uốn), dạng thấu kính,
các mạch quặng, các vỉa quặng… Hình trụ tròn nằm ngang, nằm dọc theo tâm của
hình trụ.
Nếu một đơn vị độ dài của thanh có khối lƣợng là m thì tƣơng ứng với hình
trụ ta có: =R2 với là khối lƣợng một đơn vị dài.
Trong trƣờng hợp thành phần vật chất nằm ngang ta có thể tính đƣợc giá trị
Vz trực tiếp từ công thức (1.14) mà không cần lấy tích phân, tức là:
(1.28)
Từ đó tìm đƣợc:
(1.29)
(1.30)
(1.31)
Để thuận tiện ta viết lại công thức khi đặt gốc tọa độ trên trục của hình trụ,
còn x là các tọa độ của điểm quan sát. Muốn vậy ta chỉ cần thay đổi x bởi –x trong
các công thức trên là đƣợc:
8
(1.32)
(1.33)
(1.34)
Hình 1.5: Trƣờng trọng lực của hình trụ tròn nằm ngang 1.3.3. Nửa mặt phẳng vật chất nằm ngang
Các dạng vật thể này có biên độ bé, các vùng vót nhọn, các lớp nằm ngang
có độ dày bé…
Giả sử rằng mặt phẳng vật chất nằm ngang có mật độ nằm tại độ sâu h so
với mặt đất, có đƣờng biên song song với trục y, tọa độ ngang của đƣờng biên là
đƣờng x (hình 1.7).
Sử dụng công thức tổng quát (1.14) cho trƣờng hợp này =dz, z=h, ta lấy
tích phân theo x từ x đến +, kết quả thu đƣợc:
(1.35)
9
Hình 1.6: Xác định thế và các đạo hàm của nửa
mặt phẳng vật chất nằm ngang
Từ đó, ngƣời ta tính đƣợc hàm bậc hai của thế trọng lực là:
(1.36)
Tƣơng tự:
(1.37)
Để thuận tiện cho việc tính toán sau này, ta đặt gốc tọa độ tại điểm chiếu của
cạnh bên trên trục x, còn lấy là tọa độ của điểm quan sát. Trong trƣờng hợp này ta
chỉ cần thay đổi dấu của x trong các công thức (1.34), (1.35), (1.36) là đƣợc.
(1.38)
(1.39)
(1.40)
10
Hình 1.7: Trƣờng trọng lực của nửa mặt phẳng vật chất nằm ngang
1.3.4. Hình hộp vuông góc
Nhiều vật thể địa chất gần đúng có thể đƣợc biểu diễn dƣới dạng những khối
bị giới hạn bởi những mặt phẳng, các cấu tạo địa lũy, địa hào, các khối quặng riêng
biệt, những vật thể có thể đƣợc xem là các dạng hình hộp vuông góc. Tính toán tác
dụng trọng lực do hình hộp vuông góc gây ra đƣợc dùng để nghiên cứu các vật thể
khác thƣờng gặp trong thực tế nhƣ bậc thẳng đứng, lớp thẳng đứng. Các công thức
trọng lực của hình hộp vuông góc còn đƣợc sử dụng để tính toán hiệu ứng trọng lực
do các vật thể ba chiều có hình dạng bất kỳ gây ra.
Giả sử có hình hộp vuông góc bị giới hạn bởi các mặt:
Đặt gốc tọa độ tại điểm tính toán. Xuất phát từ các công thức (1.7), (1.8) ta
đƣợc:
(1.41)
(1.42)
(1.43)
(1.44)
Trong đó:
11
Với hình hộp kéo dài ra vô cùng theo trục y, thì từ các công thức, cho y chạy
từ - +, ta thu đƣợc:
(1.45)
(1.46)
(1.47)
1.3.5. Lăng trụ thẳng đứng
Lăng trụ thẳng đứng là vật bị giới hạn bởi hai mặt phẳng thẳng đứng song
song với nhau và một mặt phẳng nằm ngang. Trong trƣờng hợp này ta xem z2=,
z1=h.
Với điều kiện này chúng ta thấy rằng g = , hơn nữa cho độ dày lớp bằng 2d
thì:
Chuyển gốc tọa độ về tạo hình chiếu của trọng điểm trên của lớp trên mặt đất
(thay đổi dấu của x). Từ công thức (1.41), chúng ta thu đƣợc các kết quả sau đây:
(1.48)
(1.49)
1.3.6. Bậc thẳng đứng
Trong lý thuyết phân tích các dị thƣờng trọng lực, ngƣời ta hiểu bậc thẳng
đứng là vật thể hai chiều bị giới hạn bởi hai mặt phẳng song song vô hạn và một
mặt thẳng đứng. Tiết diện ngang của của vật thể này có dạng dải vuông góc vô
cùng, có các cạnh song song với trục x và trục z.
12
Từ (hình 1.8a và 1.8b) ta thấy rằng bậc thẳng đứng tƣơng tự nhƣ mặt phẳng
vật chất nằm ngang nhƣng nó tổng quát và phức tạp hơn.
Hình 1.8 a, 1.8 b: Bậc thẳng đứng
Theo định nghĩa trên, ngƣời ta sử dụng mật độ dƣ của các bậc thẳng đứng là
khác không còn toàn bộ không gian là bằng không. Nhƣng nếu phần không gian
dƣới bậc có mật độ dƣ khác không thì dị thƣờng trọng lực gây ra vẫn không thay
đổi.
Trong thực tế các vật địa chất dạng này là các cấu tạo tiếp xúc với vòm muối
hoặc là các khối xâm nhập với các vùng đất đá vây quanh. Có thể nói rằng bậc
thẳng đứng là một trong những vấn đề cơ bản của lý thuyết và thực tế khi phân tích
dị thƣờng trọng lực.
Gọi h là tọa độ ngang của đƣờng biên của bậc, h1 và h2 là độ sâu đến các mép
giới hạn trên và dƣới của bậc. Để tính hiệu ứng trọng lực trong trƣờng hợp này,
ngƣời ta lấy tích phân các công thức tổng quát (1.6), (1.7), (1.8), (1.9) theo các biến
x, y, z, sau đó thay cận tích phân y=, x đến và z từ h1 đến h2. Cụ thể là:
(1.50)
(1.51)
(1.52)
(1.53)
13
Trong đó:
Trƣớc hết cho y1=-, y2=+ ta có:
(1.54)
(1.55)
(1.56)
Các công thức trên biểu thị tác dụng trọng lực của hình hộp chạy dài vô cùng
theo trục y với tiết diện ngang là hình chữ nhật có các cạnh bằng (x1, x2) và (h1, h2).
Để tính đƣợc hiệu ứng trọng lực của bậc thẳng đứng, trong các công thức trên thay
x1=x, x2=. Kết quả chúng ta thu đƣợc:
(1.57)
(1.58)
(1.59)
Dịch chuyển gốc tọa độ về điểm nằm trên biên của bậc thẳng đứng (điểm x,
0) và thay cận lấy tích phân, ta có:
(1.60)
(1.61)
(1.62)
14
Hình 1.9: Trƣờng trọng lực trên bậc thẳng đứng 1.3.7. Bậc nghiêng
Trong thực tế trƣờng hợp bậc nghiêng thƣờng gặp nhiều hơn bậc thẳng
đứng.
Trong trƣờng hợp này có thể dùng các công thức khác nhau để tính các đại
lƣợng Vz, Vxz, Vzz. Dƣới đây chúng ta sẽ xét phƣơng pháp đƣợc xem là tốt hơn cả.
Trên hình vẽ chúng ta vẽ mặt cắt của bậc nghiêng. Lấy điểm quan sát làm gốc tọa
độ. Tọa độ giao điểm của tuyến với đƣờng biên của bậc nghiêng kéo dài là x. Tại độ
sâu bất kỳ là Z( h1 độ biên của lớp nghiêng đó là (x-zcotg, z). Hình 1.10: Bậc nghiêng Dùng các công thức của nửa mặt phẳng vật chất nằm ngang (1.35), (1.36) và (1.37) chúng ta có thể tính đƣợc hiệu ứng trọng lực của bậc nghiêng. (1.63) (1.64) (1.65) Sau khi lấy tích phân và đƣa gốc tọa độ từ điểm quan sát đến điểm 0 bằng cách thay dấu biến số x, ta thu đƣợc: (1.66) (1.67) (1.68) Các công thức (1.66), (1.67) và (1.68) khá phức tạp nên rất khó đề ra đƣợc phƣơng pháp giải tích để tích phân. Trong bài toán này, để giải bài toán thuận, ngƣời ta biến đổi các công thức trên. Từ hình vẽ ta thấy các thông số , liên hệ với tọa độ vuông góc bằng biểu thức: Nếu gọi 1, 1, 2, 2 là các thông số tƣơng ứng với các góc trên và dƣới của bậc ta có thể viết lại công thức (1.66 - 1.68) ở dạng: (1.69) (1.70) (1.71) Đối với các bậc nghiêng, ngƣời ta không để ra phƣơng pháp giải tích để xác định các thông số của chúng mà thông thƣờng so sánh dạng đƣờng cong của chứng từ bậc nghiêng với các góc khác nhau. Chƣơng 2 PHƢƠNG PHÁP XÁC ĐỊNH DỊ THƢỜNG TRỌNG LỰC CỦA CÁC RANH GIỚI 2D VÀ 3D 2.1. Phƣơng pháp xác định dị thƣờng trọng lực của ranh giới trầm tích 2D trong miền không gian 2.1.1. Xác định dị thƣờng trọng lực của ranh giới trầm tích trên cơ sở phân chia nó thành các đa giác có tiết diện ngang bất kỳ Theo Talwari và Ewing [18] dị thƣờng trọng lực của vật thể có tiết diện ngang bất kỳ và mật độ dƣ thay đổi theo chiều sâu có thể đƣợc xác định bằng cách chia vật thể thành những lớp nằm ngang rồi lấy tổng dị thƣờng trọng lực của chúng. Hình 2.1: Vật thể hai chiều có tiết diện
ngang bất kỳ
và mật độ dƣ (z) là hàm của độ sâu z. Dị thƣờng trọng lực dg của lớp nằm ngang có chiều dày dz, nằm ở độ sâu z và đƣợc giới hạn bởi chu vi của vật thể đƣợc xác định bởi: dg = 2f(z )(2 - 1) dz (2.1) ở đây 2 - 1 là góc nhìn từ điểm cần tính dị thƣờng trọng lực P(0,0) tới lớp nằm ngang (Hình 2.1), (z) là mật độ dƣ của vật thể đƣợc xem nhƣ là một hàm của chiều sâu z, f là hằng số hấp dẫn. Dị thƣờng trọng lực g của toàn bộ vật thể tại điểm P(0,0) đƣợc xác định bằng cách lấy tích phân vế phải của đẳng thức (2.1) trong phạm vi từ zđỉnh đến zđáy,đó tƣơng ứng là chiều sâu tới đỉnh và đáy vật thể. g = -2f [ (z) 1dz + (z) 2 dz] = -2f (z) dz (2.2) Đẳng thức (2.2) chỉ ra rằng việc tính toán dị thƣờng trọng lực của vật thể 2 chiều có mật độ dƣ thay đổi theo độ sâu đƣợc thực hiện bằng cách lấy tích phân đƣờng dọc theo chu vi của vật thể theo chiều kim đồng hồ. Trên cơ sở công thức này, Murthy I..V.R và Bhaskara Rao D [14] đã đƣa ra các thuật toán xác định dị thƣờng trọng lực của vật thể 2 chiều có mật độ dƣ thay đổi ở dạng hàm mũ nhƣ dƣới đây. 2.1.2. Trƣờng hợp mật độ dƣ thay đổi tuyến tính theo độ sâu Trong trƣờng hợp này dị thƣờng trọng lực đƣợc tính bằng cách xấp xỉ tiết diện ngang của vật thể bởi một đa giác N cạnh ABCDEF (Hình 2.2). Tọa độ (xk, yk) của các đỉnh A, B, C... đƣợc tính trong hệ tọa độ mà gốc đặt tại điểm cần tính dị thƣờng trọng lực P(0,0). Giả sử rằng mật độ dƣ thay đổi theo quy luật: (z) = (0) + az (2.3) trong giới hạn của vật thể. Ở đây (0) là giá trị của mật độ dƣ tại mặt quan sát còn a biểu thị tốc độ biến đổi theo chiều sâu của mật độ dƣ. Trong trƣờng hợp này để tính dị thƣờng trọng lực của vật thể ta tiến hành tính tích phân trong vế phải của đẳng thức (2.1) dọc theo mỗi cạnh của vật thể (ví dụ cạnh BC) rồi sau đó lấy tổng các giá trị này. Ta có: (2.4) Thay (z) từ (2.3) vào (2.4) rồi thực hiện việc lấy tích phân ta đƣợc: kk)} (2.5) k+1k+1-z2 dgBC = 2f { (0).A[sinik ln ( rk+1 / rk) - cos ik(k+1-k)]+ (a/2)(zk+1-zk) sinik
- (a/2) A2[cos 2ikln(rk+1/rk) - sin 2ik(k+1- k)]-[ (0) (zk+1k+1-zkk)
- (a/2) (z2 Hình 2.3: Việc phân chia mỗi
cạnh đa giác
giác thành các đoạn. Hình 2.2: Xấp xỉ vật thể có tiết diện
ngang
bất kỳ bằng đa giác N cạnh. Trong đó: sin ik = (zk+1) / k = x2 k + z2 k cos ik = (xk+1-xk) /
A = xksin ik - zk cos ik; r2 k = /2 - arctg(xk/zk) với z 0 = /2(1-xk/xk với z =0 Đối với rk+1 ta có các công thức tƣơng tự. 2.1.3. Trƣờng hợp mật độ dƣ thay đổi theo quy luật hàm mũ theo chiều sâu Việc xác định hiệu ứng trọng lực do vật thể hai chiều có hình dạng bất kỳ và mật độ dƣ thay đổi theo quy luật hàm số mũ theo chiều sâu: (z ) = (o) e- z (2.6) đƣợc thực hiện nhƣ sau: . Xấp xỉ vật thể bằng một đa giác N cạnh. . Chia mỗi cạnh của đa giác thành s đoạn nhỏ và giả sử rằng trên mỗi đoạn đó mật độ dƣ thay đổi một cách tuyến tính. Nếu (xk, yk) và (xk+1, yk+1) tƣơng ứng là tọa độ hai đỉnh của một cạnh nào đó của đa giác N cạnh (ví dụ cạnh BC) (Hình 2.3) thì tọa độ (x'j, z'j) của các đoạn đƣợc chia ra trên cạnh đó là: x'j= x k + [(xk+1 - xk) / s ](j -1) với j = 1,2...s+1 (2.7) z'j= zk + [( zk+1 - zk) / s ]( j - 1) với j = 1,2...s+1 (2.8) Giả sử rằng trên mỗi đoạn chia mật độ dƣ đƣợc xác định bởi (z) = (0) +a z. Vì rằng (z'j+1) và (z'j) đƣợc tính một cách dễ dàng từ (2.6) nên các giá trị (0) và a trên mỗi đoạn đó đƣợc xác định nhƣ sau: a = [(z'j+1) - (z'j) ] / (z'j +1 - z'j ) (2.9) (0) = (z'j) - a z'j (2.10) Thay các giá trị (0), a, ( (x'j, z' j ) và ( x'j+1, z'j+1 ) vào công thức (2.5) ta tính đƣợc hiệu ứng trọng lực của từng đoạn nhỏ đƣợc chia ra trên cạnh BC của đa giác. Dị thƣờng trọng lực dg của cả cạnh BC của đa giác sẽ đƣợc tính bằng cách lấy tổng s lần tính giá trị trọng lực của các đoạn này. Quá trình đƣợc tiến hành tƣơng tự cho các cạnh khác của đa giác. Kết quả dị thƣờng trọng lực do toàn bộ vật thể gây ra là: (2.11) Thật rõ ràng trong trƣờng hợp vật thể có mật độ dƣ thay đổi theo chiều sâu theo quy luật hàm mũ thì thời gian tính dị thƣờng trọng lực sẽ s lần lớn hơn thời gian tính trong trƣờng hợp mật độ dƣ thay đổi tuyến tính theo chiều sâu. 2.1.4. Mật độ dƣ thay đổi theo dạng hàm hypepol Theo Litinsky, sự thay đổi mật độ dƣ của vật thể theo chiều sâu có thể đƣợc xấp xỉ bằng một hàm hypepol nhƣ sau: (2.12) Trong đó (z) là mật độ dƣ ở độ sâu z, o là mật độ dƣ đƣợc ngoại suy đối với mặt đất còn là tốc độ biểu diễn sự thay đổi của mật độ dƣ theo độ dài. Dị thƣờng trọng lực g(0) ở điểm P(0) đƣợc xác định nhƣ sau: (2.12a) Trong đó G: là hằng số hấp dẫn Thay thế (z) trong phƣơng trình (2.12a) bằng các giá trị của phƣơng trình (2.12) có thể viết đƣợc nhƣ sau: (2.12b) áp dụng công thức Stoke, phƣơng trình (2.12b) có thể đƣợc biến đổi: (2.12c) Thay cho việc tính tích phân này, ta tính gần đúng nó bằng tổng sau: (2.13) Với: Trong đó: và: , và i đƣợc giải thích trong hình (2.3). Các ký hiệu xk, xk+1, zk, zk+1, Hai số hạng đầu tiên trong phƣơng trình (2.13) có thể đƣợc giản ƣớc đi bởi khi xét chung cho tất cả các cạnh của vùng thì: chúng sẽ biến mất. Vì vậy, phƣơng trình (2.13) đƣợc viết lại dƣới dạng đơn giản là: (2.14) 2.1.5. Mật độ dƣ thay đổi theo dạng hàm parabolic Sự thay đổi mật độ dƣ theo độ sâu cũng có thể đƣợc biểu diễn bởi hàm có dạng parabolic sau đây: (2.15) trong đó là hằng số với đơn vị là độ dài nghịch đảo. Tƣơng tự nhƣ trƣờng hợp trƣớc, ta có: (2.16) với: (2.17) Trong đó: và: B = - 2cos cosi -2 o Phƣơng trình này có thể đƣợc viết trong dạng đơn giản bởi giới hạn của số hạng đầu tiên bởi vì tổng kết của chúng cho tất cả các mặt của vùng sẽ biến mất. Vì vậy, chúng ta có: (2.18) 2.1.6. Xác định dị thƣờng trọng lực của ranh giới trầm tích trên cơ sở phân chia nó thành các lăng trụ thẳng đứng đặt cạnh nhau Thông thƣờng trong trƣờng hợp bài toán hai chiều việc xác định dị thƣờng trọng lực do một đối tƣợng địa chất có tiết diện ngang bất kỳ gây ra đƣợc thực hiện bằng cách xấp xỉ tiết diện ngang của nó bằng một đa giác N cạnh. Nhƣ vậy thực chất của việc giải bài toán ngƣợc là xác định vị trí các đỉnh của đa giác sao cho sự sai lệnh giữa dị thƣờng quan sát và tính toán là nhỏ nhất. Với các phƣơng pháp này quá trình tính toán đòi hỏi đƣa vào các tọa độ đỉnh tiên nghiệm của các đa giác và chúng phải đủ gần với các tọa độ thật thì phƣơng pháp mới có độ hội tụ tốt. Tuy nhiên trên thực tế việc đảm bảo điều kiện này là một việc hết sức khó khăn. Trong mục này, ta áp dụng một phƣơng pháp khác để tính dị thƣờng trọng lực, đó là phƣơng pháp do Murthy I.V.R đƣa ra. Ở đây thay thế cho việc xấp xỉ bằng một đa giác, đối tƣợng gây dị thƣờng trọng lực đƣợc xấp xỉ bằng một chuỗi các lăng trụ thẳng đứng đặt cạnh nhau. Trên cơ sở phƣơng pháp này ta tiến hành việc giải bài toán ngƣợc theo phƣơng pháp lựa chọn nhằm xác định độ sâu tới ranh giới phân chia mật độ ở từng điểm quan sát trên tuyến. Ngoài độ sâu trung bình Z và mật độ dƣ quá trình tính toán không đòi hỏi đƣa vào bất kỳ một giá trị tiên nghiệm nào khác về mặt ranh giới phân chia. Hình 2.4: Ranh giới phân chia mật độ và việc xấp xỉ nó bằng
các lăng trụ
trụ thẳng đứng đặt cạnh nhau Xét một mặt phân chia mật độ có dạng nhƣ (Hình 2.4), ở đây Z là độ sâu trung bình của ranh giới mà trên đó có sự thay đổi cấu trúc địa chất gây nên dị thƣờng trọng lực, xk là điểm quan sát thứ k trên tuyến. Giả sử rằng tuyến quan sát đủ dài sao cho có thể phủ hết cả phần thay đổi cấu trúc địa chất. Để xác định hiệu ứng trọng lực do ranh giới phân chia mật độ này gây ra ta xấp xỉ nó bằng một chuỗi các lăng trụ có hình chiếu tâm lên tuyến quan sát trùng với vị trí các điểm quan sát, có độ rộng bằng khoảng cách dx giữa các điểm quan sát (Hình 2.4). Bằng cách nhƣ vậy việc xác định hiệu ứng trọng lực của ranh giới phân chia mật độ thực chất đƣợc thay thế bằng việc tính tổng hiệu ứng trọng lực của các lăng trụ có chiều cao khác nhau. Tại điểm quan sát P(xk) trên tuyến ta có: (2.19) Trong đó Axk + B là phần trƣờng khu vực giả sử có dạng tuyến tính còn quan sát thứ i. Theo Rao B.S.R và Murthy I.V.R [14] ta có: Fk (z) = 2f{z[arctg((xk+dx/2)/z)-arctg((xk-dx/2)/z)]+
+0.5[(xk+dx/2)ln((xk+dx/2)2+z2)-(xk-dx/2)ln((xk-dx/2)2+z2)]} (2.20) ở đây dx là khoảng cách giữa các điểm quan sát trên tuyến, ZT(i) là độ sâu tới mặt ranh giới phân chia mật độ nằm dƣới điểm quan sát thứ i, f là hằng số hấp dẫn, là mật độ dƣ. 2.2. Phƣơng pháp xác định dị thƣờng trọng lực của ranh giới trầm tích 3D trong miền không gian 2.2.1. Cơ sở lý thuyết Theo Bhaskara, D.,1986 sự thay đổi của mật độ dƣ theo độ sâu có thể xấp
xỉ bằng một hàm bậc hai : (z) = ao+ a1z + a2z2 (2.21) trong đó z là độ sâu ;a1,a2, là các hệ số suy giảm;ao là mật độ dƣ của lớp trầm tích bề mặt . Để xác định dị thƣờng trọng lực, bể trầm tích đƣợc chia thành các lăng trụ đặt sát nhau ở phía dƣới mỗi điểm quan sát, có độ rộng bằng khoảng cách giữa các điểm quan sát và có chiều cao bằng độ sâu của bể trầm tích tại điểm đó. Dị thƣờng trọng lực gây ra bởi bể trầm tích đƣợc tính bằng cách lấy tổng dị thƣờng trọng lực của tất cả các lăng trụ này. Hình 2.5 – Mô hình lăng trụ 3 chiều Dị thƣờng trọng lực của mỗi lăng trụ đƣợc xác định theo công thức sau (Prakash và Ramesh Babu N. ,1990 ): (2.22) ở đây f là hằng số hấp dẫn; Z1,Z2 tƣơng ứng là độ sâu đến đỉnh và đáy của lăng trụ. T và W tƣơng ứng là nửa bề rộng của lăng trụ theo các trục x và y . KÕt qu¶ cña viÖc tính tích phân này sau khi thay (z) từ công thức (2.21) là : +fa1 - , (2.23) +fa2 trong đó: X1 = x+T, X2 = x-T, Y1 = y+W, Y2 = y-W và R = Ta nhận thấy rằng công thức này rất phức tạp. Khi tính dị thƣờng trọng lực của bể trầm tích nêú dị thƣờng của tất cả các lăng trụ đều đƣợc tính theo công thức này thì sẽ tốn rất nhiều thời gian tính trên máy. Để khắc phục khó khăn này, V.D.Braskara Rao đề nghị sử dụng công thức gần đúng sau đây để tính dị thƣờng của các lăng trụ nằm xa điểm quan sát. Công thức này đạt đƣợc khi xấp xỉ hiệu ứng của lăng trụ bằng một thanh vật chất thẳng đứng đặt tại tâm của lăng trụ. +fa2xy dg(x,y)=fa0xy (2.24) Hay dg(x,y)=fa0xy+fa1xy ao là mật độ dƣ của lớp trầm tích bề mặt, a1 là hệ số suy giảm, ở đây R = còn x,y tƣơng ứng là khoảng cách giữa các điểm quan sát theo các trục x và y. Cuối cùng, dị thƣờng trọng lực của bể trầm tích đƣợc xác định theo công thức sau: (2.25) Trong đó: M là số lăng trụ đƣợc chia theo trục x. N lµ sè l¨ng trô ® îc chia theo trôc y Theo Talwani and Ewing (1960) đã chia môi trƣờng thành rất nhiều lớp mỏng và mỗi lớp này lại đƣợc xấp xỉ là một đa giác m cạnh. Do vậy ảnh hƣởng của mỗi lớp mỏng đƣợc xác định theo công thức :. ; (2.26) Trong đó V là dị thƣờng do mỗi lớp mỏng có bề dày gây ra. V đƣợc biều diễn bởi 1 tích phân mặt. Talwani và Ewing biểu diễn V nhƣ sau : (2.27) Hình 2.6. Mô hình khối đa diện Trong đó là mật độ khối của lớp mỏng, z là độ sâu của lớp, là tọa độ điểm thứ i của lớp mỏng. nếu , nếu nếu nếu (2.28) 2.2.2. Xây dựng chƣơng trình giải bài toán ngƣợc 3D và tính toán thử nghiệm trên mô hình. ViÖc gi¶i bµi to¸n ng îc ® îc thùc hiÖn theo ph ¬ng ph¸p lùa chän. §¸nh gi¸ tiªn nghiÖm vÒ ®é s©u cña bÓ trÇm tÝch ® îc x¸c ®Þnh theo ph ¬ng ph¸p ®¸nh gi¸ ®é s©u trùc tiÕp cña Bott. NÕu ký hiÖu (i,j) lµ täa ®é t©m cña l¨ng trô thø i,j th× Z (i , j) = (2.29) Sau khi t×m ® îc Z(i,j) trªn nghiÖm, theo c¸c c«ng thøc (2.23), (2.24) vµ (2.25) ta x¸c ®Þnh ® îc dÞ th êng träng lùc g (i ,j) t¹i tÊt c¶ c¸c ®iÓm (i ,j) trªn mÆt ph¼ng quan s¸t. Chªnh lÖch gi÷a dÞ th êng quan s¸t vµ dÞ th êng tÝnh to¸n ® îc sö dông ®Ó ®iÒu chØnh ®é s©u cña c¸c l¨ng trô, tøc lµ ®iÒu chØnh ®é s©u cña bÓ trÇm tÝch sau mçi lÇn lùa chän . (2.30) víi: (z) = a0 + a1z + a2z2 Qu¸ tr×nh lùa chän dõng l¹i khi sai sè b×nh ph ¬ng trung b×nh gi÷a dÞ th êng quan s¸t vµ tÝnh to¸n nhá h¬n sai sè cho phÐp hoÆc sè lÇn lùa chän v ît qu¸ sè lÇn lùa chän cho tr íc. §Ó lµm gi¶m thêi gian tÝnh trªn m¸y, qu¸ tr×nh lùa chän ® îc chia lµm nhiÒu b íc. Trong mçi b íc vïng c¸c l¨ng trô ® îc sö dông c«ng thøc chÝnh x¸c (2.23) ®Ó tÝnh dÞ th êng träng lùc t¨ng tØ lÖ víi sè thø tù cña b íc, c¸c l¨ng trô cßn l¹i (xa ®iÓm quan s¸t ) ® îc tÝnh theo c«ng thøc gÇn ®óng (2.24). ChØ riªng ë lÇn lùa chän cuèi cïng dÞ th êng cña toµn bé c¸c l¨ng trô míi ® îc tÝnh theo c«ng thøc chÝnh x¸c (2.23). Trªn c¬ së c¸c thuËt to¸n ®· tr×nh bµy, ta tiÕn hµnh viÖc tÝnh to¸n thö nghiÖm trªn mét m« h×nh bÓ trÇm tÝch ba chiÒu (h×nh 2.7) cã sù thay ®æi mËt ®é d theo chiÒu s©u. Sö dông ph ¬ng ph¸p b×nh ph ¬ng tèi thiÓu ta x¸c ®Þnh ® îc c¸c hÖ sè cña hµm suy gi¶m mËt ®é d bËc hai lµ:
(z) = - 0,4937 - 0,0749 z + 0,0042z 2 ; Rms = 0,0310 ở đây, khoảng cách giữa các điểm quan sát dx và dy đƣợc lấy là 4 km theo cả hai chiều x và y. Số lần lƣạ chọn đƣợc lấy là 10. Quá trình lựa chọn đƣợc chia làm bốn bƣớc nhƣ sau: 1. Hai lần lựa chọn đầu tiên, bán kính vùng lăng trụ đƣợc tính dị thƣờng theo công thức chính xác là 1,5 dx. 2. Năm lần lựa chọn tiếp theo, bán kính vùng các lăng trụ đƣợc tính dị thƣờng theo công thức chính xác là 3,5 dx. 3. Hai lần lựa chọn tiếp theo, bán kính vùng các lăng trụ đƣợc tính dị thƣờng theo công thức chính xác là 5,5 dx 4. Lần lựa chọn cuối cùng, tất cả các lăng trụ đƣợc tính dị thƣờng theo công thức chính xác. Kết quả tính toán thử nghiệm trên mô hình bao gồm: dị thƣờng trọng lực của bể trầm tích, độ sâu đáy bể trầm tích ở lần lặp cuối , dị thƣờng tính toán ở lần lặp cuối, sai số về độ sâu (so với mô hình) ở lần lặp cuối, tốc độ hội tụ của quá trình lựa chọn tƣơng ứng đƣợc đƣa ra trên hình 2.7 a) b) a) b) c) d) e) f) Hình 2.7 - Kết quả xác định độ sâu bể trầm tích trong miền không gian a - Mô hình bể trầm tích (km); b - Dị thường quan sát tính từ mô hình (mgal) c - Bể trầm tích ở lần lặp cuối (km); d - Dị thường tính toán ở lần lặp cuối(mgal) e - Độ lệch giữa mô hình và kết quả tính (km); f -Tốc độ hội tụ của phương pháp 2.3. Các phƣơng pháp xác định độ sâu của bể trầm tích 3D trong miền tần số Trong các phƣơng pháp xác định độ sâu bể trầm tích trong miền không gian vấn đề thay đổi mật độ dƣ của bể ở dạng hàm số mũ theo độ sâu mới chỉ đƣợc các tác giả giải quyết một cách gần đúng bằng cách qui nó về dạng thay đổi tuyến tính trên những khoảng chia nhỏ của độ sâu (Murthy) hoặc xấp xỉ nó bằng một hàm bậc hai (Rao và Ramesh Babu). Hơn nữa, đối với bài toán ba chiều thực tế tính toán trên mô hình cho thấy rằng mặc dù đã sử dụng một số thủ thuật trong quá trình giải bài toán ngƣợc, thời gian tính trên máy vẫn còn rất lớn. Những khó khăn này đƣợc Chai và Hinze khắc phục bằng cách tính toán trong miền tần số. Dị thƣờng trọng lực trong miền không gian đƣợc xác định nhờ phép biến đổi Fourier ngƣợc. Độ chính xác của phƣơng pháp này phụ thuộc vào độ chính xác của việc biến đổi giữa miền không gian và miền tần số Vì vậy để có thể áp dụng đƣợc phƣơng pháp này vào việc xử lý các tài liệu thực tế thì vấn đề đầu tiên cần giải quyết là phải tăng đƣợc độ chính xác này. Về mặt lý thuyết ta có thể khắc phục điều đó bằng cách tăng mật độ các điểm quan sát và chiều dài của tuyến. Trên thực tế cả hai cách này không phải lúc nào cũng thực hiện đƣợc. Tuy nhiên độ chính xác này cũng có thể đƣợc tăng lên đáng kể nhờ việc áp dụng phƣơng pháp "trƣợt mẫu" dựa trên cơ sở phƣơng trình độ lệch của biến đổi Fourier rời rạc. 2.3.1. Nâng cao độ chính xác của việc tính dị thƣờng trọng lực trong miền tần số bằng phƣơng pháp "trƣợt mẫu" (Shift-sampling) Giả sử X(x) là hàm xác định trong miền không gian thì phổ Fourier (f) của nó đƣợc xác định bởi: (2.31) Trên cơ sở các định lý "trƣợt" trong cả miền không gian và miền tần số và lý thuyết "lấy mẫu" (sampling theorems) cho các hàm tuần hoàn, Chai và Hinze [11] đã đƣa ra phƣơng trình độ lệch của phép biến đổi Fourier ngƣợc rời rạc nhƣ sau: (2.32) Trong đó: A là độ lệch gây nên bởi sự rời rạc hóa B là độ lệch gây nên sự giới hạn của tuyến quan sát FT là ký hiệu của phép biến đổi Fourier DFT là ký hiệu của phép biến đổi Fourier rời rạc IDFT là ký hiệu của phép biến đổi Fourier ngƣợc rời rạc ,d tƣơng ứng là khoảng cách lấy mẫu trong miền không gian và miền tần số. và tƣơng ứng là bƣớc trƣợt trong miền không gian và miền tần số. Phƣơng trình này đúng cho mọi hàm thỏa mãn điều kiện Dirichlet. Với phƣơng trình này ảnh hƣởng do sự rời rạc hóa và sự giới hạn của tuyến quan sát trong phép biến đổi Fourier ngƣợc đƣợc thể hiện một cách rõ ràng dƣới dạng một công thức toán học. Việc áp dụng phƣơng trình (2.32) đối với các dị thƣờng của trƣờng thế đƣợc căn cứ vào các tính chất sau đây của chúng: a - Phổ của các dị thƣờng của trƣờng thế giảm rất nhanh theo sự tăng lên của tần số. b - Các dị thƣờng của trƣờng thế giảm một cách đơn điệu khi khoảng cách tới nguồn gây dị thƣờng tiến đến vô cùng. c - Giá trị tuyệt đối của các đạo hàm dị thƣờng của trƣờng thế cũng giảm một cách đơn điệu khi khoảng cách tới nguồn gây dị thƣờng tiến đến vô cùng. Tính chất (a) cho thấy rằng ảnh hƣởng do sự giới hạn của tuyến quan sát có thể bỏ qua khi chiều dài tuyến đủ lớn. Khi cho = 0 phƣơng trình (2.32) trở thành: (2.33) Trong đó ak(n) =X[(n+kN)]+X[(n-kN)] còn E(n,) chính là độ lệch giữa các giá trị tính toán đƣợc nhờ phép biến đổi Fourier ngƣợc và các giá trị đúng. Vấn đề đặt ra là cần phải tìm bƣớc "trƣợt" sao cho độ lệch bình phƣơng trung bình của chúng có giá trị nhỏ nhất. Việc khảo sát hàm E(n,) căn cứ vào các tính chất (b) và (c) cho thấy sai số bình phƣơng trung bình: (2.34) có hai cực tiểu xung quanh các giá trị = 0.25 và = 0.75. Qua việc tính toán trên nhiều dị thƣờng trọng lực gây nên bởi các đối tƣợng gây dị thƣờng khác nhau Chai và Hinze đã đƣa ra bƣớc trƣợt tối ƣu là = 0.22 đối với bài toán hai chiều và = 0.26 đối với bài toán ba chiều. 2.3.2. Xác định dị thƣờng trọng lực của ranh giới 3D trọng lực trong miền tần số Xét một yếu tố khối dm đƣợc đặt tại điểm có tọa độ (,,z). Hiệu ứng trọng lực dg do nó gây ra tại điểm P có tọa độ (x,y,0) là: (2.35) còn phổ hiệu ứng trọng lực của dm là: d (u,v) = 2f dm exp (-BC) (2.36) Trong đó: f là hằng số hấp dẫn; B = 2 (iu,iv,s); C = [ ] với u, v tƣơng ứng là tần số theo các trục x và y trong miền không gian. Bằng việc lấy tích phân đẳng thức (2.36) trong phạm vi một hình lăng trụ chữ nhật thẳng đứng (lăng trụ thứ m, n) có đáy nằm ở độ sâu Zmn, đỉnh trùng với mặt quan sát z = 0, tọa độ tâm (m, n), bề rộng tƣơng ứng theo các trục x, y là 2a, 2b với giả thiết mật độ dƣ của nó thay đổi theo quy luật hàm mũ theo chiều sâu với hệ số suy giảm là m,n, Chai và Hinze đã thu đƣợc công thức tính phổ dị thƣờng trọng lực của lăng trụ này là: x (2.37) Việc tính phổ dị thƣờng trọng lực của bể trầm tích đƣợc thực hiện bằng cách lấy tổng phổ dị thƣờng trọng lực của tất cả các lăng trụ này: (2.38) Dị thƣờng trọng lực của bể trầm tích đƣợc xác định nhờ phép biến đổi Fourier ngƣợc theo công thức sau: (2.39) Trên thực tế sự suy giảm mật độ dƣ của bể trầm tích theo chiều sâu có xu hƣớng tiệm cận tới giá trị d không đổi nào đó. Bởi vậy sự suy giảm mật độ dƣ của bể trầm tích theo chiều sâu có thể đƣợc xấp xỉ bằng biểu thức: (2.40) (z) = d + (0) ở đây dị thƣờng trọng lực do phần mật độ dƣ d không đổi của lăng trụ gây ra có thể đƣợc tính bởi công thức sau: (Baneree và Gupta). (2.41) Trong đó: dg mn (x,y) = fd [xln(y+R)+yln(x+R) (2.42) với: R2 = x2+y2+z2 Cuối cùng, dị thƣờng trọng lực của bể trầm tích tại điểm (x,y,0) đƣợc xác định bởi: g(x,y) = gp(x,y) +gd(x,y) (2.43) Xét trƣờng hợp bài toán hai chiều: Đối với trƣờng hợp bài toán hai chiều, việc lấy tích phân đẳng thức (2.26) theo trục y đƣợc thực hiện từ - đến + . Kết quả thu đƣợc: x{1-exp[-2u+m)zm}[exp(-2ium)] (2.44) (2.45) (2.46) Trong đó: (2.47) gm(x)=-fd[xln(x2+z2)+2z ctg] 2.3.3 Xây dựng chƣơng trình giải bài toán ngƣợc và tính toán thử nghiệm trên các mô hình Dựa vào phƣơng pháp tính phổ dị thƣờng trọng lực của bể trầm tích hai và ba chiều với mật độ dƣ thay đổi của Chai và Hinze đã trình bày ở trên, ta tiến hành việc xây dựng chƣơng trình giải bài toán ngƣợc xác định độ sâu của đáy bể trầm tích trên mô hình 3D cụ thể có các thông số sau : - Số điểm quan sát theo trục x: 64 điểm - Số điểm quan sát theo trục y: 64 điểm - Khoảng cách dx,dy giữa các điểm quan sát: 1 km - Mật độ dƣ : (z) = -0,48e-0,15z Quá trình giải bài toán ngƣợc trong miền tần số theo phƣơng pháp lựa chọn đƣợc tiến hành theo các bƣớc sau: 1- Từ các giá trị dị thƣờng quan sát (trong trƣờng hợp này chúng đƣợc tính từ hai mô hình nêu trên) xác định nghiệm ban đầu theo phƣơng pháp xác định trực tiếp độ sâu của Bott : z(m,n) = ( )ln [1+gqs(m,n) / 2f (0)] (2.48) 2 - Tính phổ dị thƣờng trọng lực của bể trầm tích đối với bài toán hai chiều và theo công thức (2.51) đối với bài toán ba chiều với bƣớc "trƣợt " phổ tƣơng ứng là 0,22 và 0,26. 3 - Kết quả của phép biến đổi Fourier ngƣợc đƣơc nhân tƣơng ứng với các thừa số : e 2 i 0,26 (n+n) / MN đối với bài toán ba chiều. Phần thực của kết quả thu đƣợcchính là dị thƣờng gtt của bể trầm tích . 4 - Điều chỉnh độ sâu của bể trầm tích sau mỗi lần lựa chọn dựa vào độ sai lệch giữa dị thƣờng quan sát và dị thƣờng tính toán dz(m,n) = (2.49) 5 - Quay lại thực hiện các bƣớc 2,3,4 khi sai số bình phƣơng trung bình giữa dị thƣờng quan sát và tính toán chƣa nhỏ hơn sai số cho phép hoặc số lần lựa chọn chƣa bằng số lần lựa chọn cho trƣớc . Các kết quả tính toán trên mô hình 3 chiều đƣợc đƣa ra trên hình 2.8 a) b) c) d) c) d) e) f) Hình 2. 8 - Kết quả xác định đố sâu bể trầm tích trong miền tần số a - Mô hình bể trầm tích (km); b - Dị thường quan sát tính từ mô hình (mgal) c - Bể trầm tích ở lần lặp cuối (km); d - Dị thường tính toán ở lần lặp cuối(mgal) e - Độ lệch giữa mô hình và kết quả tính (km); f -Tốc độ hội tụ của phương pháp Chƣơng 3 MÔ HÌNH HÓA VIỆC GIẢI BẢI TOÁN NGƢỢC XÁC ĐỊNH SỰ PHÂN BỐ MẬT ĐỘ MÓNG KẾT TINH 3.1. Các phƣơng pháp giải bài toán ngƣợc Khi xem dị thƣờng trọng lực quan sát nhƣ là một trƣờng tổng bao gồm dị thƣờng gây ra bởi các lớp bên dƣới. Dựa vào các thuật toán xác định hiệu ứng trọng lực của ranh giới trầm tích với mật độ dƣ thay đổi và bài toán ngƣợc xác định ranh giới trầm tích của Murthy và Rao đã trình bày ở trên, chúng tôi đƣa ra thuật toán giải bài toán ngƣợc nhằm xác định sự phân bố mật độ của đá móng trên cơ sở tính "bóc lớp" dị thƣờng trọng lực gây ra bởi các ranh giới, cơ sở lý thuyết của nó đƣợc trình bày dƣới đây: 3.1.1. Xác định sự phân bố mật độ đá móng theo phƣơng pháp lựa chọn Để xác định đƣợc sự phân bố mật độ của đá móng, trƣớc hết ta xác định dị thƣờng trọng lực của các ranh giới trầm tích nằm phía trên nó mà độ sâu tới mỗi ranh giới đã đƣợc xác định bằng phƣơng các phƣơng pháp địa vật lý khác. Dị thƣờng "dƣ" đƣợc thiết lập bằng cách lấy hiệu dị thƣờng quan sát đƣợc với tổng của các dị thƣờng này tại tất cả các điểm quan sát trên bề mặt quan sát sẽ là dị thƣờng phản ánh sự bất đồng nhất của mật độ trong đá móng. Để xác định đƣợc sự phân bố mật độ này ta chia nó (đá móng) thành các lăng trụ thẳng đứng đặt cạnh nhau, mỗi lăng trụ (lăng trụ thứ (i,j)) có bề rộng . Nếu ký hiệu , i =1, 2, 3... M, j=1,3…N là mật độ dƣ của lăng trụ thứ (i,j) so với mật độ trung bình của đá móng ở lần lặp thứ k, tƣơng ứng là các hệ số của thành phần trƣờng khu vực (ở đây ta giả sử trƣờng phông khu vực là tuyến tính hàm bậc 2) thì dị thƣờng trọng lực (dị thƣờng tính toán) của móng ở lần lặp đó tại điểm quan sát thứ (I,J) trên bề mặt là: (3.1) Ký hiệu là sự sai lệch giữa dị thƣờng quan sát và dị thƣờng tính toán ở lần lặp thứ k thì: (3.2) trong đó: Trong biểu thức (3.2) x và y là tọa độ điểm quan sát thứ (I,J) đã biết, tính đƣợc từ (2.23), là độ sai lệch giữa dị thƣờng quan sát và dị thƣờng tính toán tại điểm thứ (I,J). Nhƣ vậy rõ ràng là tại M*N điểm quan sát trên diện ta sẽ xây dựng đƣợc hệ M*N phƣơng trình. Việc giải hệ M*N phƣơng trình này để tìm các giá trị đƣợc thực hiện bằng phƣơng pháp lặp thông qua việc cực tiểu hóa hàm đối tƣợng nhờ áp dụng phƣơng pháp cực tiểu hóa của Marquardt [19]. Các phƣơng trình đƣợc viết nhƣ sau: (3.3) Trong đó: I = 1 M; J=1÷N là chỉ số chạy của điểm quan sát thứ (I,J) trên bề mặt, còn i, j là chỉ số chạy lăng trụ thứ (i,j). là hằng số, nó nhận giá trị = 1 khi i = I, j=J , và = 0 khi i I, jJ; là thông số tắt dần, hay còn gọi là thông số điều chỉnh của Marquardt, giá trị đã đƣợc chứng minh là thõa mãn cho hầu hết các bài toán phức tạp với i = 2 M - 1; j=2÷N-1 chính là sự thay đổi hệ số trƣờng phông khu vực sau mỗi lần lặp là đạo hàm trƣờng theo các biến a(i,j) Quá trình lặp bao gồm 5 bƣớc sau: 1. Gán với tất cả các giá trị và các = 0 2. Tính với i = 2M -1, j=2÷N-1 sau đó giải hệ phƣơng trình (3.3) để tìm và các 3. Cộng những gia số này vào các thông số tƣơng ứng. 4. Tính theo công thức (3.1), tìm hiệu số với tất cả các giá trị I = 1M, J=1÷N. 5. Lặp lại các bƣớc 2, 3 và 4 nếu sai số bình phƣơng trung bình giữa dị thƣờng quan sát và tính toán chƣa nhỏ hơn sai số cho phép. Hoặc đơn giản hơn ta cũng có thể viết lại (3.1) dƣới dạng sau : Với ; Ở đây d là ma trận liên quan tới số liệu quan sát ( ), là ma trận chứa các tham số cần tìm một cách tối ƣu cụ thể trong bài toán này là các giá trị mật độ trong các ô, G đƣợc gọi là ma trận Jacobian mà các phần tử (i,j) của nó là các đạo hàm riêng của dị thƣờng tính toán đƣợc ở điểm thứ i theo thông số j, là ma trận chuyển vị của ma trận G. Một vấn đề thƣờng phải quan tâm tới là sự kỳ dị của ma trận trong phƣơng trình trên. Vấn đề này có thể khắc phục đƣợc, theo Marquardt, bằng cách viết phƣơng trình trên dƣới dạng : Trong đó là ma trận đơn vị. Với cách này và việc chọn thích hợp ta sẽ tránh đƣợc khả năng bằng 0 của một trong các phần tử nằm trên đƣờng chéo chính của ma trận trong phƣơng trình trên, nghĩa là tránh đƣợc kỳ dị của ma trận 3.1.2. Xác định sự phân bố mật độ đá móng theo phƣơng pháp trực tiếp Cũng giống nhƣ phƣơng pháp lựa chọn đã trình bày ở trên, theo phƣơng pháp này, để xác định sự phân bố mật độ của đá móng trƣớc hết ta cũng đi xác định dị thƣờng "dƣ" bằng cách lấy hiệu dị thƣờng quan sát đƣợc sau khi đã loại bỏ trƣờng phông khu vực với tổng của các dị thƣờng gây ra bởi các ranh giới trầm tích bên trên nó tại tất cả các điểm quan sát . Móng cũng đƣợc chia ra thành các lăng trụ thẳng đứng đặt cạnh nhau, mỗi lăng trụ (lăng trụ thứ (i,j)) có bề rộng bằng khoảng cách x, Δy giữa các điểm quan sát, có đáy trên là mặt móng và đáy dƣới là bề mặt moho và có mật độ tƣơng ứng là (i,j). Theo phƣơng pháp này quá trình tính toán đƣợc thực hiện theo các bƣớc nhƣ sau : 1. Từ các giá trị dị thƣờng quan sát (ở đây là trƣờng dƣ đá móng sau khi đã loại bỏ trƣờng phông khu vực và trƣờng do lớp trầm tích gây ra), đánh giá ban đầu về sự phân bố mật độ dƣ của đá móng đƣợc thực hiện theo phƣơng pháp xác định độ sâu trực tiếp của Bott. Theo phƣơng pháp này mật độ dƣ của mỗi lăng trụ đƣợc xác định bởi: (3.4) khi mật độ dƣ không đổi theo chiều sâu ( = 0), hoặc: khi mật độ dƣ thay đổi theo quy luật hàm mũ theo độ sâu ( 0), Trong đó: i = 1.2... M, j=1,2…N là số thứ tự các điểm quan sát trên tuyến. là dị thƣờng dƣ đá móng tại điểm quan sát thứ (i,j) : là bề dày của móng tại điểm thứ (i,j) . 2. Theo thuật toán của Murthy, Rao xác định dị thƣờng trọng lực của mỗi lăng trụ này rồi sau đó lấy tổng dị thƣờng trọng lực của M*N lăng trụ để thu đƣợc dị thƣờng của móng tại tất cả các điểm quan sát. 3. Ký hiệu là độ lệch giữa dị thƣờng và dị thƣờng tính toán của lớp tại điểm thứ (i,j) trên mặt quan sát . Độ lệch này đƣợc sử dụng để thay đổi mật độ của lớp sau mỗi lần lựa chọn: khi = 0 (3.5) hoặc khi 0 (3.6) Quá trình lựa chọn dừng lại khi sai số bình phƣơng trung bình giữa dị thƣờng "dƣ" quan sát và tính toán nhỏ hơn sai số cho phép hoặc số lần lựa chọn vƣợt quá số lần lựa chọn đã đƣợc định trƣớc. 3.2. Thuật toán và sơ đồ khối 3.2.1 Thuật toán. Chúng tôi đã xây dựng chƣơng trình giải bài toán ngƣợc xác định sự phân bố mật độ đá móng theo phƣơng pháp trực tiếp, trên cả ngôn ngữ lập trình Fortran và Matlab và bƣớc đầu chúng tôi tính toán thử nghiệm trên mô hình số để kiểm tra chƣơng trình và từ đó biết đƣợc mức độ tính đúng đắn của thuật toán. Với giả định rằng, chỉ có sự thay đổi mật độ theo diện và không có sự thay đổi mật độ theo độ sâu thì việc giải bài toán ngƣợc xác định sự phân bố mật độ của đá móng bao gồm 4 bƣớc sau : 1. Gán với tất cả các giá trị i = 1N ; j=1÷M ; 2. Giải bài toán thuận 2.23 và 2.24, i = 1 N, j=1÷M, để xác định dị thƣờng tổng các lăng kính đến từng điểm quan sát ta đƣợc tại từng điểm quan sát. 3. Tính độ lệch tại từng điểm , tính 4. Áp dụng phƣơng pháp cực tiểu hóa của Marquardt lặp lại các bƣớc 2,3 cho đến khi . Trong đó error : là sai số cho phép do ngƣời sử dụng đặt. 3.2.2. Sơ đồ khối Hình 3.1. Sơ đồ khối giải bài toán ngƣợc 3D xác định sự phân bố mật độ của đá móng Trong sơ đồ trên : lần lƣợt là trƣờng quan sát, trƣờng phông khu vực, trƣờng nâng trƣờng, trƣờng trầm tích, trƣờng đá móng, độ sâu đáy biển, độ sâu đến đáy trầm tích, độ sâu đến mặt moho và sai số do ngƣời đặt. Ở đây và là có mối tƣơng quan với nhau. Theo [8] mối tƣơng quan này đƣợc biễu diễn: (3.7) Để có đƣợc mức nâng trƣờng đƣợc chấp nhận là trƣờng phông khu vực ta thực hiện các mức nâng trƣờng khác nhau sau đó tính hệ số tƣơng quan r và lấy mức nâng trƣờng có hệ số r cao nhất làm trƣờng phông khu vực. 3.3. Tính toán thử nghiệm trên mô hình Để kiểm tra mức độ chính xác của chƣơng trình cũng nhƣ tính đúng đắn của thuật toán, dƣới đây việc giải bài toán ngƣợc xác định sự phân bố mật độ của đá móng dƣợc chúng tôi tính toán thử nghiệm trên 3 mô hình 3D có phạm vi 330 x 330 km. Khoảng cách theo cả hai chiều x và y đều là 3,3 km. Với mỗi mô hình, việc giải bài toán đƣợc thực hiện theo các bƣớc sau: 1. Giải bài toán thuận xác định các thành phần trƣờng gây ra bởi tùng ranh giới, từ đó xác định trƣờng tổng và xem nó nhƣ trƣờng quan sát. Ở đây môi trƣờng đƣợc phân chia thành 3 lớp : lớp trầm tích, lớp đá móng và lớp dƣới mặt moho. Do đó dị thƣờng quan sát đƣợc xác định bởi , với lần lƣợt là các thành phần trƣờng tƣơng ứng với từng lớp: lớp trầm tích, lớp đá móng và lớp dƣới mặt moho gây ra. 2. Tính trƣờng phông khu vực theo 2 phƣơng pháp: Phƣơng pháp xấp xỉ nó bằng đa thức bậc 3 và phƣơng pháp nâng trƣờng mà mức nâng đƣợc tối ƣu hóa bởi phƣơng pháp xấp xỉ bằng đa thức. Trƣờng móng dƣ thu đƣợc sau khi loại bỏ phần trƣờng phông này và phần trƣờng của trầm tích: đƣợc xem nhƣ dị thƣờng gây bởi sự bất đồng nhất của mật độ trong đá móng . 3. Giải bài toán ngƣợc đối với phần trƣờng dƣ này theo các phƣơng pháp kể trên, ta tìm đƣợc sự phân bố đó của mật độ trong đá móng Trong cả ba mô hình sự thay đổi mật độ dƣ của đá móng đƣợc biểu diễn trên hình vẽ: hình 3.2 Hình 3.2. Mô hình sự phân bố mật độ dƣ trong đá móng 3.3.1. Mô hình 1: 3.3.1. 1. Các tham số của mô hình a) Tham số về địa hình của các mặt ranh giới: Đối với mô hình này, địa hình mặt moho, mặt móng và bề dày trầm tích đều có dạng của một mặt cong bậc 3 (Hình 3.3) b) Tham số về mật độ - Mật độ dƣ của trầm tích suy giảm theo độ sâu theo quy luật hàm bậc hai : - Mật độ dƣ của lớp dƣới mặt Moho đƣợc lấy đồng nhất là 0,53 g/cm3 3.3.1. 2. Kết quả tính toán Kết quả tính toán đối với mô hình 1 đƣợc biểu diễn trên các hình vẽ: Hình 3.6 ÷ 3.12. Trong mỗi hình vẽ, hình bên trái biểu diễn kết quả tính toán theo phƣơng pháp tính phông khu vực bởi đa thức bậc 3 còn hình bên phải biểu diễn kết quả tính toán theo phƣơng pháp tính phông khu vực bởi phƣơng pháp nâng trƣờng. Hình 3.3. Địa hình các ranh giới và các thành phần trƣờng tƣơng ứng Hình 3.4. Trƣờng quan sát Hình 3.5. Tƣơng quan giữa trƣờng phông bậc 3 và các mức nâng trƣờng Hình 3.6. Trƣờng phông khu vực Hình 3.7. Trƣờng móng dƣ Hình 3.8. Trƣờng móng dƣ ở lần lặp cuối Hình 3.9. Sai số giữa trƣờng móng dƣ và trƣờng móng dƣ ở lần lặp cuối
(rms=0.025361 và 0.10618) Hình 3.10. Tốc độ hội tụ Hình 3.11. Kết quả tính toán sự phân bố mật độ dƣ trong đá móng Hình 3.12 : Sai số giữa sự phân bố mật độ theo mô hình và tính toán (rms=0.045527 và rms=0.038096) 3.3.1.3. Nhận xét - Việc giải bài toán ngƣợc xác định sự phân bố mật độ đá móng cho kết quả tính toán cả về hình dạng lẫn giá trị số là khá tƣơng đồng với mô hình - Sai số bình phƣơng trung bình của mật độ dƣ khi loại phông khu vực bằng cách kết hợp xấp xỉ nó bởi đa thức bậc ba với phép nâng trƣờng là nhỏ hơn đáng kể khi chỉ tính theo phông bậc 3. - Tốc độ hội tụ của việc giải bài toán ngƣợc khi loại phông khu vực bằng cách kết hợp xấp xỉ nó bởi đa thức bậc ba với phép nâng trƣờng là nhanh hơn đáng kể khi chỉ tính theo phông bậc 3. 3.3.2. Mô hình 2. 3.3.2.1. Các tham số mô hình a) Tham số về địa hình của các mặt ranh giới: Đối với mô hình này, địa hình mặt moho, mặt móng và bề dày trầm tích đều có dạng của một mặt cong bậc 3 (Hình 3.13) b) Tham số về mật độ - Mật độ dƣ của trầm tích suy giảm theo độ sâu theo quy luật hàm bậc hai : - Mật độ dƣ của lớp dƣới mặt Moho đƣợc lấy đồng nhất là 0,53 g/cm3 3.3.2.2. Kết quả tính toán Hình 3.13. Địa hình các ranh giới và các thành phần trƣờng tƣơng ứng Hình 3.14. Trƣờng quan sát Hình 3.15. Tƣơng quan giữa trƣờng phông bậc 3 và các mức nâng trƣờng Hình 3.16. Trƣờng phông khu vực Hình 3.17. Trƣờng móng dƣ Hình 3.18. Trƣờng móng dƣ ở lần lặp cuối Hình 3.19. Sai số giữa trƣờng móng dƣ và trƣờng móng dƣ ở lần lặp cuối (rms=0.023847 và rms=0.10654) Hình 3.20. Tốc độ hội tụ Hình 3.21. Kết quả tính toán sự phân bố mật độ dƣ trong đá móng Hình 3.22 : Sai số giữa sự phân bố mật độ theo mô hình và tính toán (rms=0.043996 và rms = 0.041496) 3.3.2.3. Nhận xét: - Khi địa hình mặt móng có dạng phức tạp thì sai số rms tăng lên - Tốc độ hội tụ khi lọc trƣờng phông khu vực bằng phƣơng pháp nâng trƣờng vẫn nhanh hơn khi lọc bằng phƣơng pháp xấp xỉ đa thức bậc 3. Tuy nhiên, trong thực tế địa hình các mặt ranh giới mật độ không trơn nhƣ 2 mô hình trên. Dƣới đây tác giả đƣa ra mô hình mà cả 3 mặt ranh giới đều có địa hình phức tạp. 3.3.3. Mô hình 3. 3.3.3.1. Các tham số mô hình . a) Tham số về địa hình của các mặt ranh giới: Đối với mô hình này, địa hình mặt moho, mặt móng và bề dày trầm tích đều có dạng của một mặt cong bậc 3 (Hình 3.23) b) Tham số về mật độ - Mật độ dƣ của trầm tích suy giảm theo độ sâu theo quy luật hàm bậc hai : - Mật độ dƣ của lớp dƣới mặt Moho đƣợc lấy đồng nhất là 0,53 g/cm3 3.3.3.2 Kết quả tính toán Hình 3.23: Địa hình các ranh giới và các thành phần trƣờng tƣơng ứng Hình 3.24: Trƣờng quan sát Hình 3.25: Tƣơng quan giữa trƣờng phông bậc 3 và các mức nâng trƣờng Hình 3.26 : Trƣờng phông khu vực Hình 3.27: Trƣờng móng dƣ Hình 3. 28 : Trƣờng móng dƣ ở lần lặp cuối Hình 3.29 : Sai số giữa trƣờng móng dƣ và trƣờng móng dƣ ở lần lặp cuối (rms = 0.027163 va rms=0.12242) Hình 3.30: Tốc độ hội tụ Hình 3.31: Kết quả tính toán sự phân bố mật độ dƣ trong đá móng Hình 3.32 : Sai số giữa sự phân bố mật độ theo mô hình và tính toán (rms= 0.048527 và rms = 0.050323) 3.3.3.3. Nhận xét - Tốc độ hội tụ khi giải bằng phƣơng pháp nâng trƣờng vẫn cho kết quả hội tụ nhanh hơn khi chỉ tính bằng phƣơng pháp xấp xỉ đa thức bậc 3. - Do địa hình phức tạp nên sai số rms tăng lên - Sai số khi tính bằng phƣơng pháp nâng trƣờng có gradient mạnh hơn, kết quả này thể hiện rõ sự ảnh hƣởng của địa hình mặt moho, bề dày trầm tích lên kết quả tính toán. KẾT LUẬN Qua các mô hình trên, tác giả có một vài nhận xét sau : - Thuật toán và chƣơng trình xây dựng khá đơn giản nhƣng mang lại kết quả tính khá chính xác cho tốc độ hội tụ nhanh và ổn định. - Khi giải bài toán ngƣợc xác định sự phân bố mật độ của đá móng, việc lọc
trƣờng phông cho kết quả tốt hơn khi có sự kết hợp giữa phƣơng pháp xấp xỉ nó bởi
đa thức bậc n với phƣơng pháp nâng trƣờng thông qua việc tính hệ số tƣơng quan
nhằm tìm ra mức nâng tối ƣu. - Kết quả tính trên cả ba mô hình từ đơn giản đến phức tạp cho thấy mặc dù
môi trƣờng địa chất có địa hình ranh giới phức tạp, Việc giải bài toán ngƣợc xác
định sự phân bố mật độ của đá móng vẫn cho sai số chấp nhận đƣợc (rms chỉ thay
đổi trong khoảng từ 0,038 đến 0,048 g/cm3). Điều đó cũng chứng tỏ tính ổn định
của thuật toán và chƣơng trình. Vì vậy hoàn toàn có thể áp dụng nó trong việc giải
quyết các bài toán thực tế Qua việc thực hiện đề tài mà luận văn đặt ra, chúng tôi có một số kiến nghị sau về hƣớng nghiên cứu tiếp theo của luận án: - Nghiên cứu các thuật toán giải bài toán thuận cho phép xác định nhanh hiệu
ứng trọng lực của đối tƣợng gây dị thƣờng nhằm giảm thời gian tính khi giải quyết
các bài toán 3D vốn đòi hỏi thời gian tính nhiều khi quá lớn. - Tìm hiểu thêm về các phƣơng pháp tách, lọc trƣờng để có đƣợc phần trƣờng dƣ phản ánh tốt hơn sự bất đồng nhất về mật độ trong đá móng. - Áp dụng vào thực tiễn nhằm nghiên cứu sự phân bố 3D mật độ của đá móng trong phạm vi thềm lục địa Việt nam và vùng Biển Đông kế cận TÀI LIỆU THAM KHẢO Tiếng Việt 1. Đỗ Đức Thanh (2006), Các phương pháp phân tích, xử lý tài liệu từ và trọng lực, NXB Đại học Quốc gia Hà Nội. 2. Cao Đình Triều (2000), Trọng lực và phương pháp thăm dò trọng lực, NXB Khoa học và kỹ thuật, Hà Nội. Tiếng Anh 3. Bhaskara Rao, D., Prakash, M.I., and Ramesh Babu, N.(1990), “3 and D modelling of gravity anomalies with variable density contrast”, Geophys. Prosp, Vol.38, pp. 411-422. 4. Bhaskara Rao, D., Prakash, M.I. and Ramesh Babu, N.(1993),”Gravity interpretation using Fourier transforms and simple geometrical models with exponential density contrast”, Geophysics, Vol.58, pp. 1074-1083. 5. Carlos A. Mendonca, Ahmed M.A.Meguid (2008),”Programs to compute magnetization to density ratio and the magnetization inclination from 3-D gravity and magnetic anomalies”, Computers & Geosciences, Vol.34(6), pp. 603-610. 6. Chai, Y. and Hinze, W.J., (1988),”Gravity inversion of interface above which the density contrast varies exponentially with depth”, Geophysics, Vol.53, pp. 837-845. 7. David Gómez-Ortiz Bhrigu N.P. Agarwal (2005), ”3DINVER.M: a MATLAB program to invert the gravity anomaly over a 3D horizontal density interface by Parker–Oldenburg's algorithm”, Computers & Geosciences, Volume 31(4), pp. 513–520. 8. Hualin Zeng, Deshu Xu, and Handong Tan (2007), “A model study for estimating optimum upward-continuation height for gravity separation with application to a Bougher gravity anomaly over amineral deposit, Jilin province, northeast China”, Geophysics, Vol.72(4),PP.145-150. 9. K. Mallick and K.K. Sharma (1999), “A finite element method for computation the regional gravity anomal”, Geophysics, Vol 36-07-04, PP.461-469. 10. Litinsky, V.A.,1989,”Concept of effective density : key to gravity depth determination for sedimentary basins”, Geophysics, Vol. 54, PP. 1474 - 1482. 11. P.Rama Rao, K.V. Swamy, I.V.Radhakrishna Murthy (1999),”Inversion of gravity anomalies of three-dimensional density interfaces “, Computer&geosciences, Volume 25(8), PP. 887–896. 12. Ramakrishna Murthy, I. V., Rama Rao, P., and Jagannadha Rao, S. (1990),”The density difference and generalized programs for two-and three-dimensionals gravity modeling”, Computer&geosciences, Vol.16(3), PP.277-287. 13. Richard J, Blakely (1996), Potential theory in gravity and magnetic application, Cambridge University Press. 14. R. Nagendra, P.V.S. Prasad, V.L.S. Bhimasankaram (1996), “Forward and inverse computer modeling of a gravity field resulting from a density interface using Parker-Oldenberg method”, Computers & Geosciences, Volume 22(3),PP.227–237. 15. Ya Xu, Tianyao Hao, Zhiwei Li, Qiuliang Duan and Lili Zhang (2009), “Regional gravity anomaly separation using wavelet transform and spectrum”, Journal of Geophysics and engineering, PP. 279-287.15
16
17
18
19
20
21
22
23
24
là hiệu ứng trọng lực của lăng trụ nằm dƣới điểm
25
26
27
28
29
30
31
32
33
34
35
36
a) b)
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
