PRIDE PPP AR源码解读###\src\tedit\read_rinex_file.f90

FORTRAN代码

对检查接收机钟跳解读

!! check receiver clock
  if (.not. use_brdeph .or. pclimit .eq. 0.d0) return
  iepo = 1 !!历元循环
  do while (iepo .le. nepo)
    do i0 = 1, MAXSYS   !!循环系统
      nused = 0
      pmb = 0.d0
      flg = 0
      wpmb = 0.d0
      it = 0
      ij = 0
      !!j为所有系统的卫星编号
      do j = 1, i0-1
        ij = ij + MAXSAT_SYS(j)
      end do
      do j = ij+1, ij+MAXSAT_SYS(i0)
        !!不使用北斗C01-C05
        if ('C' .eq. prn0(j)(1:1)) then
          read (prn0(j) (2:3), '(i2)') prn_int
          if (prn_int .le. 5) cycle
        end if
        !!检核数据 标记为OK可以使用
        if (istrue(flagall(iepo, j), 'ok')) then
          nused = nused + 1  !!记录使用几颗卫星
          pmb(nused) = obs(iepo, j, 6)!!无电离层组合
          flg(nused) = 0  !!当前快标记数组
          wpmb(nused) = 1.d0 !!最小二乘权阵
          it(nused) = j !!记录卫星编号
        end if
      end do
      if (nused .eq. 0) cycle  !!单系统可用卫星为0 不计算
      call get_wgt_mean(.true., pmb, flg, wpmb, nused, 30.d0, k, dt, rms, sig)

      !!当可用卫星>2并且rms大于30.0
      do while (nused - k .gt. 2 .and. rms .gt. 30.d0)
        j = k!!记录删除卫星个数
        call sign_robust(nused, pmb, flg, 60.d0, k)
        if (k .eq. j) exit
        call get_wgt_mean(.false., pmb, flg, wpmb, nused, 30.d0, k, dt, rms, sig)
      end do
      if (k .gt. 0) then
        do j = 1, nused
          if (flg(j) .le. 1) then
            obs(iepo, it(j), 6) = pmb(j) - dt
          end if
          if (flg(j) .ge. 2 .or. dabs(obs(iepo, it(j), 6)) .gt. pclimit) then
            if (dabs(obs(iepo, it(j), 6)) .gt. 2*1.d5) then
              flagall(iepo, it(j)) = set_flag(flagall(iepo, it(j)), 'pc1ms')
            else
              flagall(iepo, it(j)) = set_flag(flagall(iepo, it(j)), 'pcbad')
            end if
          end if
        end do
      end if
    end do
    iepo = iepo + 1
  end do

get_wgt_mean(l_edit,x,flg,wgt,n,mrms,ndel,mean,rms,sig)函数

!! purpose  : compute weighted mean value of an array
!!            & sign-constrained least squares
!! parameter:
!!    input : l_edit -- edit data or not
!!            x,flg  -- data array & flag array
!!            wgt    -- weight array
!!            n      -- # of data
!!    output: ndel -- # of deleted data
!!            mean -- mean value
!!            rms  -- unweighted rms of residuals
!!            sig  -- sigma of mean value
!!
!
subroutine get_wgt_mean(l_edit,x,flg,wgt,n,mrms,ndel,mean,rms,sig)
implicit none

logical*1 l_edit
real*8 x(1:*),sig,mrms,rms,mean,hp,wgt(1:*)
integer*4 flg(1:*)
integer*4 n,ndel

integer*4 i,m,ndel_new
real*8 wgt_sum
!
!! sign-constrained least squares to remove large biases
!..call sign_robust(n,x,flg,ndel)
ndel=0
ndel_new=1
do while(ndel_new.ne.0)
!
!! mean 平均值滤波 flg(i)为2 则删除
  m=0; mean=0.d0; wgt_sum=0.d0
  do i=1,n
    if(flg(i).ge.2) cycle
    mean=mean+wgt(i)*x(i); m=m+1; wgt_sum=wgt_sum+wgt(i)
  enddo
  if(m.ne.0) then
    mean=mean/wgt_sum
  else
    rms=999.d0; sig=999.d0; ndel=n; return
  endif
!
!! sigma & standard deviation 方差 标准差的实时滤波
  sig=0.d0; rms=0.d0
  do i=1,n
    if(flg(i).ge.2) cycle
    rms=rms+(x(i)-mean)**2*wgt(i)
  enddo
 
  if(m.eq.1) then
    ndel=n-1; rms=999.d0; sig=999.d0; return
  else
    rms=dsqrt(rms/(m-1))
    sig=rms/dsqrt(wgt_sum)
  endif
!
!! if achieve goal 设定l_edit 为true编辑数据 否则返回不可使用数据的个数
  if(.not.l_edit.or.3.d0*rms.lt.mrms) then
    ndel=count(flg(1:n).ge.2); return
  endif
!
!! edit data according to residuals 大于3倍rms则删除并设置flg为2
  ndel_new=0; ndel=0
  do i=1,n
    if(flg(i).ge.2) then
      ndel=ndel+1
    else if(dabs(x(i)-mean).ge.3.d0*rms/wgt(i)) then
      ndel=ndel+1
      ndel_new=ndel_new+1
      flg(i)=2
    endif
  enddo
enddo

return
end

!! purpose  : 1-dimensional sign-constrained least squares robust estimation
!!              一维符号约束最小二乘稳健估计
!! parameter:
!!    input : nxl,rxl -- element array 数组元素
!!    output: flg -- flag of each element  每个元素的标志
!!            ndl -- # of deleted
!!
!
subroutine sign_robust(nxl,rxl,flg,drang,ndl)
implicit none

integer*4 nxl,ndl,flg(1:*)
real*8    drang,rxl(1:*)

!
!! local
integer*4 i,j,k,nmed,nmd
real*8 rmed(2),rmd(2),resi(nxl),pnew(nxl),p(nxl),cost_new,cost

!
!! initialize
p = 1.d9 !!1*10的9次方
pnew = 0.d0
!
!! iterate until no elements removed 迭代直到没有元素被删除
if(count(flg(1:nxl).lt.2).gt.2) then
  cost=1.d12
!
!! locate median value(s) of original elements 定位原始元素中位数的值
  call med1(nxl,rxl,flg,nmed,rmed)
  !!循环中位数 可能有一个 或者两个
  do i=1,nmed
!
!! absolute residuals 绝对残差 相对于中位数
    do j=1,nxl
      resi(j)=0.d0
      if(flg(j).ge.2) cycle
      resi(j)=dabs(rxl(j)-rmed(i))
    enddo
!
!! locate median value(s) of absolute residuals 确定绝对残差的中位数
    call med1(nxl,resi,flg,nmd,rmd)
    do j=1,nmd
!
!! flag elements possibly contaminated & compute cost values
!!标记可能被污染的元素并计算成本值
      cost_new=0.d0
      rmd(j)=rmd(j)/0.6745  !!! chen不知道为什么
      do k=1,nxl
        pnew(k)=0.d0
        if(flg(k).ge.2) cycle
        if(resi(k).lt.drang.or.resi(k).le.3.d0*rmd(j)) then
          pnew(k)=1.d0
          cost_new=cost_new+resi(k)**2
        endif
      enddo
      cost_new=cost_new/count(pnew(1:nxl).eq.1.d0)
!
!! save flags with smallest cost value保存成本值最小的标志
      if(cost_new.lt.cost) then
        cost=cost_new
        p(1:nxl)=pnew(1:nxl)
      endif
    enddo
  enddo
!
!! make final decision
  if(count(p(1:nxl).eq.0.d0).gt.count(flg(1:nxl).ge.2)) then
    do i=1,nxl
      if(flg(i).lt.2.and.p(i).eq.0.d0) then
        flg(i)=3
      endif
    enddo
  endif
endif
ndl=count(flg(1:nxl).ge.2)
return
end
!
!! median values 中位数的值
subroutine med1(nxl,rxl,flg,nd,rd)
implicit none
integer*4 i,j,nxl,nd,flg(1:*),ix1
real*8 rxl(1:*),rd(1:*),tmp(nxl),tp
! sort elements from small to large
tmp(1:nxl)=rxl(1:nxl)
do i=1,nxl-1
  if(flg(i).ge.2) cycle
  do j=i+1,nxl
    if(flg(j).ge.2) cycle
    if(tmp(i).gt.tmp(j)) then
      tp    =tmp(i)
      tmp(i)=tmp(j)
      tmp(j)=tp
    endif
  enddo
enddo
! locate median value(s)
j=count(flg(1:nxl).lt.2)
if(mod(j,2).eq.0) then
  nd=2
  rd(1)=tmp(ix1(j/2,nxl,flg))
  rd(2)=tmp(ix1(j/2+1,nxl,flg))
else
  nd=1
  rd(1)=tmp(ix1(j/2+1,nxl,flg))
endif
return
end
!
!! return ith good data
integer*4 function ix1(ind,nxl,flg)
implicit none

integer*4 ind,nxl,flg(1:*),i
i=0
do ix1=1,nxl
  if(flg(ix1).ge.2) cycle
  i=i+1
  if(i.eq.ind) return
enddo
write(*,'(a)') '***ERROR(sign_robust): value not found'
call exit(1)
end

目前代码理解比较肤浅，有问题请及时指出！！！