逐年海温资料滤波问题

麦田_smile

我的初始资料为1850-2012共163年的逐年全球海温资料（180*89，单位为K），希望得到1-10a的波段。下面贴出程序及滤波前后的对比：
1、程序：
parameter(nx=180,ny=89,nt=163,mw1=1,mw2=10)
      dimension sst(nx,ny,nt),sst1(nx,ny,nt)
      dimension cx(nt),xy(nt),bx(nt),tt(nt)
c-----------input data------------------
      open(11,file='E:\1\78\6\ocean_cz_annual.grd',form='binary')
      read(11) (((sst(i,j,k),i=1,nx),j=1,ny),k=1,nt)
close(11)
      do 70 i=1,nx
         do 70 j=1,ny
      if (sst(i,j,1).gt.10000)then
        do k=1,nt
        sst1(i,j,k)=-9999.
        enddo
  else
   do it=1,nt
        tt(it)=sst(i,j,it)
      call filter(nt,tt,bx,cx,xy,mw1,mw2)
      sst1(i,j,it)=bx(it)
      enddo
endif
  70  continue
      open(12,file='E:\1\78\6\ocean_cz_annual_filter.grd',form='binary')  
      write(12) (((sst1(i,j,it),i=1,nx),j=1,ny),it=1,nt)
close(12)
      end
      subroutine filter(n,x,bx,cx,xy,mw1,mw2)
      dimension x(n),bx(n),cx(n),xy(n)
      w1=2*3.1415926/float(mw1)
      w2=2*3.1415926/float(mw2)
      w0=sqrt(w1*w2)
      dt=1.0
      ww=2*abs(sin(w1*dt)/(1+cos(w1*dt))
     $ -sin(w2*dt)/(1+cos(w2*dt)))
      ww0=4*sin(w1*dt)*sin(w2*dt)/((1+cos(w1*dt))
     $ *(1+cos(w2*dt)))
      a0=2*ww/(4+2*ww+ww0)
      b1=2*(ww0-4)/(4+2*ww+ww0)
      b2=(4-2*ww+ww0)/(4+2*ww+ww0)
      do 89 i0=1,n-2
      i=i0+2
      bx(1)=0.0
      bx(2)=0.0
      cx(i)=a0*(x(i)-x(i-2))-b1*bx(i-1)-b2*bx(i-2)
  89  bx(i)=cx(i)
      do 91 i=1,n
  91  x(i)=bx(n+1-i)
      do 92 i0=1,n-2
      i=i0+2
      xy(1)=x(1)
      xy(2)=x(2)
      cx(i)=a0*(x(i)-x(i-2))-b1*xy(i-1)-b2*xy(i-2)
  92  xy(i)=cx(i)
      do 93 i=1,n
  93  bx(i)=xy(n+1-i)
      end

2、滤波前后对比
前
  
后


为什么会出现这么离谱的结果呢？想了三四天还是没能想清楚其中哪一环节出错了。只好来求助各位了

麦田_smile

最后一张图是减去了273.16的滤波前的海温分布

lqouc

这个问题你弄得太复杂了吧，你不就是想要个年际变化么？
都把谐波滤波整出来了。完全没必要吧。你这个程序我就不仔细看了。
其实如果你要求不是很精确的话，拿你的原始序列减去9点或者11点滑动平均序列就好了吧。应该还是比较好的，你可以用方差验证下。
至于编程的话滑动平均什么的在grads里面就三四行的代码，不用你费这么大劲。

麦田_smile

lqouc 发表于 2013-9-15 16:29

                               
登录/注册后可看大图

这个问题你弄得太复杂了吧，你不就是想要个年际变化么？
都把谐波滤波整出来了。完全没必要吧。你这个程序 ...

非常谢谢版主的回复 我先试下

麦田_smile

lqouc 发表于 2013-9-15 16:29

                               
登录/注册后可看大图

这个问题你弄得太复杂了吧，你不就是想要个年际变化么？
都把谐波滤波整出来了。完全没必要吧。你这个程序 ...

还想请问下您 方差验证是怎么回事，原谅我，不太懂

lqouc

就是你看看滤波后两个部分的方差之和跟原始的方差量级一样不一样。这个一般不会出错，除非你的程序编错了。这也不是一个验证方法，只是为了证明你的程序没弄错。不过这么简单的程序，谁会错呢。
你真的想验证的话不妨用功率谱或者小波分析。

麦田_smile

lqouc 发表于 2013-9-15 18:13

                               
登录/注册后可看大图

就是你看看滤波后两个部分的方差之和跟原始的方差量级一样不一样。这个一般不会出错，除非你的程序编错了。 ...

好的 谢谢您

大力水手

楼主的全球海温资料在哪里找到的啊，可否分享？谢谢。

麦田_smile

大力水手 发表于 2013-9-15 23:59

                               
登录/注册后可看大图

楼主的全球海温资料在哪里找到的啊，可否分享？谢谢。

是模式转出来的，你要？

大力水手

麦田_smile 发表于 2013-9-16 09:14

                               
登录/注册后可看大图

是模式转出来的，你要？

恳请楼主发给我一份，我的电邮：zkzl@21cn.com