ĐẠI HỌC QUỐC GIA HÀ NỘI
TRƢỜNG ĐẠI HỌC KHOA HỌC TỰ NHIÊN
---------------------
PHẠM THỊ THANH HOA
NGHIÊN CỨU ÁP DỤNG PHƢƠNG PHÁP BIẾN ĐỔI TRƢỜNG
TRONG MIỀN TẦN SỐ THỰC HIỆN VIỆC BIẾN ĐỔI TRƢỜNG
TRỌNG LỰC KHU VỰC BỂ TRẦM TÍCH SÔNG HỒNG
LUẬN VĂN THẠC SĨ KHOA HỌC
Hà Nội – 2015
ĐẠI HỌC QUỐC GIA HÀ NỘI
TRƢỜNG ĐẠI HỌC KHOA HỌC TỰ NHIÊN
---------------------
PHẠM THỊ THANH HOA
NGHIÊN CỨU ÁP DỤNG PHƢƠNG PHÁP BIẾN ĐỔI TRƢỜNG
TRONG MIỀN TẦN SỐ THỰC HIỆN VIỆC BIẾN ĐỔI TRƢỜNG
TRỌNG LỰC KHU VỰC BỂ TRẦM TÍCH SÔNG HỒNG
Chuyên ngành: Vật lý địa cầu
Mã số: 60440111
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 - 2015
MỤC LỤC
DANH MỤC HÌNH VẼ
MỞ ĐẦU ..................................................................................................................... 1
CHƢƠNG 1: CÁC PHƢƠNG PHÁP BIẾN ĐỔI TRƢỜNG TRỌNG LỰC ........... 2
1.1. NHỮNG NGUYÊN LÝ CHUNG VỀ PHÉP BIẾN ĐỔI TRƢỜNG TRỌNG
LỰC. ............................................................................................................................ 2
1.2. CÁC PHƢƠNG PHÁP BIẾN ĐỔI TRƢỜNG .................................................... 3
1.2.1. Phƣơng pháp trung bình hoá ............................................................................. 3
1.2.2. Phƣơng pháp tiếp tục giải tích trƣờng ............................................................... 6
1.2.3. Phƣơng pháp tính đạo hàm bậc cao ................................................................ 12
CHƢƠNG 2: CƠ SỞ DỮ LIỆU VÀ KẾT QUẢ BIẾN ĐỔI TRƢỜNG TRỌNG
LỰC KHU VỰC BỂ TRẦM TÍCH SÔNG HỒNG .................................................. 15
2.1. NGUỒN SỐ LIỆU SỬ DỤNG ......................................................................... 15
2.2. SƠ LƢỢC VỀ ĐẶC ĐIỂM ĐỊA CHẤT - KIẾN TẠO KHU VỰC BỂ TRẦM
TÍCH SÔNG HỒNG ................................................. Error! Bookmark not defined.
2.3. ĐẶC ĐIỂM DỊ THƢỜNG TRỌNG LỰC ………………………………….. 18
2.4. KẾT QUẢ TÍNH TOÁN .................................................................................... 19
KẾT LUẬN ............................................................................................................... 31
TÀI LIỆU THAM KHẢO
PHỤ LỤC
DANH MỤC HÌNH VẼ
Hình 1.1. Cách chọn các bán kính của palet ............................................................... 5
Hình 1.2. Đánh giá hàm điều hoà tại một điểm bất kỳ trong vùng R ......................... 7
Hình 2.1. Bản đồ dị thƣờng Phai khu vực Biển Đông và kế cận .............................. 15
Hình 2.2. Bản đồ dị thƣờng Bughe khu vực Biển Đông và kế cận ........................... 16
Hình 2.3. Bản đồ dị thƣờng Bughe khu vực bể trầm tích Sông Hồng và kế cận ...... 17
Hình 2.4. Kết quả hạ trƣờng xuống độ sâu 1km ....................................................... 20
Hình 2.5. Kết quả hạ trƣờng xuống độ sâu 2 km ...................................................... 21
Hình 2.6. Kết quả tính đạo hàm ngang toàn phần ở mức z = 0................................. 22
Hình 2.7. Kết quả nâng trƣờng lên độ cao 5 km ....................................................... 23
Hình 2.8. Kết quả tính đạo hàm ngang toàn phần khi nâng trƣờng lên độ cao 5 km...
................................................................................................................................... 24
Hình 2.9. Kết quả nâng trƣờng lên độ cao 10 km ..................................................... 25
Hình 2.10. Kết quả tính đạo hàm ngang toàn phần khi nâng trƣờng lên độ cao 10 km
................................................................................................................................... 26
Hình 2.11. Kết quả nâng trƣờng lên độ cao 15 km ................................................. 257
Hình 2.12. Kết quả tính đạo hàm ngang toàn phần khi nâng trƣờng lên độ cao 15 km
................................................................................................................................... 28
Hình 2.13. Kết quả nâng trƣờng lên độ cao 20 km ................................................... 29
Hình 2.14. Kết quả tính đạo hàm ngang toàn phần khi nâng trƣờng lên độ cao 20 km
................................................................................................................................... 30
LỜI CẢM ƠN
Trong suốt quá trình học tập tại trƣờng Đại học khoa học tự nhiên em đã
nhận đƣợc sự tận tình dạy dỗ, chỉ bảo của các thầy cô trong khoa Vật Lý nói riêng
và các thầy cô trong trƣờng nói chung. Em xin gửi lời cảm ơn tới toàn thể các thầy
cô giáo đã dạy em trong thời gian qua.
Em xin chân thành cảm ơn các thầy cô giáo trong bộ môn Vật Lý Địa Cầu
đã trang bị cho em những kiến thức cơ bản trong thời gian học tập tại trƣờng. Và
đặc biệt, em xin gửi lời cảm ơn sâu sắc tới PGS.TS. Đỗ Đức Thanh, ngƣời đã trực
tiếp hƣớng dẫn em hoàn thành tốt luận văn này.
Cuối cùng em xin chân thành cảm ơn gia đình và các bạn đã quan tâm động
viên và giúp đỡ em trong quá trình học tập và trong thời gian làm luận văn.
Em mong nhận đƣợc sự quan tâm và góp ý của thầy cô và các bạn về luận
văn này.
Em xin chân thành cảm ơn!
Hà Nội, tháng 6 năm 2015
Học viên
Phạm Thị Thanh Hoa
MỞ ĐẦU
Thăm dò trọng lực là phƣơng pháp địa vật lý nghiên cứu sự phân bố
trƣờng trọng lực trên mặt đất để nghiên cứu cấu trúc bên trong của quả đất, cấu
trúc vỏ của quả đất, tìm kiếm thăm dò khoáng sản và giải quyết các nhiệm vụ địa
chất khác nhau.
Trƣờng trọng lực quan sát đƣợc là tổng hợp nhiều nguồn trƣờng của các đối
tƣợng địa chất khác nhau vì thế đặc điểm của chúng rất phức tạp và đa dạng. Để xác
định dị thƣờng trọng lực liên quan tới đối tƣợng cần nghiên cứu thì một nhiệm vụ
quan trọng là phải áp dụng các thuật toán để phân chia trƣờng trọng lực thành các
trƣờng thành phần (trƣờng khu vực, trƣờng địa phƣơng…), tách biệt các trƣờng liên
quan đến các đối tƣợng cụ thể (nâng và hạ trƣờng, trung bình trƣờng, gradien chuẩn
hoá) và nhận dạng trƣờng…
Trong luận văn này, tôi áp dụng phƣơng pháp biến trƣờng trọng lực trong
miền tần số để thực hiện việc biến đổi trƣờng bao gồm việc nâng, hạ trƣờng ở các
mức khác nhau. Đồng thời thực hiện việc tính đạo hàm bậc cao theo phƣơng nằm
ngang của thế trọng lực ở các độ cao khác nhau cũng đƣợc thực hiện nhằm xác định
vị trí của các đứt gãy sâu trong khu vực bể trầm tích Sông Hồng thuộc phạm vi
thềm lục địa Việt nam.
Luận văn đƣợc chia làm 2 chƣơng sau:
Chƣơng 1: Các phƣơng pháp biến đổi trƣờng trọng lực
Chƣơng 2: Cơ sở dữ liệu và kết quả biến đổi trƣờng trọng lực khu vực bể
trầm tích Sông Hồng
1
CHƢƠNG 1
CÁC PHƢƠNG PHÁP BIẾN ĐỔI
TRƢỜNG TRỌNG LỰC
1.1. NHỮNG NGUYÊN LÝ CHUNG VỀ PHÉP BIẾN ĐỔI TRƢỜNG TRỌNG
LỰC.
Các dị thƣờng trọng lực quan sát đƣợc, phản ánh toàn bộ các yếu tố địa chất.
Trong trƣờng tổng cộng, mỗi một yếu tố địa chất đều có đóng góp một phần nhất
định. Trong khi giải quyết các nhiệm vụ địa chất cụ thể, từ trƣờng tổng ngƣời ta
phải tách ra đƣợc các phần trƣờng riêng biệt có liên hệ trực tiếp đến đối tƣợng cần
nghiên cứu. Muốn vậy ngƣời ta phải tiến hành biến đổi trƣờng quan sát đƣợc nhằm
nhấn mạnh phần trƣờng cần thiết và làm yếu đi các phần trƣờng khác.
Các phƣơng pháp biến đổi trƣờng dị thƣờng trọng lực có nhiều điểm chung
với quá trình lọc nhiễu trong lý thuyết thông tin, mặc dù chúng có những đặc điểm
riêng của mình.
Phần cơ bản của phép biến đổi trƣờng trọng lực bao gồm việc tách trƣờng
quan sát ra thành các thành phần tƣơng ứng với các đối tƣợng địa chất nằm ở các độ
sâu khác nhau.
Các phép biến đổi trƣờng chỉ nhấn mạnh phần này và làm yếu phần khác các
thành phần có các đặc điểm khác nhau nằm trong trƣờng tổng .
Hiện nay có rất nhiều phƣơng pháp tính trong dị thƣờng trọng lực [1, 2]. Phụ
thuộc vào phép biến đổi mà hàm sau khi đƣợc biến đổi có thể là hàm thứ nguyên
của hàm xuất phát (nhƣng thuộc về mức khác) hoặc là các đạo hàm của hàm xuất
phát. Các đạo hàm này có thể thuộc mức xuất phát. Các hàm đã đƣợc biến đổi đôi
khi có thứ nguyên là tích của hàm xuất phát với toạ độ. Tất cả các phép biến đổi
trƣờng trọng lực cũng nhƣ phƣơng pháp lọc nhiễu trong lý thuyết thông tin về mặt
toán học đƣợc biểu diễn dƣới dạng tích phân chập:
- Trong trƣờng hợp bài toán ba chiều:
(1.1) Vbđ(x0,y0,z0)=
- Trong trƣờng hợp bài toán hai chiều:
2
(1.2) Vbđ(x0,y0)=
trong đó Vbđ(x0,y0,z0) và Vbđ(x0,z0) là các tham số đã đƣợc biến đổi, còn
là các hàm xuất phát (trƣờng tổng), K( và và Vxp Vxp(
là các nhân biến đổi. Các nhân biến đổi nhiều khi còn đƣợc gọi là các K
hàm trọng số. Các hàm này gọi là các hàm tuyến tính nên tất cả các biến đổi tƣơng
ứng đƣợc gọi là các biến đổi tuyến tính.
Bằng cách qui ƣớc ngƣời ta có thể chia các phép biến đổi trƣờng ra thành ba
nhóm lớn:
- Trung bình hoá.
- Tiếp tục giải tích các dị thƣờng trọng lực nhƣ là các hàm điều hoà.
- Tính các đạo hàm bậc cao của thế trọng lực.
Ta sẽ lần lƣợt xét đến các nhóm phƣơng pháp này.
1.2. CÁC PHƢƠNG PHÁP BIẾN ĐỔI TRƢỜNG
1.2.1. Phƣơng pháp trung bình hoá [2]
Việc phân chia các dị thƣờng trọng lực ra thành các thành phần khu vực và
địa phƣơng nhờ phƣơng pháp trung bình hoá đƣợc sử dụng rộng rãi trong thực tế.
Bản chất của phƣơng pháp trung bình hoá nhƣ sau: Xem trƣờng trọng lực quan sát
đƣợc gồm hai thành phần, thành phần khu vực Vr và thành phần địa phƣơng Vl.
V=Vr+Vl. (1.3)
Lấy trung bình thƣờng quan sát đƣợc trong phạm vi của vòng tròn bán kính
R. Giá trị trung bình đó đƣợc biểu diễn bằng tích phân sau:
(1.4)
Đối chiếu các công thức tổng quát (1.1) ta thấy rằng trong trƣờng hợp này:
Vbđ(x0,y0,z0)=
.
Còn Vxp(
3
Ngƣời ta chọn bán kính R sao cho nó lớn hơn nhiều so với các kích thƣớc
của các dị thƣờng địa phƣơng và nhỏ hơn nhiều so với kích thƣớc của dị thƣờng
khu vực. Khi điều kiện đó đƣợc thoả mãn thành phần khu vực đƣợc tách ra từ
trƣờng quan sát và các dị thƣờng địa phƣơng (dƣơng và âm) bù trừ lẫn nhau, còn
các thành phần khu vực lại ít bị thay đổi. Kết quả là . Trƣờng hợp đặc biệt
nếu trƣờng khu vực thay đổi theo quy luật tuyến tính thì nó hoàn toàn không bị thay
đổi sau phép trung bình, tức là:
Sau khi xác định đƣợc trƣờng khu vực Vr , trƣờng dị thƣờng địa phƣơng
đƣợc tính theo công thức:
(1.5)
Vl = V-
Để làm sáng tỏ ý nghĩa vật lý của dị thƣờng đã đƣợc trung bình hoá ngƣời ta
đƣa vào khái niệm mức độ trung bình hoá. Đó là tỷ số giữa trƣờng đã đƣợc trung
bình hoá với trƣờng xuất phát
(1.6)
Mức độ trung bình hoá đồng thời đặc trƣng cho mức độ chính xác của việc
tách trƣờng địa phƣơng.
Giả sử trƣờng dị thƣờng Vz do hình cầu có khối lƣợng M nằm ở độ sâu h gây
ra. (1.7)
Ta hãy tìm giá trị trung bình trong phạm vi vòng tròn bán kính R với tâm
trùng với hình chiếu của tâm quả cầu trên mặt đất.
(1.8)
Đặt . Lúc đó công thức (1.7) tƣơng tự với công thức của lực hấp
dẫn của đĩa vật chất tròn với mật độ mặt nằm ở độ sâu h có khối lƣợng bằng
khối lƣợng của quả cầu M và bán kính bằng bán kính của vòng tròn lấy trung bình
gây ra. Nhƣ vậy là phƣơng pháp trung bình hóa dị thƣờng Vz của chất điểm nằm ở
4
độ sâu h trong vòng tròn bán kính R tƣơng đƣơng với việc phân phối lại khối đó
thành đĩa vật chất nằm ở cùng độ sâu có bán kính bằng bán kính trung bình hoá
(1.9)
gọi R/h=Rh, từ công thức (1.8) ta có
(1.10)
Từ công thức (1.9) ta có thể chọn đƣợc bán kính trung bình hoá khi cho trƣớc
mức độ chính xác xác định dị thƣờng địa phƣơng và khi biết trƣớc độ sâu h.
Hình 1.1a Hình 1.1b
Hình 1.1. Cách chọn các bán kính của palet
Trong thực tế phần lớn bán kính trung bình hóa đƣợc chọn bằng phƣơng
pháp thực nghiệm theo cách lấy trung bình trƣờng trọng lực cho trƣớc. Muốn vậy
tại một số điểm khác nhau của trƣờng ta áp dụng phƣơng pháp trung bình hoá với
các bán kính trung bình khác nhau. Tiếp theo ngƣời ta vẽ đồ thị biểu diễn sự phụ
thuộc giữa trƣờng đã đƣợc trung bình hoá và bán kính trung bình. Theo đồ thị ngƣời
ta chọn các bán kính trung bình hoá tối ƣu ( H1.1.a.b.).
Trong trƣờng hợp (a) bán kính tối ƣu tƣơng ứng với điểm mà bắt đầu từ đó
trƣờng đƣợc trung bình hoá bắt đầu không thay đổi khi R thay đổi, còn trong trƣờng
hợp (b) bán kính tối ƣu tƣơng ứng với điểm uốn của đƣờng cong.
Ngoài vòng tròn, trong phƣơng pháp trung bình hoá ngƣời ta còn lấy các
hình khác nhau để làm miền trung bình. Đặc biệt trong thực tế ngƣời ta dùng miền
5
có dạng hình vuông (Palet vuông). Nhờ có palet vuông mà khối lƣợng tính toán
đƣợc giảm đi nhiều. Trên bản đồ các đƣờng đẳng trị ngƣời ta vẽ các đƣờng thẳng
đứng và nằm ngang cách đều nhau. Các đƣờng này tạo thành mạng lƣới các ô
vuông. Tại mỗi điểm nút (giao điểm của đƣờng thẳng đứng và nằm ngang) từ bản
đồ các đƣờng đẳng trị ngƣời ta nội suy giá trị trọng lực. Các giá trị này đƣợc sử
dụng liên hoàn trong quá trình tính toán trung bình hoá theo miền vuông.
Bản đồ dị thƣờng khu vực vẽ theo kết quả của phƣơng pháp trung bình
hoá đƣợc sử dụng rộng rãi trong phân vùng kiến tạo và giải quyết các nhiệm vụ
địa chất khác.
1.2.2. Phƣơng pháp tiếp tục giải tích trƣờng
1.2.2.1. Tiếp tục giải tích trường lên nửa không gian trên
Tiếp tục giải tích trƣờng lên nửa không gian trên là phép biến đổi trƣờng thế
đo đƣợc trên một mặt nào đấy thành trƣờng thế ở một mặt khác xa các nguồn
hơn. Nhƣ ta đã biết, phép biến đổi này làm suy yếu các dị thƣờng, tuỳ theo bƣớc
sóng của chúng. Dị thƣờng có bƣớc sóng càng ngắn càng bị suy yếu mạnh. Theo
nghĩa này, quá trình tiếp tục trƣờng lên trên là một quá trình làm suy biến các số
liệu đo đƣợc. Chúng có một ứng dụng rất to lớn trong thực tế. Thật vậy, khi phải
so sánh hoặc thống nhất các tài liệu từ hàng không đo đƣợc ở các độ cao khác
nhau, việc tiếp tục giải tích trƣờng lên trên cho phép biến đổi các đo đạc riêng rẽ
về một mặt phù hợp. Hơn nữa, nó có xu hƣớng làm nhấn mạnh phần dị thƣờng
gây bởi các nguồn sâu và làm yếu đi phần dị thƣờng gây bởi các nguồn nông. Ví
dụ, trong tài liệu đo đạc từ ở những vùng có đất đá phun trào trẻ, phần dị thƣờng
có bƣớc sóng ngắn liên quan tới các đá phun trào gần bề mặt luôn chiếm ƣu thế;
khi đó việc tiếp tục giải tích trƣờng lên trên sẽ làm yếu đi phần dị thƣờng này để
làm nổi rõ phần dị thƣờng gây bởi các nguồn nằm sâu hơn mà ta quan tâm, đó là
các lớp đá nằm phía dƣới.
Cơ sở của phép tiếp tục giải tích trƣờng lên trên là đẳng thức thứ ba của
Green. Theo đẳng thức này, nếu hàm U là điều hoà, liên tục và có đạo hàm liên tục
6
trên một vùng có biên đều đặn R, thì giá trị của U tại điểm P nằm phía trong R đƣợc
cho bởi phƣơng trình:
U(P)= (1.11)
Trong đó S là biên của R, n là hƣớng pháp tuyến ngoài còn R là khoảng cách
từ P tới điểm tích phân trên S (Hình 1.2). Phƣơng trình (1.11) minh hoạ nguyên lý
cơ bản của việc tiếp tục lên trên: một trƣờng thế có thể tính đƣợc tại một điểm bất
kỳ trong vùng theo các giá trị của trƣờng trên một bề mặt bao quanh vùng nó. Ở đây
không cần điều kiện gì về nguồn của trƣờng trừ yêu cầu là nó không đƣợc có mặt
trong vùng R.
1.2.2.1.1. Các biến đổi trong miền không gian [2]
Dạng tiếp tục giải tích trƣờng đơn giản nhất là dạng tiếp tục giải tích đƣợc
thực hiện đối với các trƣờng thế đo đƣợc trên một mặt mức phẳng nào đó. Trong hệ
toạ độ vuông góc với trục z hƣớng xuống dƣới, ta giả sử rằng trƣờng thế đo đƣợc ở
mặt mức z=z0 và vấn đề đặt ra là cần phải xác định trƣờng tại một điểm P(x,y,z0-∆z)
nằm phía trên mặt mức này (∆z>0) mặt S gồm cả mặt mức cộng với nửa hình cầu
bán kính α, nhƣ đƣợc chỉ ra trên Hình 1.2 giả thiết tất cả các nguồn nằm ở z>z0.
Khi cho α trở nên rất lớn ta dễ dàng chỉ ra rằng phần tích phân phƣơng trình
(1.11) trên nửa mặt cầu trở thành rất nhỏ. Vì vậy khi α→∞ thì:
Hình 1.2. Đánh giá hàm điều hoà tại một điểm bất kỳ trong vùng R
7
(1.12) U(x,y,z0-∆z)=
Trong đó:
r =
với lƣu ý rằng ở đây ∆z>0.
Theo phƣơng trình (1.12) để xác định đƣợc trƣờng tại điểm P, ta cần biết
không chỉ các giá trị của U trên bề mặt mà còn phải biết cả các giá trị gradient thẳng
đứng của U, một tổ hợp mà không phải bao giờ cũng có thể đáp ứng đƣợc trong
thực tế. Vì vậy, ta cần một cách nào đó để loại trừ số hạng đạo hàm trong phƣơng
trình (1.12). Điều này có thể thực hiện đƣợc nhờ đẳng thức Green thứ hai. Nếu V là
một hàm khác cũng điều hoà trên R, thì đẳng thức thứ hai của Green cho phép viết:
Cộng thêm kết quả này vào phƣơng trình (1.11) ta có:
U(P)= (1.13)
Để loại bỏ số hạng thứ nhất của hàm dƣới dấu tích phân, hàm điều hoà V cần
chọn sao cho V+1/r=0 tại mọi điểm trên S. Muốn vậy, ta lấy P’ là ảnh phản chiếu
gƣơng của P tại (x,y,z0+∆z) và đặt V= −1/ρ, trong đó:
ρ =
Chú ý rằng V xác định theo cách này thoả mãn các điều kiện cần thiết V+1/r
= 0 trên mặt phẳng nằm ngang còn V+1/r sẽ triệt tiêu trên nửa mặt cầu khi α trở
nên rất lớn và V luôn luôn điều hoà vì ρ không bao giờ triệt tiêu. Vậy phƣơng trình
(1.13) trở thành:
8
U(P) =
Khi nửa mặt cầu trở nên rất lớn, số hạng thứ nhất triệt tiêu tại mỗi điểm trên
S còn số hạng thứ hai triệt tiêu trừ những điểm thuộc mặt phẳng nằm ngang.
U(x,y,z0-∆z) = -
Thực hiện đạo hàm để z’ di chuyển tới mặt phẳng nằm ngang ta đƣợc:
, Δz>0 (1.14) U(x,y,z0-∆z) =
Phƣơng trình (1.14) biểu diễn sự tiếp tục lên trên của trƣờng thế. Phƣơng
trình chỉ ra cách làm thế nào để tính đƣợc giá trị của trƣờng thế tại một điểm bất kỳ
trên mặt mức nằm ngang ở độ cao ∆z từ các giá trị đã biết của trƣờng trên bề mặt
quan sát (z=0). Dĩ nhiên trong các ứng dụng thực tế một số điều kiện gần đúng cần
phải đƣợc chấp nhận vì ta không thể biết một cách chính xác trƣờng thế tại mọi
điểm trên một mặt phẳng vô hạn.
1.2.2.1.2. Các biến đổi trong miền tần số
Phƣơng trình (1.14) có thể đƣợc dùng để tiếp tục giải tích tài liệu đo đƣợc
trên một mặt mức tới mặt mức khác. Với việc áp dụng phƣơng pháp này, đối với
mỗi điểm của mặt mức mới ta đều phải thực hiện một cách có hiệu quả hơn nếu việc
tính toán đƣợc thực hiện trong miền tần số.
Chú ý rằng phƣơng trình (1.14) đơn giản là một tích phân chập hai chiều:
U(x,y,z0-∆z) =
Trong đó:
(1.15)
9
Nếu trƣờng thế U đo đƣợc trên miền z=z0 trong phạm vi đủ rộng so với kích
thƣớc của nguồn, thì tồn tại biến đổi Fourier F của nó. Biểu diễn trong miền tần
số của (1.14) tìm đƣợc bằng cách biến đổi cả hai vế của phƣơng trình (1.14) qua
miền tần số và áp dụng lý thuyết tích chập:
(1.16) F
với F là biến đổi Fourier của trƣờng đã tiếp tục lên trên. Ở đây, điều cần
thiết là tìm biểu diễn giải tích của F . Chú ý rằng
(1.17)
trong đó: r = và biết rằng
F ,z’>z0 ,
khi đó biến đổi Fourier của (1.17) là
(Δz>0) (1.18) F
) (
Vì vậy, việc tiếp tục trƣờng thế từ mức này sang một mức khác có thể thực
hiện đƣợc bằng cách biến đổi Fourier tài liệu đã đo đƣợc, nhân với số hạng hàm mũ
của phƣơng trình (1.18) rồi sau đấy biến đổi Fourier ngƣợc tích số vừa thu đƣợc.
Từ phƣơng trình (1.18) ta thấy rằng quá trình tiếp tục lên trên làm yếu dần tất
cả các số sóng trừ |k=0|. Sóng có bƣớc sóng càng ngắn càng bị làm yếu càng nhiều
và mức độ làm suy yếu cũng tăng theo gia số ∆z. Phƣơng trình (1.18) là một hàm
thực, không có thành phần pha, và do đó không có sự thay đổi pha đối với trƣờng
đƣợc tiếp tục lên trên.
Hàm U đƣợc miêu tả trong các công thức trên là một hàm thế bất kỳ nên các
phƣơng trình (1.14) và (1.16) có thể áp dụng đƣợc cho mọi thành phần của trƣờng
trọng lực cũng nhƣ trƣờng từ đo đƣợc trên một mặt nằm ngang. Nó cũng áp dụng
đƣợc cho cả dị thƣờng từ ∆T.
1.2.2.2. Tiếp tục giải tích trường xuống nửa không gian dưới
10
Tất cả các lập luận trƣớc đây đều dựa trên cơ sở thừa nhận rằng tất cả các
nguồn gây dị thƣờng đều định xứ phía dƣới mặt quan sát còn tất cả các điểm mà ta
cần tiếp tục giải tích tới đều ở bên trên mặt quan sát, tức là, việc tiếp tục là theo
hƣớng đi ra xa từ phía các nguồn. Dƣờng nhƣ sẽ là hợp logic nếu ta thử tiếp tục giải
tích tài liệu đo đƣợc vào vùng gần các nguồn hơn, dĩ nhiên, với điều kiện là thực sự
không có nguồn tồn tại trong vùng cần tiếp tục. Việc tính toán này gọi là tiếp tục
giải tích trƣờng xuống dƣới, sẽ rất hữu ích trong việc minh giải các tài liệu trọng lực
và từ vì nó có xu hƣớng làm nổi bật các chi tiết của phân bố nguồn, đặc biệt là các
đối tƣợng nằm nông.
Tuy nhiên quá trình tiếp tục xuống dƣới là một quá trình không ổn định.
Trong khi việc tiếp tục lên trên là một toán tử làm trơn, điều này dễ dàng thấy đƣợc
qua phƣơng trình (1.14) trong đó U(x,y,z0-Δz) tại điểm bất kỳ là trung bình có trọng
số tất cả các giá trị của U(x,y,z0), thì tiếp tục xuống dƣới là việc tính các giá trị
U(x,y,z0) từ U(x,y,z0-Δz),quá trình ngƣợc với (1.14) nên đây là toán tử “không làm
trơn”, và nhƣ ta đã biết, các tính toán nhƣ vậy không ổn định. Những thay đổi nhỏ
của của U(x,y,z0-Δz) có thể gây ra những biến đổi lớn và không thực trong các
U(x,y,z0) tính toán đƣợc. Điều này đƣợc chỉ ra bằng cách viết nghịch đảo phƣơng
trình (1.16)
F =F F = F e+
Trong trƣờng hợp này, F là biến đổi Fourier của trƣờng quan sát đƣợc
còn F là trƣờng muốn tìm đƣợc tiếp xúc xuống dƣới một khoảng Δz. Rõ ràng là
các thành phần của dị thƣờng có bƣớc sóng càng ngắn trong tài liệu đo đƣợc sẽ bị
khuếch đại càng mạnh trong quá trình này tới một mức nào đấy phụ thuộcvào giá trị
Δz và khoảng lấy mẫu số liệu (khoảng cách giữa các điểm quan sát ). Các sai số
ngẫu nhiên có mặt và không đƣợc phát hiện ra trong số liệu đo đạc có thể làm xuất
hiện trong trƣờng đã đƣợc tính toán những biến thiên lớn và không thực. Tuy phức
tạp nhƣ vậy nhƣng với những lợi thế riêng của mình, việc tiếp tục giải tích trƣờng
xuống dƣới vẫn đƣợc sử dụng rộng rãi trong thực tế.
11
1.2.3. Phƣơng pháp tính đạo hàm bậc cao
Xét một đại lƣợng vô hƣớng biến đổi trơn đo đƣợc trên một mặt nằm
ngang. Các đạo hàm ngang của dễ dàng đƣợc đánh giá bằng việc sử dụng
phƣơng pháp sai phân hữu hạn và các giá trị đo đƣợc của . Nếu các giá trị ,
với i=1,2………;j=1,2…….., là các giá trị đo đƣợc của trên một lƣới đều
đặn với bƣớc Δx và Δy tại điểm (i,j) đƣợc xấp xỉ bởi:
Các đạo hàm ngang cũng dễ dàng đƣợc thực hiện trong miền tần số. Theo lý
thuyết sai phân, các đạo hàm ngang của đƣợc xác định nhƣ sau:
(1.19) F
(1.20) F
Vì vậy, (ikx)n và (iky)n là các bộ lọc biến đổi một hàm đo đƣợc trên mặt nằm
ngang thành các đạo hàm đối với x và y một cách tƣơng ứng.
Nếu là một hàm thế, thì ta cũng có thể tính đƣợc các Gradient thẳng
đứng. Thực vậy các đạo hàm thẳng đứng bậc hai là hệ quả trực tiếp của phƣơng
trình Laplace, vì nếu là hàm thế, thì =0, tức là:
(1.21)
Nếu đo đƣợc trên một mặt nằm ngang thì phƣơng trình Laplace có thể
đƣợc biến đổi vào miền tần số nhờ (1.19) và (1.20), tức là:
(1.22) F
Vì vậy đạo hàm thẳng đứng bậc hai của trƣờng thế đo đƣợc trên một mặt
nằm ngang đƣợc xác định nhƣ là một toán tử lọc ba bƣớc: biến đổi Fourier trƣờng
12
thế, nhân với rồi biến đổi Fourier ngƣợc tích số vừa thu đƣợc.Vì vậy đạo hàm
thẳng đứng bậc hai của trƣờng thế đo đƣợc trên một mặt nằm ngang đƣợc xác định
nhƣ là một toán tử lọc ba bƣớc: biến đổi Fourier trƣờng thế, nhân với rồi biến
đổi Fourier ngƣợc tích số vừa thu đƣợc.
Phƣơng pháp tính đạo hàm bậc hai theo phƣơng thẳng đứng z là chỗ dựa đầu
tiên của kỹ thuật minh giải các số liệu đo đạc từ và trọng lực, vì nó là một phƣơng
pháp đơn giản nhƣng rất có hiệu quả trong việc giúp ta định vị và làm nổi bật các
nguồn nông. Để thấy đƣợc tại sao lại nhƣ vậy, hãy xét hai đơn cực từ, một nằm khá
nông ở độ sâu d1 và một nằm ở độ sâu lớn hơn d2. Trƣờng của mỗi đơn cực tại điểm
quan sát P tỷ lệ ngƣợc với bình phƣơng khoảng cách. Vì vậy, khi P lại gần các đơn
cực, trƣờng gây bởi đơn cực nông sẽ tăng nhanh hơn trƣờng tạo bởi đơn cực sâu.
Đạo hàm bậc hai cũng có cùng hiệu ứng nhƣ vậy. Ngoài ra, đạo hàm thẳng đứng
bậc hai còn là nổi rõ các biên của các nguồn từ và trọng lực.
Các đặc trƣng này của đạo hàm bậc hai cũng có thể đƣợc suy ra từ phƣơng
trình (1.22). Thật vậy, theo phƣơng trình này khi nhân trƣờng thế với thì rõ ràng
các thành phần có bƣớc sóng ngắn của trƣờng sẽ bị khuếch đại trong khi các thành
phần có bƣớc sóng dài sẽ bị làm mờ đi.
Đạo hàm bậc hai thẳng đứng đƣợc suy ra trực tiếp từ phƣơng trình Laplace
còn các đạo hàm thẳng đứng bậc bất kỳ cũng có thể thu đƣợc từ một trƣờng thế.
Điều này suy ra từ lập luận trƣớc đây về việc tiếp tục trƣờng lên trên. Dùng các qui
ƣớc thông thƣờng , z hƣớng xuống phía dƣới và với , đạo hàm thẳng bậc nhất
đƣợc cho bởi:
việc biến đổi sang miền tần số tạo ra:
F
13
Tƣơng tự ta có thể chỉ ra rằng gradient thẳng đứng bậc n bằng biến đổi
Fourier của thế nhân với hoặc tổng quát ta có:
(1.23) F
14
CHƢƠNG 2
CƠ SỞ DỮ LIỆU VÀ KẾT QUẢ BIẾN ĐỔI TRƢỜNG TRỌNG LỰC
KHU VỰC BỂ TRẦM TÍCH SÔNG HỒNG
Trên cơ sở lý thuyết về các phƣơng pháp biến đổi trƣờng trọng lực đã
trình bày ở chương 1, trong phần này tôi thực hiện lập chƣơng trình trên máy tính
để biến đổi trƣờng trọng lực bao gồm việc nâng, hạ trƣờng và tính đạo hàm bậc cao
theo phƣơng nằm ngang của thế trọng lực khu vực bể trầm tích Sông Hồng tại các
độ cao khác nhau, nhằm phát hiện các đứt gãy địa chất sâu trong phạm vi khu vực.
Chƣơng trình đƣợc viết bằng ngôn ngữ Matlab.
2.1. NGUỒN SỐ LIỆU SỬ DỤNG
Số liệu đƣợc chúng tôi sử dụng để tính toán đƣợc lấy từ nguồn số liệu dị thƣờng
trọng lực vệ tinh (free air anomaly) tỉ lệ 1 điểm đo/1 phút, vesion V.201 từ địa chỉ http://topex.ucsd.edu/cgi-bin/get_data.cgi, trải rộng trong phạm vi từ kinh tuyến 900 E đến 1300 E, từ vĩ tuyến -100 N đến 300 N.
Hình 2.1. Bản đồ dị thƣờng Phai khu vực Biển Đông và kế cận
15
Từ nguồn số liệu ban đầu này chúng tôi tiến hành lƣới hóa lại với khoảng
mắt lƣới cỡ 4km trong hệ tọa độ Đề các rồi sau đó thực hiện việc tính hiệu chỉnh
lớp giữa với mật độ lớp trung gian để xây dựng nên tờ bản đồ dị
thƣờng Bughe khu vực Biển Đông (Hình 2.2). Việc vẽ bản đồ đƣợc thực hiện bằng
phần mềm GMT version 2010.
Hình 2.2. Bản đồ dị thƣờng Bughe khu vực Biển Đông và kế cận
Trên cơ sở tờ bản đồ dị thƣờng trọng lực Bughe này, chúng tôi đã thực hiện việc
biến đổi trƣờng trọng lực trong miền tần số đối với khu vực nghiên cứu nằm trong phạm vi bể trầm tích Sông Hồng có kinh độ từ 105.0086 0E đến 110.8038 0E và vĩ độ kéo dài từ 17.0214 0N đến 22.8181 0N ( Hình 2.3). Ở đây việc nâng trƣờng trong
miền tần số đã đƣợc chúng tôi thực hiện theo thuật toán đã trình bày ở chƣơng 1.
16
Hình 2.3. Bản đồ dị thƣờng Bughe khu vực bể trầm tích Sông Hồng và kế cận
2.2. SƠ LƢỢC VỀ ĐẶC ĐIỂM ĐỊA CHẤT - KIẾN TẠO KHU VỰC BỂ TRẦM
TÍCH SÔNG HỒNG
Nằm trên rìa Tây bắc của Biển Đông, bể trầm tích Sông Hồng (hay còn gọi là bể
trầm tích Bắc Bộ ) là bồn trũng Kainozôi chính với diện tích trải rộng khoảng 100.00 km2 bao gồm đồng bằng hạ lƣu sông Hồng, vịnh Bắc Bộ và phần thềm lục
địa Miền trung Việt Nam, bể trầm tích này bao rìa thụ động ở góc đông bắc của
khối lục địa Đông dƣơng.
Những đơn vị kiến tạo chính của bể trầm tích Sông Hồng bao gồm :
- Miền võng Hà nội: Nằm ở phía bắc của bể trầm tích sông Hồng, hầu hết là
phần tách dãn trên lục địa của bể với trầm tích Kainozôi đạt tới chiều dày 10 km.
Trũng nhỏ này là một graben đƣợc phát triển giữa hai cánh chuyển dịch ly tâm (
tách dãn ) của đứt gãy sông Hồng và biến dạng bị khống chế chủ yếu bởi nén ép vào
Đệ tam muộn.
17
- Bồn Bắc Bộ Trung tâm: gồm tâm bồn trũng chính. Chiều dày trầm tích vƣợt
quá 10 km, có lẽ lớn hơn 15 km, bao gồm 3 km trầm tích Đệ tứ gần trục bồn
trũng . Bồn trũng có dạng bất đối xứng với sƣờn phía tây trải rộng, các đứt gãy
cắt thành khối. Biến dạng gần trục bồn đƣợc đặc trƣng bởi điapia phiến sét xảy ra
vào Đệ tứ.
- Graben Quảng Ngãi: úng với phần tách dãn phía nam của bồn Bắc Bộ. Đó là
phụ bồn trũng hẹp với trầm tích Kainozôi dày 6 km. Một móng lồi đã chia graben
với bồn trũng Bắc Bộ Trung tâm về phía bắc. Trục graben đôi nơi bị các đại mạch
grabio tuổi Mĩoxen muộn chia cắt.
- Khối nhô Tri tôn: Đó là một khối nhô dạng hình nêm, nghiêng về phía bắc ở
rìa tây của bồn trũng Bắc Bộ. Dãy núi kéo về phía nam tạo một nền bằng nông, rộng
bị phủ bởi Cacbonate dày. Khối nhô Tri tôn bị phân cách với bồn Bắc Bộ bởi đứt
gãy sông Hồ
- Đứt gãy sông Hồng: Về nguồn gốc có lẽ là một đới khâu Paleozoi hoặc
Mezozoi trong quá trình phát triển Đông dƣơng chủ yếu chuyển dịch trƣợt bằng
trái trong suốt Đệ tam. Nghiên cứu tân kiến tạo ở Trung Quốc cho thấy đới này
chuyển dịch phải trong giai đoạn hiện nay .
2.3. ĐẶC ĐIỂM DỊ THƢỜNG TRỌNG LỰC
Trên phạm vi vịnh Bắc Bộ trƣờng dị thƣờng Bughe (Hình 2.2) có cấu trúc khá
phức tạp. Giá trị trƣờng biến đổi trong phạm vi từ - 35 mgal đến +50 mgal, là nơi
có giá trị trƣờng âm lớn nhất trên thềm lục địa Việt Nam.
Có hai đới dị thƣờng âm lớn nhất thấy khá rõ trên bản đồ dị thƣờng trọng lực
Bughe khu vực vịnh Bắc Bộ, một theo hƣớng bắc -tây bắc và một có hƣớng kinh
tuyến. Đới thứ nhất bắt đầu từ phía ngoài bờ biển Đồ Sơn - Hải Phòng chạy theo
kinh tuyến 107 thẳng đến bờ biển Quảng Bình. Bề rộng của dải này tới 50-60 km
với biên độ dị thƣờng lớn nhất là -35 mgal. Đới thứ hai hẹp hơn có hƣớng bắc -tây
bắc xuất phát từ bờ biển Tam Kỳ - Quảng Ngãi kéo dài về phía trung tâm vịnh Bắc
Bộ, biên độ trung bình là -10 mgal. Phía đông nam đảo Bạch Long Vĩ có một đới dị
thƣờng âm phía đông bắc. Đới này hẹp nhƣng có biên độ lớn đến -30 mgal và chạy
về phía trung tâm nhập vào đới âm thứ nhất đã nêu.
18
Giữa các đới dị thƣờng âm hình thành 3 đới dị thƣờng dƣơng nhỏ hơn với biên
độ 5 - 10 mgal phân bố ở các khu vực: Phía bắc đảo Bạch Long Vĩ có hƣớng vĩ
tuyến, phía ngoài bờ biển Ninh Bình, Thanh Hóa hƣớng kinh tuyến và khá phân dị,
phía ngoài bờ biển vùng bán đảo Sơn Trà theo hƣớng tây bắc.
Phía đông nam của vịnh Bắc Bộ còn có dị thƣờng dƣơng lớn biên độ đạt tới 50
mgal và có hƣớng đông bắc.
2.4. KẾT QUẢ TÍNH TOÁN
Việc biến đổi trƣờng trọng lực khu vực bể trầm tích Sông Hồng nhằm xác định
vị trí các đứt gãy địa chất trong khu vực này đƣợc thực hiện trong miền tần số. Toàn bộ vùng nghiên cứu có kinh độ từ 105.0086 0E đến 110.8038 0E và vĩ độ kéo dài từ 17.0214 0N đến 22.8181 0N đƣợc chia thành mạng lƣới 256 x 256 ô vuông mà nút
mạng là các điểm quan sát, cách đều nhau các khoảng cách 2,5 km. Việc biến đổi
trƣờng bao gồm các bƣớc sau:
a) Biến đổi Furier trƣờng dị thƣờng Bughe thu đƣợc trong vùng nghiên cứu
b) Thực hiện phép nhân chập trƣờng dị thƣờng đã đƣợc biến đổi sang miền tần số
này với các nhân biến đổi khác nhau, trong đó:
với số sóng khi thực hiện việc nâng trƣờng -
với số sóng khi thực hiện việc hạ trƣờng -
và khi thực hiện đạo hàm ngang bậc n -
c) Biến đổi Furier ngƣợc kết quả thu đƣợc để thu đƣợc trƣờng biến đổi trong miền
không gian
Ở đây, kết quả biến đổi trƣờng thu đƣợc đƣợc thực hiện qua việc nâng trƣờng
lên các độ cao z = 0 km, z = 5 km, z = 10 km, z = 15 km, z = 20 km và hạ trƣờng
xuống các độ sâu z = - 0.5 km, z = - 2 km. Ở mỗi mức, các đạo hàm ngang toàn
phần cũng đƣợc chúng tôi tính toán. Kết quả đƣợc thể hiện trên các hình vẽ: hình
2.4; 2.5,….2.14.
19
Hình 2.4. Kết quả hạ trƣờng xuống độ sâu 1km
20
Hình 2.5. Kết quả hạ trƣờng xuống độ sâu 2 km
21
Hình 2.6. Kết quả tính đạo hàm ngang toàn phần ở mức z = 0
22
Hình 2.7. Kết quả nâng trƣờng lên độ cao 5 km
23
Hình 2.8. Kết quả tính đạo hàm ngang toàn phần khi
nâng trƣờng lên độ cao 5 km
24
Hình 2.9. Kết quả nâng trƣờng lên độ cao 10 km
25
Hình 2.10. Kết quả tính đạo hàm ngang toàn phần khi
nâng trƣờng lên độ cao 10 km
26
Hình 2.11. Kết quả nâng trƣờng lên độ cao 15 km
27
Hình 2.12. Kết quả tính đạo hàm ngang toàn phần khi
nâng trƣờng lên độ cao 15 km
28
Hình 2.13. Kết quả nâng trƣờng lên độ cao 20 km
29
Hình 2.14. Kết quả tính đạo hàm ngang toàn phần khi
nâng trƣờng lên độ cao 20 km
30
KẾT LUẬN
Qua việc nghiên cứu áp dụng phƣơng pháp biến đổi trƣờng trong miền tần số
thực hiện việc biến đổi trƣờng trọng lực khu vực bể trầm tích Sông Hồng thuộc
phạm vi thềm lục địa Việt nam, tôi rút ra một số nhận xét và kết luận sau:
Việc hạ trƣờng xuống các độ sâu khác nhau chỉ nhấn mạnh đƣợc các dị
thƣờng địa phƣơng nhƣng không cho phép xác định đƣợc vị trí của các đứt gãy sâu
trong khu vực nghiên cứu.
Do chƣa tách đƣợc các dị thƣờng địa phƣơng trong trƣờng trọng lực quan
sát nên việc tính và xây dựng sơ đồ đạo hàm ngang toàn phần ở mức z = 0 cũng
không cho phép xác định đƣợc vị trí của các đứt gãy sâu trong khu vực nghiên cứu.
Việc tính toán và xây dựng sơ đồ đạo hàm ngang toàn phần của trƣờng trọng
lực đã đƣợc nâng lên ở các mức z = 5 và 10 km phản ánh khá rõ vị trí của các đứt
gãy địa chất trong khu vực nghiên cứu. Kết quả tính toán cho thấy các đứt gãy này
định hƣớng theo hai hƣớng chính là tây bắc – đông nam và đông bắc - tây nam.
Với trƣờng dị thƣờng trọng lực đƣợc nâng lên ở các mức z = 15 và 20 km,
chỉ còn các đứt gãy sâu của khu vực là còn đƣợc phản ánh trong sơ đồ đạo hàm
ngang toàn phần của trƣờng.
Việc sử dụng phần mềm Matlab để thực hiện việc biến đổi trƣờng rất thuận
tiện cho ngƣời sử dụng trong việc tính toán và biểu diễn kết quả. Đặc biệt đối với
bài toán ba chiều.
31
TÀI LIỆU THAM KHẢO
Tiếng Việt
1. Tôn Tích Ái (2003), Trọng lực và thăm dò trọng lực, NXB Đại học quốc gia
Hà Nội.
2. Đỗ Đứ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 ĐHQG Hà Nội.
3. Nguyễn Nhƣ Trung (2014), “Xác định độ sâu móng trầm tích khu vực bể Bắc Bộ
theo phân tích ngƣợc trực tiếp 3D số liệu trọng lực vệ tinh”, Tạp chí các khoa
học về trái đất, 36 (3CD), tr.315-320.
4. V.M. Nikiforov, Phùng Văn Phách, R.G. Kulinich, V.V. Khakhinov,
I.V. Dmitriev, G.N. Skabarnya, M.G. Valitov, Hoàng Văn Vƣợng, Lê Đức Anh,
Nguyễn Văn Điệp, Trần Văn Khá (2014), “Một số kết quả thử nghiệm phƣơng
pháp từ telua nghiên cứu vỏ Trái Đất khu vực Bắc vịnh Bắc Bộ”, Tạp chí các
khoa học về trái đất, 36 (3CD), tr.306-314.
5. Hoàng Văn Vƣợng, Trần Văn Khá, Đào Thị Hà, Nguyễn Kim Dũng (2014),
“Nghiên cứu đặc điểm cấu trúc và mật độ trung bình đất đá trầm tích khu vực
trũng sâu Biển Đông - Quần đảo Trƣờng Sa và kế cận theo tài liệu địa vật lý”,
Tạp chí các khoa học về trái đất, 36 (3CD), tr.321-328.
6. Đoàn Văn Tuyến và nnk (1999), “Đặc điểm cấu trúc sâu đới Sông Hồng trên
khu vực tây bắc vùng trũng Hà Nội theo kết quả phân tích tài liệu từ - telur”,
Tạp chí các khoa học về Trái Đất, 21(1), tr. 31-35
Tiếng Anh
7. Richard J.Brackeyly (1992), Potential theory in gravity and magnetic
application, Cambridge University Press.
PHỤ LỤC
(Listing chương trình máy tính)
close all; clear;clc;
data=importdata('BACBObughe.XYZ');
g=data(:,3);
dd=[' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n',...
' %12.5f %12.5f %12.5f %12.5f %12.5f %12.5f\n'];
fid = fopen('tam.txt', 'wt');
fprintf(fid,dd,g);
fclose(fid);
%edit tam.txt;
fid = fopen('tam.txt', 'r');
df=['%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n', ...
'%g %g %g %g %g %g \n'];
a = fscanf(fid,df, [256 256]);
nx=256;
ny=256;
dlamda=0.0227;
dphi=0.0227;
dx=dphi*110;
dy=dlamda*110;
dz=2.0;
dkx=(2*pi)/(nx*dx);
dky=(2*pi)/(ny*dy);
nyqx=nx/2+1;
nyqy=ny/2+1;
for j=1:nx
if j<=nyqx
kx(j)=(j-1)*dkx;
else
kx(j)=(j-nx-1)*dkx;
end
end
dz=-20.0;
du=(2*pi)/(nx*dx);
dv=(2*pi)/(ny*dy);
nyqx=nx/2+1;
nyqy=ny/2+1;
for j=1:nx
if j<=nyqx
u(j)=(j-1)*du;
else
u(j)=(j-nx-1)*du;
end
end
nyqx=nx/2+1;
nyqy=ny/2+1;
for j=1:nx
if j<=nyqx
u(j)=(j-1)*du;
else
u(j)=(j-nx-1)*du;
end
end
for i=1:ny
if i <= nyqy
v(i)=(i-1)*dv;
else
v(i)=(i-ny-1)*dv;
end
end
for j=1:ny
for i=1:nx
k(i,j)=sqrt(u(j).^2+v(i).^2);
nhan(i,j)=exp(k(i,j)*dz);
for i=1:ny
if i <= nyqy
v(i)=(i-1)*dv;
else
v(i)=(i-ny-1)*dv;
end
end
for j=1:ny
for i=1:nx
k(i,j)=sqrt(u(j).^2+v(i).^2);
nhan(i,j)=exp(k(i,j)*dz);
for i=1:ny
if i <= nyqy
ky(i)=(i-1)*dky;
else
ky(i)=(i-ny-1)*dky;
end
end
for j=1:nx
for i=1:ny
k(i,j)=sqrt(kx(j).^2+ky(i).^2);
nhan(i,j)=exp(k(i,j)*dz);
end
end
h = fft2(a);
for j=1:nx
for i=1:ny
dg(i,j)=h(i,j)*nhan(i,j);
end
end
nangtruong=ifft2(dg);
%dg1=dg1./(nx*ny);
[gzx gzy]=gradient(nangtruong,dx,dy);
daoham=sqrt(gzx.^2+gzy.^2);
x=linspace(17.0214,22.8181,256);
y=linspace(105.0083,110.8038,256);
figure(1);mesh(x,y,nangtruong)
figure(2);mesh(x,y,daoham)
figure(3);contour(x,y,nangtruong,10)
figure(4);contour(x,y,daoham,10)
dg1=nangtruong(:);
data1=[data(:,1) data(:,2) dg1];
data11=data1';
data111=data11(:);
dg2=daoham(:);
data2=[data(:,1) data(:,2) dg2];
data22=data2';
data222=data22(:);
%fprintf(dd,data111);
dd=[' %12.5f %12.5f %12.5f\n'];
fid = fopen('ha2.xyz', 'wt');
% fprintf(dd,h);
for j=1:nx
for i=1:ny
dg(i,j)=h(i,j)*nhan(i,j);
%fprintf(dd,dg);
end
