地球物理反演作业

Posted on | In

作业1

fortran计算：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
program ch1_filter
implicit none
real, dimension(:), allocatable :: Xin, Xout1_average, Xout2_average,  &
Xout1_weighted, Xout2_weighted, Xout1_middle, Xout2_middle
integer :: n=0, i,j, STAT
real :: wn
integer, parameter :: window1=7, window2=13
real, dimension(window1) :: w1, input1
real, dimension(window2) :: w2, input2


open(1, file='homework.txt', form='formatted', status='old' )
! 读取文件行数
DO WHILE (.TRUE.)
    READ(1,*,IOSTAT=STAT)
    IF(STAT .NE. 0) EXIT
    n=n+1
ENDDO

open(2, status='replace',action='write', file='result.dat')
allocate(Xin(n), Xout1_average(int(n-2*(window1-1)/2)),Xout2_average(int(n-2*(window2-1)/2)))
allocate(Xout1_weighted(n-2*(window1-1)/2),Xout2_weighted(n-2*(window2-1)/2))
allocate(Xout1_middle(n-2*(window1-1)/2),Xout2_middle(n-2*(window2-1)/2))
close(unit=1)


open(3, file='homework.txt', form='formatted', status='old' )
read(3,*) (Xin(i) , i=1,n) 

do i=1,n
print *, Xin(i)
end do

! 平均处理
do i=1, n-2*(window1-1)/2
    Xout1_average(i)=0.
    do j=i, i+window1-1
        Xout1_average(i)=Xout1_average(i)+Xin(j)
    enddo
    Xout1_average(i)=Xout1_average(i)/window1
enddo

do i=1, n-2*(window2-1)/2
    Xout2_average(i)=0.
    do j=i, i+window2-1
        Xout2_average(i)=Xout2_average(i)+Xin(j)
    enddo
    Xout2_average(i)=Xout2_average(i)/window2
enddo

! 加权平均1
do i=1, window1
    wn=(window1+1)/2
    if(i.le.wn) then 
    w1(i) = i/(wn**2)
    else
    w1(i) = (2*wn-i)/(wn**2)
    endif
enddo

do i=1, n-2*(window1-1)/2
    Xout1_weighted(i)=0.
    do j=i, i+window1-1
        Xout1_weighted(i)=Xout1_weighted(i)+Xin(j)*w1(j-i+1)
    enddo
enddo
! 加权平均2
do i=1, window2
    wn=(window2+1)/2
    if(i.le.wn) then 
    w2(i) = i/(wn**2)
    else
    w2(i) = (2*wn-i)/(wn**2)
    endif
enddo

do i=1, n-2*(window1-1)/2
    Xout2_weighted(i)=0.
    do j=i, i+window1-1
        Xout2_weighted(i)=Xout2_weighted(i)+Xin(j)*w2(j-i+1)
    enddo
enddo

! 中值滤波1
do i=1, n-2*(window1-1)/2

    do j=i, i+window1-1
        input1(j-i+1) = Xin(j)
    enddo
    call sort(input1,Xout1_middle(i),window1)
enddo

! 中值滤波2
do i=1, n-2*(window2-1)/2
    
    do j=i, i+window2-1
        input2(j-i+1) = Xin(j)
    enddo
    call sort(input2,Xout2_middle(i),window2)
enddo


do i=1, 132
    write(2, 100) i, Xin(i+10), Xout1_average(i+7), Xout2_average(i+4), Xout1_weighted(i+7), Xout2_weighted(i+4), &
            Xout1_middle(i+7), Xout2_middle(i+4)
enddo
 100 format(I4, 7(x, F15.4))
end program ch1_filter

SUBROUTINE sort(a,b,n)
dimension :: a(n)
real b
real temp
do i=1,n-1
  do j=i+1,n
   if (a(i).gt.a(j)) then
   temp = a(i)
   a(i) = a(j)
   a(j) = temp
   endif
  enddo
enddo 
b=a((n+1)/2)
return
end

matlab绘制图像

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
load result.dat

figure(1)
hold on;
subplot(4,2,1:2)
plot(result(:,1), result(:,2), 'r-')
title('原始数据')
subplot(4,2,3)
plot(result(:,1), result(:,3), 'r-')
title('平均值，窗口7')
subplot(4,2,4)
plot(result(:,1), result(:,4), 'r-')
title('平均值，窗口13')
subplot(4,2,5)
plot(result(:,1), result(:,5), 'r-')
title('权重，窗口7')
subplot(4,2,6)
plot(result(:,1), result(:,6), 'r-')
title('权重，窗口13')
subplot(4,2,7)
plot(result(:,1), result(:,7), 'r-')
title('中值滤波，窗口7')
subplot(4,2,8)
plot(result(:,1), result(:,8), 'r-')
title('中值滤波，窗口13')