Giáo trình Xử lý số liệu trắc địa: Phần 2 - PGS.TS Đặng Nam Chinh (Chủ biên)
lượt xem 6
download
Nối tiếp phần 1, phần 2 của giáo trình "Xử lý số liệu trắc địa" do PGS.TS Đặng Nam Chinh chủ biên tiếp tục trình bày những nội dung về: phân tích thống kê và xấp xỉ hàm; xác định các đặc trưng thống kê của các dãy số liệu quan trắc; khái niệm về phân tích hồi quy; nội suy và các ứng dụng; nội suy theo khoảng cách; nội suy theo hàm đa thức; nội suy Lagrange; nội suy Collocation; nội suy Kriging;.. Mời các bạn cùng tham khảo!
Bình luận(0) Đăng nhập để gửi bình luận!
Nội dung Text: Giáo trình Xử lý số liệu trắc địa: Phần 2 - PGS.TS Đặng Nam Chinh (Chủ biên)
- Chương 2 PHÂN TÍCH THỐNG KÊ VÀ XẤP XỈ HÀM 2.1. DÃY SỐ LIỆU QUAN TRẮC Trong thực tế, người ta thường tiến hành đo đạc hoặc quan trắc một đại lượng hay yếu tố nào đó nhiều lần (đo góc, đo cạnh, ...) hoặc một loại đại lượng hay yếu tố nào đó được phân bố trong không gian (đo trọng lực, đo cường độ từ trường, độ cao bề mặt nào đó) hoặc quan trắc đại lượng có tính biến đổi theo thời gian (trị đo pha, đo khoảng cách giả trong công nghệ GPS, giá trị quan trắc lún, giá trị quan trắc dịch chuyển, biến dạng, ...). Các số liệu quan trắc đó được gọi chung là dãy số liệu hay bộ số liệu (Data Samples). Một đặc điểm chung của đo đạc hay quan trắc là được tiến hành trong cùng điều kiện hoặc không cùng điều kiện bằng thiết bị quan trắc để nhận được các giá trị quan trắc luôn kèm theo sai số. Nếu xét trong không gian, các trị đo chỉ có thể được quan trắc tại những vị trí (điểm) nhất định với số lượng quan trắc là hữu hạn mà không thể quan trắc tất cả vị trí trong không gian đó. Nếu xét theo thời gian, các trị quan trắc cũng chỉ có thể thực hiện tại những thời điểm nhất định với tần suất cao hoặc thấp chứ không thể quan trắc trên toàn bộ trục thời gian. Như vậy đặc điểm chung của dãy số liệu quan trắc là tập hợp các số liệu rời rạc. Ký hiệu u là véc tơ tọa độ trong hệ thống tọa độ, (thường là x, y trên một mặt nào đó) và ký hiệu z là giá trị đo (quan trắc) xác định tại vị trí tọa độ u, khi đó có thể coi z như là “hàm” của tọa độ u và được viết là z(u). Xử lý dãy số liệu rời rạc không chỉ là nhiệm vụ của ngành trắc địa mà là của nhiều lĩnh vực khác như môi trường, địa vật lý, địa chất công trình, địa chất thủy văn,... Để xử lý dãy số liệu quan trắc trước hết phải xác định được một số đặc trưng thống kê của chúng. Đối với dãy số liệu phân bố trong không gian, cần phải xác định một số đặc trưng địa thống thống kê (Geostatistics) như: tính chất của dãy số liệu (đẳng hướng hay không đẳng hướng, phân bố đều hay không đều), đặc trưng phân bố của dữ liệu (chuẩn hay không chuẩn), trị trung bình, độ lệch chuẩn, mối liên hệ không gian của các cặp giá trị, ... Liên quan đến các tính chất này, phải xác định các giá trị như phương sai (Variance), hiệp phương sai (Covariance), hệ số tương quan (Correlation), bán phương sai (Semivariance), bước gián đoạn (trễ) (Lag). Trong đó người ta sử dụng các thuật ngữ như “Variogram”, “Semivariogram” hoặc “Correlogram” để chỉ đặc trưng tương 73
- quan không gian của dãy số liệu. Các tham số đặc trưng cho mối liên hệ không gian của dãy số liệu được người ta sử dụng để ước lượng địa thống kê (Geostatistical Estimation) như để nội suy, ngoại suy giá trị tại những vị trí không quan trắc và sử dụng để mô phỏng địa thống kê (Geostatistical Simulation) thể hiện quy luật biến đổi trong không gian của số liệu, lập mô hình,... Đối với dãy số liệu quan trắc theo chuỗi thời gian, cần xác định đặc trưng thống kê của chúng như tính chất dừng hay không dừng (Stationary or Nonstationary), xu thế biến thiên của dãy số liệu, ... Nếu như chuỗi thời gian đó là dừng, cũng sẽ xác định được một số đặc trưng thống kê tương tự như trị trung bình, phương sai, độ trễ, hiệp phương sai hay tự phương sai, ... Nếu có dãy số liệu Z(u) và dãy số liệu G(u), khi đó sẽ có các khái niệm về hiệp phương sai chéo (Cross-covariance), hệ số tương quan chéo (Cross- correlation) cũng như variogram chéo (CrossVariogram), ... 2.2. XÁC ĐỊNH CÁC ĐẶC TRƯNG THỐNG KÊ CỦA CÁC DÃY SỐ LIỆU QUAN TRẮC Như đã nói ở trên, xác định các đặc trưng thống kê của dãy số liệu quan trắc phân bố trong không gian là cần thiết cho công tác xử lý dãy số liệu quan trắc. Sau đây giới thiệu một số công thức tính toán phân tích, trong đó có một số công thức chỉ sử dụng mà không chứng minh. 2.2.1. Tính chất phân bố của dãy số liệu Hình 2.1. Phân bố các điểm quan trắc độ rỗng 74
- Với tập hợp số liệu quan trắc (số lượng khá lớn) có thể sử dụng đồ thị để biểu thị mật độ phân bố các giá trị của tập hợp đó bằng cách chia các khoảng giá trị và xác định số lượng giá trị trong các khoảng đó. Ví dụ dựa trên 85 điểm đo xác định giá trị độ rỗng trung bình (Averaged Porosity Values) có giá trị trong khoảng từ 12% đến 17% phân bố trên một diện tích 20 kmx16 km (hình 2.1), có thể biểu thị đặc tính phân bố chuẩn của bộ số liệu đo trên hình 2.2 [21, 22]. Hình 2.2. Quy luật phân bố của các giá trị quan trắc Nhận thấy rằng, giá trị độ rỗng của 85 điểm quan trắc có quy luật phân bố rất gần với luật phân bố chuẩn. Trong trắc địa vật lý (Physical Geodesy), các giá trị trọng lực g được đo ở các điểm khác nhau với một mật độ nào đó cũng hình thành nên một dãy số liệu phân bố trong không gian. 2.2.2. Trị trung bình và phương sai Từ các giá trị quan trắc, dễ dàng tính được trị trung bình theo công thức: 1 n Z(u ) tb z( u i ) (2.2.1.) n i 1 trong đó n là số trị quan trắc. Phương sai (Variance) của dãy số liệu được tính theo công thức: 1 n z(u i ) Z(u ) tb 2 Var (z(u )) (2.2.2) n i 1 Thứ nguyên của phương sai bằng bình phương của thứ nguyên z(u). 75
- 2.2.3. Hiệp phương sai và hệ số tương quan Để tiếp cận với khái niệm hiệp phương sai, cần xét tới biến mới là độ trễ (lag) hay khoảng cách (xét tương quan) giữa các cặp điểm, ký hiệu là s. Ứng với giá trị s khác nhau sẽ tính được giá trị hiệp phương sai khác nhau. Lúc này hiệp phương sai được coi là hàm của khoảng cách s, ký hiệu là C(s). Khái niệm hiệp phương sai nêu trên được xét trong một tập hợp số liệu cho nên còn được gọi là tự hiệp phương sai (Auto-Covariance). Nếu như dãy số liệu phân bố trong không gian thỏa mãn tính đẳng hướng (Non-trend) và có tính chất dừng (Stationary), thì hiệp phương sai (Covariance) giữa các giá trị z(u) và z(u+s) được tính theo công thức: 1 N (s ) C(s) [z(u i ) Z(u ) tb ][z(u i s) Z(u s) tb ] (2.2.3) N(s) i 1 1 n với Z(u s) tb z ( u i s) (2.2.4) n i1 Trong (2.2.3), N(s) là số cặp điểm có khoảng cách là s. Số cặp điểm N(s) càng lớn thì giá trị hiệp phương sai tính được C(s) càng có độ tin cậy cao. Từ (2.2.3) có thể chứng minh được rằng: 1 N (s ) C(s) z(u i )z(u i s) Z(u ) tb Z(u s) tb (2.2.5) N (s) i 1 Vì cùng xuất phát từ một dãy giá trị z(u i ) , do đó trị trung bình tính theo (2.2.1) và trị trung bình tính theo (2.2.4) sẽ như nhau, khi đó biểu thức (2.2.5) sẽ được viết dưới dạng: 1 N (s ) 2 C(s) z( u i ) z(u i s) [ Z(u ) tb ] (2.2.6) N(s) i 1 Vì vậy phương sai Var (z(u )) s tính theo z(u i s) cũng bằng giá trị phương sai Var(z(u )) tính theo z(u i ) ở công thức (2.2.2). Theo định nghĩa, hệ số tương quan ứng với khoảng cách s được tính: C(s) C(s) (s) (2.2.7) Var (z( u )).Var (z(u )) s Var (z(u )) Hiệp phương sai tính theo (2.2.3) hoặc (2.2.6) là các giá trị hiệp phương sai thực nghiệm. Quy luật chung của hiệp phương sai thực nghiệm thể hiện trên hình 2.3. Dựa vào quy luật của các giá trị hiệp phương sai thực nghiệm, có thể xác định được dạng của hàm hiệp phương sai lý thuyết tương ứng theo nguyên 76
- tắc xấp xỉ hàm. Hàm hiệp phương sai lý thuyết có vai trò quan trọng trong nội suy theo phương pháp Collocation và nội suy theo phương pháp Kriging. Có thể nhận thấy rằng, hiệp phương sai tính theo (2.2.3) khi s = 0 chính là phương sai tính theo công thức (2.2.2). 2.2.4. Bán phương sai Tương tự như hiệp phương sai, bán phương sai (Semivariance) là giá trị đặc trưng cho mối phụ thuộc của các giá trị z(u) theo khoảng cách s, tuy nhiên nó cũng có điểm khác với hiệp phương sai. Hiệp phương sai và bán phương sai đều là các chỉ tiêu thống kê về mức độ (đo) sự tương quan không gian (Spatial Autocorrelation). Bán phương sai còn được sử dụng cho thuật ngữ Semi- variogram hoặc Variogram. Bán phương sai ở khoảng cách s được tính theo công thức kỳ vọng: 1 (s ) E[(z(u ) z(u s)) 2 ] (2.2.8) 2 hoặc viết ở dạng công thức thực dụng: 2 1 N (s ) (s) z(u i ) z(u i s) (2.2.9) 2 N(s) i 1 Bán phương sai tính theo công thức trên gọi là bán phương sai thực nghiệm. Quy luật chung của bán phương sai là khi s=0 giá trị (s) 0 và tăng tới một giới hạn nào đó (Sill), giới hạn tăng sẽ có khoảng cách tương ứng gọi là độ trễ tương quan. Đồ thị chung của bán phương sai thực nghiệm, phương sai thực nghiệm và hệ số tương quan có dạng như hình 2.3. Hình 2.3. Đồ thị phương sai, bán phương sai và hệ số tương quan 77
- Dựa vào quy luật của các giá trị bán phương sai để chọn hàm bán phương sai lý thuyết phù hợp và dựa trên nguyên tắc xấp xỉ hàm để xác định các tham số của hàm bán phương sai lý thuyết. Hàm bán phương sai lý thuyết có vài trò quan trọng trong nội suy Kriging. Mối liên hệ giữa bán phương sai và hiệp phương sai được thể hiện qua công thức sau: (s) C(0) C(s) (2.2.10) trong đó C(0) là hiệp phương sai tương ứng với s = 0 và cũng chính là phương sai. 2.3. KHÁI NIỆM VỀ PHÂN TÍCH HỒI QUY 2.3.1. Khái niệm chung Đối với dãy số liệu quan trắc theo chuỗi thời gian, có thể gặp những trường hợp là chuỗi ngẫu nhiên không dừng (Non-Stationary). Trong trường hợp này, cần phải xác định quy luật biến đổi của chuỗi số liệu hay xác định ma trận trạng thái của quy luật đó. Hồi quy có liên quan với phân tích hệ số tương quan, dựa trên hệ số tương quan người ta đưa ra (lựa chọn) phương pháp hồi quy phù hợp. Đối với các số liệu đo thay đổi liên tục theo thời gian như các trị đo pha trị đo khoảng cách giả trong công nghệ GPS, chúng ta cần xem xét quy luật biến thiên theo thời gian của các trị đo đó đồng thời có thể xem xét cả quy luật biến thiên của vận tốc trị đo và gia tốc trị đo. Ký hiệu các trị đo pha sóng tải tại thời điểm t i1 và tại thời điểm t i là L( t i1 ) và L( t i ) , vận tốc trị đo pha tại thời điểm t i ký hiệu là VL ( t i ) được tính theo công thức: dL( t i ) L( t i ) L( t i 1 ) L( t i ) VL ( t i ) (2.2.11) dt t i t i 1 t trong đó t là tần suất ghi tín hiệu GPS. Ký hiệu a L ( t i ) là gia tốc trị đo pha tại thời điểm t i , được tính như sau:: d 2 L( t i ) VL ( t i ) VL ( t i 1 ) a L (t i ) (2.2.12) dt 2 t Nhờ khảo sát sự biến thiên của vận tốc trị đo pha và gia tốc trị đo pha chúng ta có thể phát hiện được hiện tượng trượt chu kỳ ... 78
- Ví dụ: Xét một chuỗi số liệu trị đo pha sóng tải L1 của vệ tinh PRN-19 bằng máy thu 1 tần số Trimble 4600LS, tần suất ghi tín hiệu 15s, bắt đầu từ 7h 23m 30s đến 8h 20m 45s ngày 15 tháng 10 năm 2006 tại Hà Nội (gồm 230 trị đo pha), tính được các giá trị thay đổi trị đo pha như vận tốc trị đo pha (dL/dt) và gia tốc trị đo pha (d2L/dt2). Kết quả trị đo pha và các giá trị tính toán được thể hiện trên bảng sau: TT Thời gian L1 dL/dt d2L/dt2 1 7 23 30.0 223519.116 5206.061 -.317 2 7 23 45.0 301610.037 5201.300 -.318 3 7 24 .0 379629.540 5196.531 -.315 4 7 24 15.0 457577.500 5191.803 -.320 5 7 24 30.0 535454.543 5187.002 -.300 6 7 24 45.0 613259.580 5182.500 -.288 7 7 25 .0 690997.078 5178.175 -.270 8 7 25 15.0 768669.705 5174.123 -.248 9 7 25 30.0 846281.548 5170.401 -.267 10 7 25 45.0 923837.570 5166.402 -.246 11 7 26 .0 1001333.595 5162.713 -.225 12 7 26 15.0 1078774.294 5159.333 -.241 13 7 26 30.0 1156164.288 5155.715 -.241 14 7 26 45.0 1233500.009 5152.103 -.215 15 7 27 .0 1310781.555 5148.878 -.205 16 7 27 15.0 1388014.722 5145.805 -.192 17 7 27 30.0 1465201.803 5142.926 -.195 18 7 27 45.0 1542345.698 5140.007 -.221 19 7 28 .0 1619445.800 5136.698 -.174 20 7 28 15.0 1696496.263 5134.087 -.223 21 7 28 30.0 1773507.571 5130.739 -.243 22 7 28 45.0 1850468.657 5127.088 -.190 23 7 29 .0 1927374.980 5124.238 -.159 24 7 29 15.0 2004238.553 5121.849 -.169 25 7 29 30.0 2081066.285 5119.314 -.180 26 7 29 45.0 2157855.988 5116.614 -.168 27 7 30 .0 2234605.199 5114.095 -.196 28 7 30 15.0 2311316.618 5111.153 -.178 29 7 30 30.0 2387983.910 5108.485 -.174 30 7 30 45.0 2464611.188 5105.875 -.178 31 7 31 .0 2541199.309 5103.201 -.172 32 7 31 15.0 2617747.324 5100.614 -.174 33 7 31 30.0 2694256.527 5098.006 -.170 34 7 31 45.0 2770726.614 5095.457 -.204 35 7 32 .0 2847158.463 5092.401 -.177 36 7 32 15.0 2923544.471 5089.746 -.167 37 7 32 30.0 2999890.656 5087.245 -.172 38 7 32 45.0 3076199.329 5084.665 -.167 39 7 33 .0 3152469.309 5082.164 -.179 40 7 33 15.0 3228701.767 5079.483 -.186 41 7 33 30.0 3304894.016 5076.696 -.185 42 7 33 45.0 3381044.453 5073.914 -.149 43 7 34 .0 3457153.168 5071.681 -.154 44 7 34 15.0 3533228.386 5069.377 -.190 45 7 34 30.0 3609269.038 5066.520 -.155 46 7 34 45.0 3685266.835 5064.190 -.178 79
- 47 7 35 .0 3761229.680 5061.514 -.135 48 7 35 15.0 3837152.384 5059.489 -.164 49 7 35 30.0 3913044.720 5057.027 -.165 50 7 35 45.0 3988900.131 5054.553 -.150 51 7 36 .0 4064718.432 5052.301 -.173 52 7 36 15.0 4140502.953 5049.701 -.160 53 7 36 30.0 4216248.461 5047.304 -.177 54 7 36 45.0 4291958.017 5044.643 -.161 55 7 37 .0 4367627.658 5042.232 -.156 56 7 37 15.0 4443261.143 5039.885 -.163 57 7 37 30.0 4518859.415 5037.440 -.178 58 7 37 45.0 4594421.017 5034.763 -.169 59 7 38 .0 4669942.460 5032.226 -.158 60 7 38 15.0 4745425.849 5029.851 -.167 61 7 38 30.0 4820873.616 5027.353 -.162 62 7 38 45.0 4896283.904 5024.920 -.170 63 7 39 .0 4971657.704 5022.370 -.145 64 7 39 15.0 5046993.251 5020.194 -.146 65 7 39 30.0 5122296.159 5017.999 -.190 66 7 39 45.0 5197566.143 5015.142 -.160 67 7 40 .0 5272793.269 5012.736 -.160 68 7 40 15.0 5347984.304 5010.331 -.151 69 7 40 30.0 5423139.268 5008.065 -.154 70 7 40 45.0 5498260.248 5005.759 -.151 71 7 41 .0 5573346.633 5003.492 -.148 72 7 41 15.0 5648399.008 5001.275 -.151 73 7 41 30.0 5723418.127 4999.017 -.169 74 7 41 45.0 5798403.386 4996.487 -.158 75 7 42 .0 5873350.692 4994.116 -.164 76 7 42 15.0 5948262.433 4991.651 -.164 77 7 42 30.0 6023137.194 4989.193 -.144 78 7 42 45.0 6097975.083 4987.035 -.128 79 7 43 .0 6172780.610 4985.109 -.132 80 7 43 15.0 6247557.248 4983.129 -.191 81 7 43 30.0 6322304.186 4980.264 -.163 82 7 43 45.0 6397008.152 4977.813 -.161 83 7 44 .0 6471675.344 4975.398 -.155 84 7 44 15.0 6546306.310 4973.074 -.146 85 7 44 30.0 6620902.414 4970.882 -.164 86 7 44 45.0 6695465.643 4968.415 -.157 87 7 45 .0 6769991.869 4966.058 -.145 88 7 45 15.0 6844482.743 4963.888 -.138 89 7 45 30.0 6918941.059 4961.816 -.137 90 7 45 45.0 6993368.293 4959.767 -.179 91 7 46 .0 7067764.800 4957.080 -.160 92 7 46 15.0 7142121.004 4954.683 -.130 93 7 46 30.0 7216441.250 4952.739 -.143 94 7 46 45.0 7290732.340 4950.599 -.151 95 7 47 .0 7364991.322 4948.337 -.142 96 7 47 15.0 7439216.383 4946.205 -.165 97 7 47 30.0 7513409.456 4943.725 -.156 98 7 47 45.0 7587565.335 4941.389 -.138 99 7 48 .0 7661686.169 4939.319 -.142 100 7 48 15.0 7735775.949 4937.184 -.150 101 7 48 30.0 7809833.713 4934.939 -.132 102 7 48 45.0 7883857.791 4932.956 -.169 103 7 49 .0 7957852.130 4930.420 -.132 80
- 104 7 49 15.0 8031808.424 4928.434 -.148 105 7 49 30.0 8105734.930 4926.209 -.149 106 7 49 45.0 8179628.066 4923.969 -.162 107 7 50 .0 8253487.596 4921.537 -.152 108 7 50 15.0 8327310.657 4919.261 -.122 109 7 50 30.0 8401099.577 4917.431 -.156 110 7 50 45.0 8474861.049 4915.096 -.168 111 7 51 .0 8548587.489 4912.572 -.155 112 7 51 15.0 8622276.067 4910.253 -.102 113 7 51 30.0 8695929.856 4908.717 -.184 114 7 51 45.0 8769560.614 4905.950 -.151 115 7 52 .0 8843149.869 4903.690 -.135 116 7 52 15.0 8916705.223 4901.663 -.152 117 7 52 30.0 8990230.166 4899.379 -.138 118 7 52 45.0 9063720.851 4897.304 -.145 119 7 53 .0 9137180.412 4895.125 -.169 120 7 53 15.0 9210607.280 4892.585 -.147 121 7 53 30.0 9283996.057 4890.387 -.164 122 7 53 45.0 9357351.858 4887.932 -.162 123 7 54 .0 9430670.838 4885.504 -.165 124 7 54 15.0 9503953.394 4883.033 -.134 125 7 54 30.0 9577198.890 4881.020 -.150 126 7 54 45.0 9650414.195 4878.773 -.140 127 7 55 .0 9723595.787 4876.673 -.148 128 7 55 15.0 9796745.884 4874.447 -.154 129 7 55 30.0 9869862.594 4872.143 -.168 130 7 55 45.0 9942944.738 4869.622 -.148 131 7 56 .0 10015989.062 4867.402 -.125 132 7 56 15.0 10089000.095 4865.528 -.162 133 7 56 30.0 10161983.018 4863.104 -.147 134 7 56 45.0 10234929.571 4860.902 -.147 135 7 57 .0 10307843.107 4858.698 -.148 136 7 57 15.0 10380723.578 4856.473 -.150 137 7 57 30.0 10453570.670 4854.228 -.149 138 7 57 45.0 10526384.083 4851.987 -.164 139 7 58 .0 10599163.887 4849.531 -.137 140 7 58 15.0 10671906.856 4847.478 -.144 141 7 58 30.0 10744619.022 4845.322 -.138 142 7 58 45.0 10817298.849 4843.258 -.141 143 7 59 .0 10889947.723 4841.141 -.159 144 7 59 15.0 10962564.844 4838.760 -.142 145 7 59 30.0 11035146.249 4836.630 -.121 146 7 59 45.0 11107695.705 4834.808 -.126 147 8 0 .0 11180217.824 4832.917 -.149 148 8 0 15.0 11252711.581 4830.678 -.180 149 8 0 30.0 11325171.744 4827.976 -.123 150 8 0 45.0 11397591.385 4826.135 -.128 151 8 1 .0 11469983.414 4824.210 -.148 152 8 1 15.0 11542346.568 4821.988 -.116 153 8 1 30.0 11614676.395 4820.254 -.154 154 8 1 45.0 11686980.199 4817.939 -.143 155 8 2 .0 11759249.282 4815.797 -.145 156 8 2 15.0 11831486.243 4813.620 -.129 157 8 2 30.0 11903690.544 4811.682 -.148 158 8 2 45.0 11975865.768 4809.461 -.156 159 8 3 .0 12048007.676 4807.117 -.116 160 8 3 15.0 12120114.424 4805.377 -.134 81
- 161 8 3 30.0 12192195.082 4803.366 -.127 162 8 3 45.0 12264245.579 4801.466 -.125 163 8 4 .0 12336267.570 4799.588 -.163 164 8 4 15.0 12408261.387 4797.144 -.142 165 8 4 30.0 12480218.545 4795.013 -.124 166 8 4 45.0 12552143.745 4793.147 -.125 167 8 5 .0 12624040.950 4791.276 -.151 168 8 5 15.0 12695910.089 4789.017 -.139 169 8 5 30.0 12767745.344 4786.926 -.142 170 8 5 45.0 12839549.231 4784.799 -.138 171 8 6 .0 12911321.211 4782.734 -.138 172 8 6 15.0 12983062.222 4780.659 -.129 173 8 6 30.0 13054772.110 4778.720 -.150 174 8 6 45.0 13126452.914 4776.465 -.158 175 8 7 .0 13198099.882 4774.098 -.130 176 8 7 15.0 13269711.356 4772.142 -.127 177 8 7 30.0 13341293.482 4770.236 -.110 178 8 7 45.0 13412847.025 4768.583 -.139 179 8 8 .0 13484375.768 4766.504 -.131 180 8 8 15.0 13555873.335 4764.536 -.153 181 8 8 30.0 13627341.381 4762.242 -.152 182 8 8 45.0 13698775.006 4759.957 -.141 183 8 9 .0 13770174.354 4757.849 -.108 184 8 9 15.0 13841542.094 4756.225 -.131 185 8 9 30.0 13912885.473 4754.254 -.144 186 8 9 45.0 13984199.290 4752.099 -.165 187 8 10 .0 14055480.774 4749.619 -.134 188 8 10 15.0 14126725.063 4747.607 -.099 189 8 10 30.0 14197939.172 4746.122 -.140 190 8 10 45.0 14269131.009 4744.015 -.154 191 8 11 .0 14340291.233 4741.710 -.114 192 8 11 15.0 14411416.886 4740.001 -.119 193 8 11 30.0 14482516.908 4738.210 -.151 194 8 11 45.0 14553590.063 4735.946 -.139 195 8 12 .0 14624629.254 4733.862 -.155 196 8 12 15.0 14695637.183 4731.542 -.136 197 8 12 30.0 14766610.310 4729.509 -.133 198 8 12 45.0 14837552.947 4727.507 -.145 199 8 13 .0 14908465.551 4725.331 -.140 200 8 13 15.0 14979345.511 4723.231 -.137 201 8 13 30.0 15050193.975 4721.178 -.144 202 8 13 45.0 15121011.638 4719.012 -.143 203 8 14 .0 15191796.820 4716.870 -.145 204 8 14 15.0 15262549.865 4714.691 -.148 205 8 14 30.0 15333270.234 4712.478 -.154 206 8 14 45.0 15403957.399 4710.169 -.131 207 8 15 .0 15474609.932 4708.198 -.145 208 8 15 15.0 15545232.896 4706.027 -.135 209 8 15 30.0 15615823.295 4703.996 -.149 210 8 15 45.0 15686383.233 4701.767 -.161 211 8 16 .0 15756909.732 4699.351 -.157 212 8 16 15.0 15827400.000 4696.998 -.145 213 8 16 30.0 15897854.971 4694.826 -.119 214 8 16 45.0 15968277.361 4693.047 -.143 215 8 17 .0 16038673.065 4690.907 -.145 216 8 17 15.0 16109036.670 4688.738 -.140 217 8 17 30.0 16179367.739 4686.636 -.152 82
- 218 8 17 45.0 16249667.276 4684.359 -.162 219 8 18 .0 16319932.662 4681.929 -.163 220 8 18 15.0 16390161.596 4679.483 -.133 221 8 18 30.0 16460353.848 4677.494 -.154 222 8 18 45.0 16530516.263 4675.190 -.149 223 8 19 .0 16600644.117 4672.962 -.167 224 8 19 15.0 16670738.553 4670.461 -.132 225 8 19 30.0 16740795.472 4668.485 -.138 226 8 19 45.0 16810822.746 4666.413 -.148 227 8 20 .0 16880818.940 4664.194 -.148 228 8 20 15.0 16950781.846 4661.974 -.147 229 8 20 30.0 17020711.452 4659.771 230 8 20 45.0 17090608.013 Từ các số liệu ở bảng trên, có thể nhận thấy một số đặc tính thống kê của dãy số liệu trên qua các đồ thị hình 2.4, hình 2.5 và hình 2.6. Hình 2.4. Đồ thị biến thiên trị đo pha L1 Hình 2.5. Đồ thị biến thiên vận tốc trị đo pha (dL/dt) Hình 2.6. Đồ thị biến thiên gia tốc trị đo pha (d2L/dt2) 83
- Nhận thấy trong đoạn số liệu trên, trị đo pha tăng dần (hình 2.4) còn vận tốc trị đo pha lại giảm dần (hình 2.5), các giá trị gia tốc trị đo pha có thể hiện yếu tố ngẫu nhiên chứa trong các trị đo pha (hình 2.6). Cũng từ hình 2.6 có thể nhận thấy rằng khi mới bật máy thu, gia tốc pha chưa ổn định, nhưng sau đó giá trị gia tốc pha khá ổn định mặc dù có sự biến thiên, thể hiện tính chất của một biến ngẫu nhiên dừng. Trong phần này, chúng ta sẽ xem xét 2 vấn đề thường được áp dụng đối với các chuỗi số liệu quan trắc là hồi quy (Regression) và tự hồi quy (Auto- regression). 2.3.2. Hệ số tương quan thực nghiệm Xét hai biến ngẫu nhiên X và Y, được cho bởi dãy các giá trị thực nghiệm rời rạc như sau: Bảng 2.1. Các giá trị của hai biến X và Y X x1 x2 x3 ... xn Y y1 y2 y3 ... yn Phương sai của X và Y ký hiệu là V(X) và V(Y) được tính: 2 2 1 n 1 n Var (X) i ( x x TB ) ; Var ( Y ) i ( y y TB ) (2.3.1) (n 1) i1 (n 1) i1 trong đó: n n xi yi i 1 i 1 x TB ; y TB (2.3.2) n n Hiệp phương sai chéo giữa X và Y được tính: 1 n COV (X, Y) ( x i x TB )( y i y TB ) (2.3.3) (n 1) i 1 Hệ số tương quan (mẫu, hay thực nghiệm) giữa biến X và Y được xác định bởi công thức sau: COV(X, Y) xy (2.3.4) Var (X).Var (Y) Thay (2.3.1) và (2.3.3) vào (2.3.4) ta được: 84
- n ( x i x TB )( y i y TB ) i 1 xy (2.3.5) n n 2 2 ( x i x TB ) ( y i y TB ) i 1 i 1 Hệ số tương quan xy xác định theo (2.3.5) có giá trị nằm trong khoảng - 1 đến +1 tức là: 1 1 (2.3.6) Nếu xy 1 , khi đó X và Y có tương quan tuyến tính dương tuyệt đối ( đồng biến) Nếu xy 1 , khi đó X và Y có tương quan tuyến tính âm tuyệt đối ( nghịch biến) Nếu xy 0 , khi đó X và Y không tương quan tuyến tính. Từ các giá trị xi, yi có thể triển vẽ lên đồ thị, từ đó có thể nhận dạng được quy luật của mối quan hệ X, Y . Đồ thị trên hình 2.7 thể hiện mối quan hệ giữa X và Y là tuyến tính. Hình 2.7. Quan hệ tuyến tính 2.3.3. Hồi quy tuyến tính Nếu xác định được hệ số tương quan thực nghiệm có giá trị tuyệt đối bằng 1 hoặc xấp xỉ 1, khi đó có thể biểu diễn X, Y bởi hàm hồi quy tuyến tính như sau: Y = A + B.X (2.3.7) Trên thực tế chỉ có thể ước lượng được các tham số của hàm là a và b theo tiêu chuẩn độ lệch nhỏ nhất (ước lượng không chệch), khi đó hàm hồi quy sẽ là: 85
- y i a b.x i i (2.3.8) trong đó i là sai số của hàm hồi quy. Dựa vào các dãy giá trị xi, yi sẽ xác định được giá trị các tham số a, b của hàm hồi quy tuyến tính theo nguyên lý bình phương nhỏ nhất. Phương pháp xác định a, b thực hiện theo nguyên lý bình phương nhỏ nhất, sẽ được trình bày ở phần “Xấp xỉ hàm”. Ngoài hồi quy tuyến tính, có thể sử dụng hàm hồi quy phi tuyến (bậc hai hoặc bậc 3,...) nếu như giá trị tuyệt đối của hệ số tương quan thực nghiệm không gần với 1. 2.3.4. Phân tích tự hồi quy Phân tích tự hồi quy AR (Auto-Regression Analysis) là phương pháp xử lý số liệu chỉ dựa trên các giá trị của một chuỗi số liệu thời gian, nhằm đưa ra dự báo được coi là tốt nhất tại thời điểm t dựa trên các giá trị trước đó. Với chuỗi số liệu thời gian Y(Y1 , Y2 ,..., Yn ) , mô hình tự hồi quy có dạng tổng quát như sau: Yt f (Yt 1 , Yt 2 ,...., Yt p , ) (2.3.9) trong đó là nhiễu trắng (ngẫu nhiên) Mô hình tự hồi quy thường sử dụng có dạng hàm tuyến tính sau: p Yt b0 bi Yt i t (2.3.10) i 1 trong đó: Yt giá trị biến Y ở thời điểm t; Yt i (i = 1, 2, ..., p) là giá trị biến Y ở thời điểm t-i; bi là hệ số hồi quy ; p là bậc tự hồi quy; t là nhiễu. Trong trường hợp này cần xác định các hệ số bi (i 0,1,2,... p) thỏa mãn điều kiện bình phương tối thiểu, tức là bài toán ước lượng bình phương nhỏ nhất. Từ biểu thức (2.3.10) viết được hệ phương trình số hiệu chỉnh ở dạng ma trận như sau: V AX L (2.3.11) trong đó: 86
- 1 Y p Y p 1 ... Y1 Y p 1 V1 1 Y b o Y V p 1 Y p ... Y2 b 1 p2 V 2 ; A . . . . . ; X ; L . (2.3.12) .. .. . . . . . . Vn p 1 Yn1 Yn 2 ... Yn p b p Yn Theo điều kiện V T V min , véc tơ X được xác định: X ( A T A ) 1 A T L (2.3.13) Có hai dạng tự hồi quy thường gặp đó là tự hồi quy bậc nhất (p = 1) và tự hồi quy bậc 2, (p = 2) Mô hình tự hồi quy bậc nhất (p = 1) có dạng: Yt b o b1Y( t 1) t (2.3.14) Mô hình tự hồi quy bậc hai (p = 2) có dạng: Yt b o b1Y( t 1) b 2 Y( t 2 ) t (2.3.15) 2.4. XẤP XỈ HÀM VÀ ỨNG DỤNG 2.4.1. Khái niệm chung về xấp xỉ hàm Xấp xỉ hàm (Functional Approximation) là phương pháp tính toán dựa trên nguyên lý bình phương nhỏ nhất để xác định các tham số của một hàm lý thuyết (đã xác định) dựa trên các số liệu thực nghiệm rời rạc, nó còn được gọi là phương pháp ước lượng tham số. Để xấp xỉ hàm giữa hai đại lượng X và Y, trong đó coi Y là hàm của X, trước hết phải xem xét các dãy số liệu xi, yi và lựa chọn hàm số lý thuyết mô tả mối liên hệ (hàm số) giữa hai đại lượng đó dưới dạng: Y f (X) (2.4.1) Trong thực tế có những quy luật có thể biết trước dạng của hàm số f(X) nhưng cũng có nhiều trường hợp không biết được mối quan hệ hàm số của các dãy số liệu thực nghiệm. Trong trường hợp này có thể sử dụng hàm đa thức bậc cao (bậc k) để biểu thị quan hệ cần tìm. Đa thức bậc k có dạng chung là: y i a o a 1 x i a 2 x i2 a 3 x 3i .... a k x ik (2.4.2) Trong đa thức trên, các tham số (hệ số) a o , a 1 , a 2 ....a k cần được xác định theo nguyên lý bình phương nhỏ nhất, khi đó ta đã thiết lập được một hàm lý thuyết để mô phỏng mối quan hệ giữa hai dãy số liệu quan trắc (dãy số liệu 87
- thực nghiệm). Trên hình 2.8 thể hiện các các giá trị quan trắc rời rạc và đồ thị của hàm trong hệ tọa độ XOY. Sử dụng nguyên lý bình phương nhỏ nhất để xác định các tham số nói chung và hàm xấp xỉ nói riêng còn được gọi là phương pháp tính toán làm khớp tốt nhất (Best-Fit Computing) [25]. Cần lưu ý tới một nguyên tắc bắt buộc là dãy số liệu quan trắc phải đủ lớn để có thể xác định được các tham số. Nếu sử dụng hàm đa thức một biến số dạng (2.4.2), số lượng trị quan trắc n trong dãy xi và yi phải thỏa mãn: n k 1 (2.4.3) Hình 2.8. Giá trị quan trắc và đồ thị của hàm xấp xỉ Trong nghiên cứu các quy luật giá trị của các đại lượng phân bố trong không gian (như trọng trường, từ trường, ...) có thể phải xác định (ước lượng) mối quan hệ giữa đại lượng Z với các giá trị tọa độ X và Y dưới dạng hàm số sau: Z g ( X, Y ) (2.4.4) Nếu chưa xác định được quan hệ hàm g(X,Y) nói trên, có thể lựa chọn hàm đa thức bậc cao đối với X, Y để mô tả quy luật giá trị của Z. Hàm đa thức bậc cao có dạng: z i b 0 b1 x i b 2 y i b 3 x i2 b 4 y i2 b 5 x i y i ..... (2.4.5) Các tham số b o , b1 , b 2 , b 3 ..... cũng được xác định theo nguyên lý bình phương nhỏ nhất. 88
- Nếu áp dụng dạng hàm đa thức 2 biến (2.4.5), số giá trị quan trắc phải lớn hơn hoặc bằng số lượng tham số bi có trong đa thức đó. Chất lượng của xấp xỉ hàm phụ thuộc vào các yếu tố sau: - Dạng hàm số lựa chọn (phù hợp hay không phù hợp). - Số lượng trị quan trắc trong các dãy số liệu. - Phân bố của các giá trị quan trắc. - Độ chính xác của các giá trị quan trắc. 2.4.2. Phương pháp tính xấp xỉ hàm Mục tiêu của phương pháp xấp xỉ hàm là xác định các tham số của một hàm số (lý thuyết) dựa trên các số liệu thực nghiệm rời rạc. Dạng của hàm số lý thuyết này có thể xác định trước (chọn trước), trong đó có các tham số cần xác định. Cũng có thể gặp trường hợp chưa xác định được dạng của hàm số mà phải sử dụng một hàm đa thức bậc cao để làm hàm xấp xỉ với quy luật của các dãy số liệu thực nghiệm. Các tham số của hàm xấp xỉ có vai trò trong dự báo, nội suy hoặc làm trơn kết quả quan trắc. 2.4.2.1. Trường hợp chưa biết trước dạng hàm số Trong một số trường hợp, chúng ta có các dãy số liệu quan trắc nhưng chưa biết mối quan hệ giữa chúng. Trong trường hợp này cần phải xác định mối quan hệ giữa các dãy số liệu để có thể biểu diễn giá trị dãy số liệu Y như là các giá trị (gần đúng) của một hàm toán học tương ứng với các biến số X. Trong đó phải xác định các tham số của hàm. Để đơn giản, xét trường hợp xác định mối quan hệ giữa hai dãy số liệu quan trắc X và Y, với n cặp giá trị dạng số. Trước hết cần áp dụng phương pháp phân tích hồi quy để đánh giá một vài tính chất của hai dãy số liệu đó. Các tính toán phân tích hồi quy gồm: - Tính hệ số tương quan thực nghiệm theo công thức (2.3.4). Nếu như trị tuyệt đối của hệ số tương quan thực nghiệm xy có giá trị xấp xỉ 1 thì dạng hàm số cần xác định sẽ là hàm tuyến tính: Y a. b.X (2.4.6) Việc xác định 2 tham số a,b của hàm tuyến tính (2.4.6) khá dễ dàng. Nếu trị tuyệt đối của hệ số tương quan thực nghiệm xy có giá trị khác với 1 thì mối quan hệ giữa X và Y không thể là dạng tuyến tính, trong trường hợp này có thể chọn hàm đa thức bậc k dạng biểu thức (2.4.2) để làm hàm xấp xỉ. 89
- Một cách tổng quát, trong trường hợp này chúng ta phải xác định các tham số của một đa thức có thể là hàm bậc nhất (hàm tuyến tính, k = 1) hoặc là đa thức bậc cao (k > 1). y i a 0 a 1 x i a 2 x i2 a 3 x 3i .... a k x ik (2.4.7) trong đó có k + 1 ẩn số cần xác định là ao, a1, a2,....,ak. Từ (2.4.7) lập được phương trình số hiệu chỉnh (phương trình sai số): v i a 0 a 1 x i a 2 x i2 a 3 x 3i .... a k x ik y i (2.4.8) Các phương trình số hiệu chỉnh (2.4.8) được viết ở dạng ma trận: V A.X L (2.4.9) trong đó: 1 x1 x 12 ... x 1k a o v1 y1 1 x2 x 22 ... x k2 a v y A 1 x3 x 32 k ... x 3 ; X 1 ; V 2 ; L 2 (2.4.10) .. .. .. .. .. .. .. .. ak v v y 1 xn x 2n .. x kn n Khi số lượng phương trình số hiệu chỉnh n nhiều hơn số lượng ẩn số, tức là n > k + 1 và coi các phương trình (2.4.8) cùng độ chính xác, theo nguyên lý [vv] = min, sẽ lập được hệ phương trình chuẩn: A T A.X A T L 0 (2.4.11) Giải hệ phương trình chuẩn, sẽ nhận được véc tơ ẩn số X: X ( A T A ) 1 A T L (2.4.12) Sau khi xác định được các tham số của hàm xấp xỉ là ao, a1, a2....ak , chúng ta sẽ nhận được giá trị ước lượng của hàm ứng với biến x i là y i được xác định theo công thức: y i a 0 a 1 x i a 2 x i2 a 3 x 3i .... a k x ik (2.4.13) Giữa giá trị ước lượng y i , giá trị quan trắc y i và số hiệu chỉnh v i có mối quan hệ sau: yi vi yi (2.4.14) Để đánh giá độ chính xác các tham số, chúng ta áp dụng công thức tính sai số trung phương của ẩn số: 90
- m ai Q ai (2.4.15) trong đó Qai là phần tử trên đường chéo của ma trận nghịch đảo Q (A T A) 1 tương ứng với tham số ai, là sai số trung phương đơn vị trọng số được tính theo công thức: vv (2.4.16) n k 1 Độ lớn của giá trị phản ánh chất lượng (mức phù hợp) của bài toán xấp xỉ hàm. 2.4.2.2. Trường hợp đã có trước dạng hàm số Trong trường hợp giữa hai dãy số liệu X và Y đã biết mối quan hệ dưới dạng một hàm số có các tham số chưa xác định: Y=f(X,a,b,c) (2.4.17) trong đó a,b,c là các tham số chưa xác định của hàm (trường hợp này chỉ xét 3 tham số, theo đó có thể xét cho các trường hợp khác) Nhiệm vụ là phải xác định các tham số a, b, c của hàm dựa trên chuỗi số liệu quan trắc của biến X và hàm Y là x i và y i với i = 1, 2,...,n. Trước hết phải xác định trị gần đúng của các tham số a, b, c. Ký hiệu trị gần đúng của các tham số a, b, c là ao, bo, co và số hiệu chỉnh tương ứng là da, db, dc. Giữa chúng có quan hệ: a = ao + da b = bo + db (2.4.18) c = co + dc Từ hàm (2.4.17) và mối quan hệ (2.4.18) ta có thể viết: Y f (X, a o da , b o db, c o dc) (2.4.19) Nếu các trị gần đúng ao, bo và co được xác định rất gần với giá trị cần xác định, khi đó da, db, dc là những giá trị nhỏ, từ đó áp dụng khai triển chuỗi Taylor để đưa (2.4.19) về dạng phương trình số hiệu chỉnh sau: Y Y Y vi da db dc l i (2.4.20) a 0 b 0 c 0 trong đó số hạng tự do li được tính: li f (x i , a 0 , b 0 , c 0 ) y i (2.4.21) 91
- Y Y Y Nếu ký hiệu: i ; i ; i (2.4.22) a 0 b 0 c 0 phương trình số hiệu chỉnh (2.4.20) được viết: v i i da i db i dc l i (2.4.23) Hoặc viết ở dạng ma trận: V A.X L (2.4.24) trong đó: v1 1 1 1 l1 v da 2 2 2 l V ; A 2 ; X db ; L 2 (2.4.25) .. .. .. .. .. dc v n n n n l n Nếu số giá trị quan trắc (n) lớn hơn 3 (n > 3), ta sẽ xác định da, db, dc theo nguyên lý bình phương nhỏ nhất. Hệ phương trình chuẩn sẽ là: A T A.X A T L 0 (2.4.26) Từ đó ta có lời giải hệ phương trình chuẩn: X (A T A) 1 A T L (2.4.27) Việc giải hệ phương trình chuẩn (2.4.6) và tính các tham số a, b, c của hàm xấp xỉ được thực hiện theo quy trình tính toán như đã nêu ở phần trước. Việc tính toán xác định các tham số a, b, c phải tính lặp (nhích dần) cho đến khi trị tuyệt đối của các số hiệu chỉnh da, db, dc khá nhỏ, không lớn hơn (gọi là điều kiện dừng lặp). Việc xác định giá trị gần đúng đầu tiên của các tham số a, b, c không ảnh hưởng đến kết quả xác định tham số mà chỉ ảnh hưởng đến số lần tính lặp. Độ chính xác của tham số phụ thuộc vào giá trị được chọn khi tính lặp. 2.4.3. Ứng dụng của phương pháp xấp xỉ hàm Phương pháp xấp xỉ hàm được sử dụng trong nghiên cứu các quy luật tự nhiên dựa trên các số liệu do đạc hay quan trắc. Các giá trị của hàm được xác định theo phương pháp xấp xỉ (theo nguyên lý bình phương nhỏ nhất) là cơ sở có độ tin cậy cần thiết để sử dụng cho tính toán nội suy, mô phỏng hoặc dự báo các quy luật tự nhiên đó. Trong trắc địa, phương pháp xấp xỉ hàm được ứng dụng trong xử lý các số liệu đo đạc như: - Xác định các tham số của hàm nội suy hoặc hàm hiệp phương sai, hàm bán phương sai trong tính toán nội suy. 92
CÓ THỂ BẠN MUỐN DOWNLOAD
-
Giáo trình Xử lý nước thải
322 p | 1822 | 666
-
Giáo trình Xử lý số liệu thực nghiệm - TS. Mai Xuân Trung
154 p | 955 | 212
-
Giáo trình: Xử lý bức xạ và cơ sở của công nghệ bức xạ
98 p | 196 | 52
-
Kỹ thuật xử lý số liệu bằng thống kê toán học trên máy tính: Phần 1
89 p | 191 | 45
-
Xử lý số liệu thực nghiệm - TS Mai Xuan Trung
154 p | 212 | 25
-
Giáo trình Xử lý số liệu thực nghiệm bằng Toán học thống kê - TS Tạ Thị Thảo
93 p | 241 | 18
-
Giáo trình Xử lý nước thải - TS. Nguyễn Trung Việt
83 p | 98 | 17
-
Giáo trình Hệ thống định vị toàn cầu GPS và đo GPS (Nghề: Trắc địa công trình - Cao đẳng) - Trường Cao đẳng nghề Xây dựng
78 p | 27 | 9
-
Giáo trình Xử lý bức xạ và cơ sở của công nghệ bức xạ - GS. TS. Trần Đại Nghiệp Phần 5
0 p | 95 | 8
-
Giáo trình Định vị vệ tinh: Phần 2 - Trường ĐH Công nghiệp Quảng Ninh
87 p | 19 | 8
-
Giáo trình Xử lý bức xạ và cơ sở của công nghệ bức xạ - GS. TS. Trần Đại Nghiệp Phần 1
0 p | 89 | 6
-
Giáo trình Xử lý số liệu trắc địa: Phần 1 - PGS.TS Đặng Nam Chinh (Chủ biên)
72 p | 50 | 6
-
Giáo trình Thực tập trắc địa cơ sở: Phần 2 - Trường ĐH Tài nguyên và Môi trường Hà Nội
68 p | 16 | 5
-
Giáo trình Xử lý số liệu khí tượng và dự báo thời tiết bằng phương pháp thống kê vật lý: Phần 2
59 p | 32 | 5
-
Giáo trình Xử lý số liệu khí tượng và dự báo thời tiết bằng phương pháp thống kê vật lý: Phần 1
76 p | 31 | 5
-
Khảo sát một số phương pháp xác định vị trí điểm giao cắt trong xử lý số liệu đo cao vệ tinh
4 p | 88 | 4
-
Giáo trình Xử lý nước cấp (Ngành: Công nghệ kỹ thuật tài nguyên nước - Cao đẳng) - Trường Cao đẳng Xây dựng số 1
122 p | 5 | 2
Chịu trách nhiệm nội dung:
Nguyễn Công Hà - Giám đốc Công ty TNHH TÀI LIỆU TRỰC TUYẾN VI NA
LIÊN HỆ
Địa chỉ: P402, 54A Nơ Trang Long, Phường 14, Q.Bình Thạnh, TP.HCM
Hotline: 093 303 0098
Email: support@tailieu.vn