Tài liệu Bài giảng xử lý số liệu địa vật lý - phan thiên hương

Tài liệu tham khảo 1. Bài giảng cơ sở lý thuyết XLSL ĐVL, Phạm Năng Vũ, 2002 2. XLSLĐVL, Phạm Năng Vũ, Nguyễn Huy Ngọc, 1997 3. Digital signal processing, Proaskis J.G.; Manolakis D.G., 1996 4. Fundamentals of Geophysical Data processing, Claerbout J.F, 1976 Xử lý số liệu Địa vật lý 5. Spectral analysis and filter theory in applied Geophysics, Buttkus B., 2000 6. Time series analysis and inverse theory for geophysics, 2004 7. Seismic data processing, Ylmaz O., Doherty S.M., 1987 Phan Thiên Hương-2013 Phan Thiên Hương-2013 Phan Thiên Hương-2013 Phan Thiên Hương-2013 Mở đầu • Vai trò của XLSL • 1a) mô hình địa chất theo các vết lộ 0.1b và 0.1c mô hình địa chất suy ra từ tài liệu địa vật lý. • Tuy nhiên, khi nghiên cứu một cách gián tiếp thì các phương pháp địa vật lý sẽ có những nhược điểm như tính đa trị của kết quả (hình 1a và 1b), sự không phù hợp giữa các mô hình tính toán đơn giản với thực tế địa chất phức tạp, kết quả có thể phụ thuộc vào nhiều yếu tố như độ nhạy, độ ổn định và độ chính xác của máy móc. Phan Thiên Hương-2013 Phan Thiên Hương-2013 1 Phan Thiên Hương-2013 Phan Thiên Hương-2013 Phan Thiên Hương-2013 Phan Thiên Hương-2013 Môi trường Địa chất Trường ĐVL Số liệu ĐVL Minh giải Xử lý Nguồn Kết quả Phân tích Phan Thiên Hương-2013 Phan Thiên Hương-2013 2 • Số liệu • Tín hiệu • Nhiễu n Fj = ∑ fi i =1 Trường vật lý đo được là trường tổng gồm các phần đóng góp của nhiều đối n tượng tạo ra trường lên điểm quan sát. Tại điểm j hoặc thời điểm j nào đó ta Fj = ∑ fi có giá trị trường quan sát được Fj i =1 fi là giá trị của phần trường do đối tượng thứ i gửi về điểm quan sát thứ j. Fj = s j + n j Phan Thiên Hương-2013 Phan Thiên Hương-2013 Nội dung của khóa học này là giới thiệu các thuật toán dùng trong xử lý giúp làm Môi trường Địa chất Trường ĐVL Số liệu ĐVL rõ những thông tin có ích hoặc làm cho quá trình minh giải sau đó trở nên dễ dàng hơn. Nói cách khác 2 nhiệm vụ chính được đề cập trong giáo trình này là: •Xác định tín hiệu đo được có chứa tín hiệu có ích không? Nguồn •Nếu tồn tại tín hiệu thì phải tách tín hiệu ra khỏi giá trị quan sát Khóa học này sẽ gồm các phần chính sau: •Số hóa Kết quả Phân tích Minh giải Xử lý •Biến đổi Fourier và Z •Lọc tuyến tính: - lọc không tối ưu • - lọc tối ưu •Phát hiện tín hiệu yếu: áp dụng lý thuyết toán xác suất, thống kê, các hàm ngẫu nhiên để xây dựng các chỉ tiêu định nghiệm thống kê. Dựa trên các định nghiệm thống kê này để xác định tồn tại tín hiệu hay không tại vị trí quan sát. Phan Thiên Hương-2013 Phan Thiên Hương-2013 3 Đường analog Chương I : Số hóa tín hiệu 1. Tín hiệu và xử lý số liệu: • Tín hiệu • nguồn tín hiệu. • hệ thống thiết bị • xử lý số liệu 2. Tính chất của tín hiệu • Tín hiệu có thể phân làm 2 loại: tín hiệu liên tục và tín hiệu rời rạc • Tín hiệu ngẫu nhiên • Tần số của tín hiệu 3. Số hóa hàm liên tục ( analog to digital) Phan Thiên Hương- 2013 Phan Thiên Hương- 2013 Tín hiệu và xử lý số liệu Tín hiệu và xử lý số liệu Tín hiệu : Tín hiệu : • là một đại lượng vật lý thay đổi theo thời gian, không gian hoặc bất kỳ một hoặc nhiều bíến số khác. • là một đại lượng vật lý thay đổi theo thời gian, không gian hoặc bất kỳ một hoặc nhiều bíến số khác. • Trong toán học có thể biểu diễn tín hiệu như hàm của một hoặc nhiều biến độc lập. Thí dụ như s (t ) = 7 t s( x, y) = 3x + 2 xy + 10 y 2 Phan Thiên Hương- 2013 Tín hiệu và xử lý số liệu Phan Thiên Hương- 2013 Tín hiệu và xử lý số liệu Tín hiệu : • là một đại lượng vật lý thay đổi theo thời gian, không gian hoặc bất kỳ một hoặc nhiều bíến số khác. • Trong toán học có thể biểu diễn tín hiệu như hàm của một hoặc nhiều biến độc lập. Thí dụ như s (t ) = 7t s( x, y ) = 3x + 2 xy + 10 y 2 • Trong thực tế Tín hiệu lời nói phụ thuộc biên độ, tần số và pha N ∑ A (t ) sin[2πF (t )t + θ (t )] i i i i =1 Phan Thiên Hương- 2013 Phan Thiên Hương- 2013 1 Tín hiệu và xử lý số liệu Tín hiệu và xử lý số liệu nguồn tín hiệu. nguồn tín hiệu. •Tín hiệu tự nhiên hệ thống thiết bị (system). •Tín hiệu nhân tạo Thí dụ như các bộ lọc. Phan Thiên Hương- 2013 Tín hiệu và xử lý số liệu Phan Thiên Hương- 2013 Tín hiệu và xử lý số liệu hệ thống thiết bị nguồn tín hiệu. hệ thống thiết bị (system). Thí dụ như các bộ phận thu phát tín hiệu, máy tính. • Môi trường địa chất • Trường Địa vật lý xử lý số liệu. Thí dụ như ta dùng bộ lọc để loại bỏ nhiễu đồng nghĩa ta đã tiến hành xử lý số liệu. Các phần mềm ứng dụng như geoframe, vista, coscad,… Phan Thiên Hương- 2013 Tín hiệu và xử lý số liệu Phan Thiên Hương- 2013 Tín hiệu và xử lý số liệu Xử lý số liệu Phan Thiên Hương- 2013 Phan Thiên Hương- 2013 2 Tính chất của tín hiệu Tính chất của tín hiệu Tín hiệu có thể phân làm 2 loại: liên tục và rời rạc Tín hiệu : xác định ngẫu nhiên: Phan Thiên Hương- 2013 Tính chất của tín hiệu Phan Thiên Hương- 2013 Tính chất của tín hiệu Tần số của tín hiệu liên tục và rời rạc A1: Ứng với một giá trị F cho trước, xa(t) là hàm tuần hoàn. xa (t ) = A cos(Ωt + θ ) −∞ 2 Fmax: FN = Fmax FN –Nyquist rate (tỷ số . Phan Thiên Hương- 2013 Tín hiệu và xử lý số liệu Phan Thiên Hương- 2013 Tính chất của tín hiệu Tín hiệu : Phan Thiên Hương- 2013 xác định ngẫu nhiên: Phan Thiên Hương- 2013 Tín hiệu và xử lý số liệu Phan Thiên Hương- 2013 6 CHƯƠNG II §1.Thuật toán biến đổi Fourier của hàm tuần hoàn PHÂN TÁCH TÍN HIỆU ĐỊA VẬT LÝ §2. Thuật toán biến đổi Fourier của hàm không tuần hoàn §3 Đặc trưng của phổ : §4. Biến đổi Fourier của tín hiệu rời rạc : Biến đổi Fourier §5. Biến đổi Fourier 2 chiều §6 Hàm Dirac Delta DVL- K54 Phan Thiên Hương DVL- K54 Phan Thiên Hương DVL- K54 Phan Thiên Hương DVL- K54 Phan Thiên Hương DVL- K54 Phan Thiên Hương DVL- K54 Phan Thiên Hương biến đổi Fourier • tín hiệu có thể phân tích thành tổng của các hàm sin ( tần số khác nhau). • Thuật toán dùng để phân tích tín hiệu thành các hàm sin này được gọi là biến đổi Fourier 1 Tín hiệu khác nhau bởi biên độ, tần số và pha. DVL- K54 Phan Thiên Hương §1.Thuật toán biến đổi Fourier của hàm tuần hoàn DVL- K54 Phan Thiên Hương §1.Thuật toán biến đổi Fourier của hàm tuần hoàn Một hàm tuần hoàn x(t) có thể đựoc chuyển đổi sang chuỗi Fourier nếu nó thỏa mãn những điều kiện sau: 1. Là hàm đơn trị, tuần hoàn của biến t với chu kỳ T Một hàm tuần hoàn x(t) có thể đựoc chuyển đổi sang chuỗi Fourier nếu nó thỏa mãn những điều kiện sau: 1. Là hàm đơn trị, tuần hoàn của biến t với chu kỳ T 2. Thỏa mãn điều kiện Dirichlet (có nghĩa là hàm liên tục loại trừ một số hữu hạn điểm không liên tục và chỉ có một hữu hạn cực đại và cực tiểu) DVL- K54 Phan Thiên Hương §1.Thuật toán biến đổi Fourier của hàm tuần hoàn DVL- K54 Phan Thiên Hương §1.Thuật toán biến đổi Fourier của hàm tuần hoàn x(t) có thể được biểu diễn dưới dạng chuỗi Fuorier như sau: Một hàm tuần hoàn x(t) có thể đựoc chuyển đổi sang chuỗi Fourier nếu nó thỏa mãn những điều kiện sau: 1. Là hàm đơn trị, tuần hoàn của biến t với chu kỳ T 2. Thỏa mãn điều kiện Dirichlet (có nghĩa là hàm liên tục loại trừ một số hữu hạn điểm không liên tục và chỉ có một hữu hạn cực Phép biến đổi này được gọi là biến đổi Fourier. 1 Với f 0 = T a0 giá trị trung bình của tín hiệu trên chu kỳ T T DVL- K54 n =1 Ở đây f0 hay ω0=2π/T là tần số cơ bản ứng với chu kỳ T. đại và cực tiểu) 3. Là hàm có giới hạn, có nghĩa ∞ x (t ) = a 0 + ∑ (a n cos( 2πnf 0 t ) + bn sin( 2πnf 0 t )) ∫ x(t ) dt ≤ c < ∞ 0 Phan Thiên Hương DVL- K54 Phan Thiên Hương 2 §1.Thuật toán biến đổi Fourier của hàm tuần hoàn §1.Thuật toán biến đổi Fourier của hàm tuần hoàn Mặt khác ta có thể biểu diễn chuỗi Fourier theo dạng mũ phức: Với an va bn là hệ số Fourier, được xác định theo x(t) : 1 i 2πnf 0t (e + e −i 2πnf 0t ) 2 1 sin( 2πnf 0t ) = (e i 2πnf 0t − e −i 2πnf 0t ) 2 cos( 2πnf 0t ) = T an = 2 x(t ) cos(2πnf0t )dt T ∫0 n=0,1,2,... T T 2 bn = ∫ x(t) sin(2πnf0t)dt T0 Xn = n=0,1,2,... 1 x(t )e −i 2πnf0t dt T ∫0 Xn được gọi là phổ của x(t) DVL- K54 Phan Thiên Hương §2. Thuật toán biến đổi Fourier của hàm không tuần hoàn DVL- K54 Phan Thiên Hương §2. Thuật toán biến đổi Fourier của hàm không tuần hoàn Nếu tín hiệu x(t) là hàm không tuần hoàn thì ta không thể biểu diễn nó dưới dạng phổ rời rạc như đối với hàm tuần hoàn. Tuy nhiên nếu hàm x(t) vẫn tuân theo điều kiện Dirichlet trong 1 khoảng bất kỳ nào và tích ∞ X(f ) = ∫ x(t )e − i 2πft dt −∞ phân vẫn hội tụ (converged) nghĩa là tích phân này có thể tính được thì hàm x(t) có thể được biẻu diễn dưới dạng tích phân Với ω = 2πf (radians per second), Fourier: x(t ) = ∞ x(t ) = ∫ X ( f )e i 2πft 1 2π ∞ ∫ X (ω )e iωt dω −∞ df −∞ DVL- K54 Phan Thiên Hương §2. Thuật toán biến đổi Fourier của hàm không tuần hoàn DVL- K54 Phan Thiên Hương DVL- K54 Phan Thiên Hương §2. Thuật toán biến đổi Fourier của hàm không tuần hoàn DVL- K54 Phan Thiên Hương 3 §2. Thuật toán biến đổi Fourier của hàm không tuần hoàn §2. Thuật toán biến đổi Fourier của hàm không tuần hoàn 1 x(t ) =  0 −ε /2 ≤ t ≤ ε /2 t < −ε / 2hayt > ε / 2 ε /2 ∞ X(f ) = ∫ x(t )e −∞ DVL- K54 − i 2πft dt = ∫ −ε / 2 e −i 2πft dt = e i 2πf ε −e i 2πf 2 −i 2πf ε 2 =ε sin( 2πf 2πf ε ε 2 ) Phan Thiên Hương DVL- K54 Phan Thiên 2 Hương Phan Thiên Hương DVL- K54 Phan Thiên Hương §3 Đặc trưng của phổ : Tín hiệu khác nhau bởi biên độ, tần số và pha. DVL- K54 §3 Đặc trưng của phổ : §3 Đặc trưng của phổ : Phổ được biểu diễn dưới dạng số phức : X(f)=U(f)+iV(f) Trong đó U(f) là phần thực còn V(f) là phần ảo. Tương đương với : X(f) = |X(f)| eiθ(f) Với |X(f)| = (U2(f) + V2(f))1/2 Và θ(f) = arctan (V(f)/U(f)) khi U(f)≠0 trong khoảng -π đến π |X(f)| - biên độ của phổ (amplitude spectrum of x(t)) of amplitude spectral density) and θ(f)- pha của phổ (phase spectrum of x(t)) (còn gọi là phổ biên độ và phổ pha) DVL- K54 Phan Thiên Hương DVL- K54 Phan Thiên Hương 4 §3 Đặc trưng của phổ : §3 Đặc trưng của phổ : Nếu x(t) trong miền thời gian là tổng của 2 tín hiệu x(t ) = s1 (t ) + s 2 (t ) Thì trong miền tần số ta có : X ( f ) e iθ x ( f ) = S1 ( f ) e iθ1 ( f ) + S 2 ( f ) e iθ 2 ( f ) Khi đó phổ biên độ : Và phổ pha : DVL- K54 Phan Thiên Hương DVL- K54 Phan Thiên Hương §4. Biến đổi Fourier của tín hiệu rời rạc : §3 Đặc trưng của phổ : Độ rộng của phổ Xi = {X-r, X-r+1,......X-i,...X-1,X0, X1,.....Xi,.....Xr-1} r −1 X i = A0 + 2∑ Am cos 2πmf1t + Bm sin 2πmf 0 t m =1 DVL- K54 Phan Thiên Hương DVL- K54 Phan Thiên Hương §4. Biến đổi Fourier của tín hiệu rời rạc : §4. Biến đổi Fourier của tín hiệu rời rạc : Am, Bm –hệ số Fourier, chúng là các giá trị có thể xác định được trên cơ sở thỏa mãn đẳng thức (4.1) Am = 1 r −1 mi ∑ X i cos(2π n ) n i =− r Bm = 1 r −1 mi ∑ X i sin(2π n ) n i =− r Với n=2r |X(f)| = (Am2 +Bm2)1/2 θ(f) = arctg(-Bm/Am) DVL- K54 Phan Thiên Hương DVL- K54 Phan Thiên Hương 5 §5. Biến đổi Fourier 2 chiều §6 Hàm Dirac Delta Biến đổi Fourier còn có thể dùng cho hàm 2 biến. Gọi g(x,y) là hàm không tuần hoàn của 2 biến x,y. Nếu Thông thường ta có tại giá trị t=t0 thì hàm x(t)=x(t0). Đối với hàm Delta ta có : δ(t-t0) =1 for t= 0 δ(t-t0)=0 for t≠0 X(f) = Khi đó biến đổi Fourier (Fourier transform) được : X(f) Và ngược lại : DVL- K54 Phan Thiên Hương DVL- K54 Phan Thiên Hương §6 Hàm Dirac Delta §6 Hàm Dirac Delta DVL- K54 Phan Thiên Hương DVL- K54 Phan Thiên Hương DVL- K54 Phan Thiên Hương 6 §1.Khái niệm chung: §2. Phương pháp thực hiện quá trình lọc. Xung (impulse), sự đáp ứng của xung trong miền thời gian và miền tần số. Chương III §3. Các bộ lọc số tuyến tính không tối ưu LÝ THUYẾT CƠ BẢN CỦA CÁC BỘ LỌC Bộ lọc tần thấp Bộ lọc tần cao Bộ lọc dải Bộ lọc khe (hình chữ V) ĐVL K54- Mỏ Địa chất XLSL ĐVL-Phan Thiên Hương LÝ THUYẾT CƠ BẢN CỦA CÁC BỘ LỌC Định nghia Tóan tử lọc được phân loại theo tính chất của nó: 1. Toán tử lọc được gọi là tuyến tính nếu nó có tính chồng chập (additive) Homogeneous (đồng nhất): 2. Toán tử lọc được gọi là bất biến theo thời gian (timeinvariant) nếu mối quan hệ trong toán tử lọc luôn đúng khi thời gian được dich chuyển với time shift τ: ĐVL K54- Mỏ Địa chất XLSL ĐVL-Phan Thiên Hương Để thực hiện quá trình lọc đòi hỏi phải giải quyết 2 nhiệm vụ ĐVL K54- Mỏ Địa chất XLSL ĐVL-Phan Thiên Hương • Bộ lọc luôn phải thỏa mãn tính ổn định. Bộ lọc tuyến tính, bất biến được gọi là ổn định nếu như tín hiệu đầu vào x(t) bị giới hạn (bounded) thì tín hiệu đầu ra y(t) cũng bị chặn (bounded too) • Nếu như tính ổn định của bộ lọc không được thỏa mãn, thì chúng ta phải đặt điều kiện để nó thỏa mãn tính ổn định. • Toán tử lọc còn được phân lọai causal và acausal (tính nhân quả): Trong trường hợp causal, tín hiệu ra của bộ lọc tại thời gian to chỉ phụ thuộc vào các giá trị của tín hiệu đầu vào tại thời gian t< to ĐVL K54- Mỏ Địa chất XLSL ĐVL-Phan Thiên Hương Trong miền thời gian • Thiết kế các bộ lọc – thuật toán lọc đảm bảo tốt nhất mục đích đặt ra: a) • lọc tốt nhất nhiễu; b) làm méo ít nhất có thể tín hiệu; c) thực hiện nhanh • tín hiệu đầu vào là x(t), gọn, hiệu quả trên máy tính. • tín hiệu đầu ra y(t) Xây dựng các phương pháp thực hiện quá trình lọc trên máy tính. • qua bộ lọc h(t). phương pháp thực hiện quá trình lọc có thể được thực hiện • trong miền thời gian • trong miền tần số. ĐVL K54- Mỏ Địa chất XLSL ĐVL-Phan Thiên Hương ĐVL K54- Mỏ Địa chất XLSL ĐVL-Phan Thiên Hương 1 Trong miền thời gian Trong miền thời gian Đối với bộ lọc tuyến tính, bất biến (time-invariant), tín hiệu đầu ra có thể được biểu diễn dưới dạng rời rạc: y(t) luôn luôn tính đựơc cho tín hiệu đầu vào x(t) khi biết h(t). y j = ∑ hi x j −i xj-i – giá trị rời rạc của tín hiệu đầu vào hi – giá trị rời rạc của hàm lọc hay còn gọi là hàm trọng số. Trong trường hợp nếu ta dùng toán tử lọc chính là hàm Dirac delta 1; i = 0 0; i ≠ 0 δi =  Thì đầu ra chính là hàm xi . ĐVL K54- Mỏ Địa chất XLSL ĐVL-Phan Thiên Hương 2. Trong miền tần số: x (t ) = ĐVL K54- Mỏ Địa chất 1 2π ĐVL K54- Mỏ Địa chất XLSL ĐVL-Phan Thiên Hương Bộ lọc tần thấp ∞ ∫ X ( f )e i 2 πft df −∞ XLSL ĐVL-Phan Thiên Hương ĐVL K54- Mỏ Địa chất XLSL ĐVL-Phan Thiên Hương Bộ lọc tần thấp Bộ lọc tần cao H(ω) Hl Hh fc f Hình 4.5: a)quan hệ giữa bộ lọc tần số thấp và tần số cao; b) Phổ của bộ lọc tần cao ĐVL K54- Mỏ Địa chất XLSL ĐVL-Phan Thiên Hương ĐVL K54- Mỏ Địa chất XLSL ĐVL-Phan Thiên Hương 2 Bộ lọc dải ĐVL K54- Mỏ Địa chất Bộ lọc khe XLSL ĐVL-Phan Thiên Hương ĐVL K54- Mỏ Địa chất XLSL ĐVL-Phan Thiên Hương Trong miền thời gian X = [2, -1, 3, 2, 4]; h = [1, 0, 0]; figure a subplot(3,1,1) stem(X', 'filled'); title('x'); subplot(3,1,2) stem(h', 'filled'); title('h'); nx = length(X); b nh = length(h); ny = nx + nh - 1; y = zeros(1, ny); y=conv(X, h); % convolution of X and h subplot(3,1,3) stem(y, 'filled'); title('x*h'); ĐVL K54- Mỏ Địa chất XLSL ĐVL-Phan Thiên Hương 3 MỘT SỐ KHÁI NIỆM CƠ BẢN BỘ LỌC GIẢI TÍCH MỘT SỐ KHÁI NIỆM CƠ BẢN •BỘ LỌC CHEBYSHEV •BỘ LỌC BUTTERWORTH •BỘ LỌC TUYẾN TÍNH TRONG MIỀN KHÔNG GIAN •THUẬT TOÁN TRUNG BÌNH TRƯỜNG •Dải thông tần (passband): •Dải dốc (slopband-roll off): •Dải cắt (stopband): •Tần số góc (hay còn gọi là tần số cắt) cutoff or corner frequency •Octave OR Decade •dB (decibel) = 10 log10 ( A12 / A22 ) •THUẬT TOÁN NÂNG HẠ TRƯỜNG BỘ LỌC GIẢI TÍCH Bộ lọc chebyshev dB Tỷ số biên độ 1. Biên độ của hàm trọng số (Amplitude of the frequency response Tỷ số năng lượng 1 2 H chebyshev (ω ) = ω  1 + ε C N    ωc  2 -120 10-6 10-12 -80 10-4 10-8 -40 0.01 10-4 BỘ LỌC GIẢI TÍCH ω C N   ωc   đa thức Chebyshev bậc N  C0(ω)=1; C1(ω)= ω; CN+1(ω)+CN(ω)= 2ω CN (ω) 2.Cực trị được tính theo công thức: δ = 1 − (1 − ε 2 ) −1 / 2 BỘ LỌC GIẢI TÍCH Bộ lọc Butterworth Biên độ của hàm trọng số (Amplitude of the frequency response) 1 2 H b (ω ) = ω  1 +    ωc  2n Độ dốc của bộ lọc được tính như sau K (ω ) = 20 lg H (ω 2 ) H (ω1 ) BỘ LỌC GIẢI TÍCH Bộ lọc Butterworth •Thuật toán trung bình trường Bộ lọc tuyến tính trong miền không gian Bộ lọc tuyến tính trong miền không gian • Thuật toán trung bình trường ∆g = H (ω ) = 1 N N ∑ ∆g i =1  ωN  sin    2  ωN 2 i