Giáo trình Xác định độ cao GEOID và dị thường trọng lực từ các hệ số hàm điều hòa cầu - Nguyễn Văn Sáng

Nội dung text: Giáo trình Xác định độ cao GEOID và dị thường trọng lực từ các hệ số hàm điều hòa cầu - Nguyễn Văn Sáng

1. T¹p chÝ KHKT Má - §Þa chÊt, sè 53, 01-2016, tr.58-62 TRẮC ĐỊA – BẢN ĐỒ & QUẢN LÝ ĐẤT ĐAI (trang 58÷67) XÁC ĐỊNH ĐỘ CAO GEOID VÀ DỊ THƯỜNG TRỌNG LỰC TỪ CÁC HỆ SỐ HÀM ĐIỀU HÒA CẦU NGUYỄN VĂN SÁNG, Trường Đại học Mỏ - Địa chất PHẠM VĂN TUYÊN, Công ty cổ phần Dịch vụ và thương mại 568 Tóm tắt: Bài báo trình bày chi tiết các công thức toán học, để tính độ cao geoid và dị thường trọng lực trên cơ sở sử dụng các hệ số hàm điều hòa cầu của các mô hình thế trọng trường và được lập thành chương trình máy tính Geomat2015 bằng ngôn ngữ lập trình Matlab. Các tính toán thực nghiệm được thực hiện với các hệ số hàm điều hòa cầu của mô hình thế trọng trường toàn cầu EGM2008 và vùng thực nghiệm là trên vùng biển Vịnh Bắc Bộ - Việt Nam được biểu diễn ở dạng lưới ô vuông có kích thước (6’ x 6’) với 874 điểm lưới. Kết quả tính toán được so sánh với kết quả được cung cấp bởi tổ chức The International Centre for Global Earth Models (ICGEM) cho thấy sự đúng đắn của kết quả tính toán với các thống kê: độ lệch lớn nhất, nhỏ nhất và độ lệch chuẩn đạt được của độ cao geoid tương ứng là 0,0082m; -0,0030m và ±0,0015m và của dị thường trọng lực tương ứng là 0,0588mgal; -0,2607 mgal và ±0,0264mgal. 1. Đặt vấn đề cầu của các mô hình thế trọng trường khác nhau Trên thế giới, đo cao vệ tinh được ứng dụng nhưng kết quả tính này chỉ ở dưới dạng mắt lưới rất hiệu quả trong nhiều lĩnh vực trong đó có việc ô vuông tùy theo kích thước người sử dụng lựa xác định dị thường trọng lực biển, đã có nhiều chọn. Trong bài toán xác định dị thường trọng quốc gia ứng dụng kết quả đo cao vệ tinh để xác lực biển bằng số liệu đo cao vệ tinh thì số liệu định dị thường trọng lực cho vùng biển của mình. tính toán lại ở dạng các điểm rời rạc. Như vậy, Trong bài toán xác định dị thường trọng lực không thể sử dụng phần mềm sẵn có của tổ chức biển từ số liệu đo cao vệ tinh, có một bước rất trên để tính toán dị thường trọng lực và độ cao quan trọng đó là loại bỏ phần độ cao geoid geoid được. Để khắc phục điều đó trong bài báo (NEGM) trong số liệu đo cao vệ tinh và khôi phục này chúng tôi sẽ giới thiệu chi tiết các công thức lại dị thường trọng lực (∆gEGM) bằng mô hình toán học, để tính ra độ cao geoid và dị thường trường trọng lực toàn cầu theo kỹ thuật “remove trọng lực ở điểm bất kỳ trên cơ sở sử dụng các hệ - restore”. Hiện nay có một số tổ chức trên thế số hàm điều hòa cầu của các mô hình thế trọng giới như ICGE cho phép tính các đại lượng nêu trường và được lập thành chương trình máy tính trên với độ chính xác cao từ hệ số hàm điều hòa để tính toán. 2. Các công thức tính độ cao geoid và dị thường trọng lực từ các hệ số hàm điều hòa cầu Công thức tổng quát xác định độ cao geoid và dị thường trọng lực [1,5,6]: n GM Nmax a n NEGM N0,  Cn,, mcos( m ) S n m sin( m  ) P n m (sin ') , (1) .rr nm 20 n GMa Nmax n gnCm Sm P ( 1)cos(n,, mn )sin( m  ) (sin ') , (2) EGMn m 2  , rr nm 20 trong đó: GM là hằng số trọng trường địa tâm; r là bán kính địa tâm của điểm xét;  là gia tốc lực trọng trường chuẩn trên mặt elipsoid; a là bán trục lớn của ellipsoid; 58
2. ', là vĩ độ và kinh độ địa tâm của điểm xét; Cnm, , S nm, là hệ số điều hòa cầu chuẩn hóa đầy đủ bậc n, hạng m; Pnm, ( sin ' ) là hàm Legendre kết hợp đã chuẩn hóa; N0 là đại lượng mức 0 (zero-degree term). a) Tính bán kính r và vĩ độ địa tâm ’ của từng điểm xét Bán kính địa tâm của điểm xét được tính bằng công thức: ee222(1)sin rxyza()1 222 ’ (3) 1sin e22 acos .cos  a cos .sin  a (1 e2 )sin với: x ;; y z . (4) 1 e2 sin 2 1 e 2 sin 2 1 e 2 sin 2 Độ vĩ địa tâm của điểm xét được tính bằng công thức: 2 zb 'arctanarctantan , (5) 22 xy a trong đó: b là bán trục bé của ellipsoid; , là vĩ độ và kinh độ địa lý của điểm xét; ab22 e2 là tâm sai thứ nhất. a2 b) Tính trọng lực chuẩn trên mặt elipsoid Trọng lực chuẩn trên mặt ellipsoid được tính theo công thức: 1sin k 2 ba  () với k pe, (6) e 22 1sin e a e trong đó: ep, là trọng lực chuẩn trên xích đạo và cực của elipsoid. c) Tính hàm Legendre kết hợp đã chuẩn hóa Đặt tsin' ;uc os' , hàm Legendre kết hợp đã chuẩn hóa có thể được tính theo công thức truy hồi như sau: +) Xét trường hợp n > m: Pn, m(t) an,, m . t . P n 1, m (t) b n m . P n 2, m (t), (7) (2n1)(2n1) (2n 1)(n m 1)(n m 1) với a ; b . (8) nm, (nm)(nm) nm, (n m)(n m)(2n 3) +) Xét trường hợp n = m: - Khi: m 1 ta có: Pm,m(t) u Pmm 1, 1 (t) , 2m m m 21i hoặc Pum,m (t) 3 . (10) i 2 2i d) Tính cos(m) và sin(m) Sin(m ) 2cos sin((m 1) ) sin((m 2) ) . (11) Cos(m ) 2cos  cos((m 1)  ) cos((m 2)  ), (12) +) Với m=0: Sin(m ) 0;cos(m )=1. +) Với m=1: Sin(m ) Sin (  );cos(m  ) cos(  ) . 59
3. e) Tính đại lượng mức 0 (zero-degree term): WGS84 thì đại lượng mức 0 (zero-degree term) N0 Đại lượng này sinh ra khi [2,5]: có giá trị: - Hằng số trọng trường địa tâm toàn cầu (3.986004415 3.986004418).1014 N GM. được sử dụng trong mô hình EGM không 0 6371000x 9.7976432222 bằng hằng số trọng trường địa tâm GM. của 62636855.6693 62636851.7146 0 0.4084m ellipsoid quốc tế. 9.7976432222 - Sự khác nhau của tham số hình học của 3. Xây dựng chương trình tính độ cao geoid và ellipsoid trọng lực khác với ellipsoid quốc tế. dị thường trọng lực từ các hệ số điều hòa cầu GMGMU W Trên cơ sở các công thức từ (1) đến (13), ta có: N 000 , (13) 0 chúng tôi đã tiến hành xây dựng chương trình cho R0. phép tính toán của độ cao geoid và dị thường trong đó: R0: bán kính trung bình của trái đất (R0 trọng lực của điểm bất kỳ khi cho biết các thành = 6371km ); phần tọa độ địa lý ( , ) của các điểm cần tính. W0: thế trọng trường thực của mặt elipsoid toàn cầu; Chương trình có tên là Geomat2015. Sơ đồ U0 : thế trọng trường chuẩn của mặt khối của chương trình được trình bày trên hình 1. elipsoid quốc tế; Chương trình Geomat2015 là một tổ hợp của 10 chương trình con được viết bằng ngôn  : gia tốc lực trọng trường chuẩn trung ngữ lập trình Matlab: EGM_ReadCnmSnm. m; 2 bình (  9.7976432222.ms ). radgra.m; sinmlcosml.m; legfdn.m; hundu.m; Như vậy, áp dụng công thức (13) khi sử dụng Undulation.m; GRA.m; DeltaG.m; N_EGM.m; mô hình EGM2008 (được gắn với elipsoid trọng G_EGM.m. Chương trình được chạy trực tiếp lực TFS2008) để tính toán với ellipsoid quốc tế trên nền của phần mềm Matlab. Geomat201 5 Bắt đầu Đọc file:( ,) Đọc: ; Đọc tham số: GM;a,b,e,p Tính: Tính: r,, ’ Tính : Tính: NGM ; ∆gGM Kết thúc Hình 1. Sơ đồ khối của chương trình Geomat2015 60
4. 4. Tính toán thực nghiệm trong bài báo này, cùng với các tham số khi tính Các tính toán thực nghiệm được thực hiện toán với ellipsoid WGS84: 14 3 -2 với chương trình Geomat2015 và sử dụng các GM = 3.986004418.10 m /s ; bán trục lớn của hệ số hàm điều hòa cầu chuẩn hóa đầy đủ của mô ellipsoid a = 6378137.0 m. hình thế trọng trường toàn cầu EGM2008 với các Các thống kê về độ lệch lớn nhất, nhỏ nhất, hệ số mở rộng tới số bậc 2190 và số hạng 2159 trung bình, độ lệch trung phương và độ lệch do ICGEM cung cấp[3]. Vùng thực nghiệm là chuẩn, được thể hiện trên bảng 1. trên vùng biển Vịnh Bắc Bộ - Việt Nam trong Kết quả tính toán thống kê ở bảng 1 cho thấy, phạm vi vĩ độ (16.97010 ÷ 21.47010) và kinh độ lệch giữa kết quả tính được so với kết quả độ (105.61670 ÷ 108.31670) được biểu diễn ở cung cấp bởi ICGEM là rất nhỏ. Độ lệch trung dạng lưới ô vuông có kích thước (6’ x 6’) với 874 bình nhỏ chứng tỏ trong kết quả tính không có điểm lưới. Các kết quả tính toán được so sánh với chứa sai số hệ thống. Kết quả này thể hiện sự kết quả tính toán do ICGEM cung cấp. chính xác của các công thức và chương trình thiết Lý thuyết và các công thức dùng để tính độ lập được. Độ lệch độ cao geoid và dị thường cao geoid và dị thường trọng lực của ICGEM [1] trọng lực có đồ thị tuân theo luật phân bố chuẩn hoàn toàn tương tự với các công thức đã nêu (hình 2,3). Bảng 1. Thống kê độ lệch độ cao geoid và dị thường trọng lực được tính từ kết quả của chương trình Geomat2015 với kết quả được cung cấp bởi ICGEM Các chỉ tiêu Stt N (m) g (mgal) so sánh E G M E G M 1 Độ lệch lớn nhất Max 0.0082 0.0588 2 Độ lệch nhỏ nhất Min -0.0030 -0.2607 3 Độ lệch trung bình 0.0027 -0.0672 4 Độ lệch trung phương  ±0.0031 ±0.0722 5 Độ lệch tiêu chuẩn  ±0.0015 ±0.0264 Hình 2. Biểu đồ phân bố độ lệch chuẩn giữa kết quả tính toán độ cao geoid với kết quả được cung cấp bởi ICGEM 61
5. Hình 3. Biểu đồ phân bố độ lệch chuẩn giữa kết quả tính toán dị thường trọng lực với kết quả được cung cấp bởi ICGEM 5. Kết luận TÀI LIỆU THAM KHẢO - Các kết quả tính toán của độ cao geoid và dị thường trọng lực bằng chương trình [1]. Franz Barthelmes, 2013. Definition of Geomat2015được so sánh với kết quả cung Functionals of the Geopotential and Their cấp bởi The International Centre for Global Calculation from Spherical Harmonic Models, Earth Models (ICGEM) đã khẳng định sự đúng GFZ German research centre for geosciences. đắn cả về cơ sở lý thuyết và tính toán thực [2]. Hà Minh Hòa, 2014. Lý thuyết và thực tiễn nghiệm của chương trình tính. của trọng lực trắc địa. Nhà xuất bản khoa học và - Chương trình Geomat2015 có thể dùng kỹ thuật, Hà Nội. để khảo sát độ cao geoid và dị thường trọng lực [3]. /ICGEM/ shms/ từ các hệ số điều hòa cầu của nhiều mô hình egm2008.gfc trường trọng lực toàn cầu khác nhau cho các [4]. Martin Vermeer, 2015. Physical Geodesy điểm bất kỳ khi biết các thành phần tọa độ địa lý Maa-6.3271. Toukokuuta. ( , ) của các điểm cần tính. [5]. NIMA,2000. Department of Defense World - Các công thức và chương trình Geodetic System 1984. National Imagery and Geomat2015 có thể được dùng để loại bỏ phần Mapping Agency, America. độ cao geoid trong số liệu đo cao vệ tinh và khôi [6]. Nico Sneeuw, 2006. Physical Geodesy. phục lại dị thường trọng lực bằng mô hình trường Institute of Geodesy, Univesity Stuttgart. trọng lực toàn cầu theo kỹ thuật “remove – [7]. Wolfgang Torge,2001. Geodesy. Third restore” trong bài toán xác định dị thường trọng completely revised and extended edition, Walter lực từ số liệu đo cao vệ tinh. de Gruyter, Berlin, New York. ABSTRACT Determination of geoid height and gravity anomaly from spherical harmonic coefficients Nguyen Van Sang, Hanoi University of Mining and Geology Pham Van Tuyen, JSC service and commercial 568 The article presents the detailed mathematical formulas to calculate the geoid height and gravity ״Geomat2015״ anomaly by using spherical harmonic coefficients and established a computer program by Matlab programming language. The geoid heights and gravity anomalies have been calculated by using spherical harmonic coefficients of the Earth Gravitational Model EGM2008 in Gulf of Tonkin- Vietnam is represented in the form of grid squares of size (6 ' x 6 ') with 874 points grid. The results are compared with the results provided by The International Centre for Global Earth Models (ICGEM) shows the correctness of calculation results with statistics: The maximum, minimum deviation and standard deviation of the geoid heights being 0.0082m; -0.0030m and ± 0.0015m respectively and gravity anomalies being 0.0588mgal; -0.2607 mgal and ±0.0264mgal respectively. 62