以前用过Manner-Kendall M-K突变检验，代码是matlab编写的。

最近在用R语言，把以前的相关工作进行了代码迁移。

网上有些Manner-Kendall M-K突变检验的R代码，但在应用中会出现问题。而且，本人能力有限，有些代码不能很好地理解。因此，现将个人整理后的R代码列出如下：

仅参考，版本归原作者。

------------------------

## 以这个网址的代码为主 http://www.r-china.net/forum.php?mod=viewthread&tid=723
## 数据和突变点输出，参考https://blog.csdn.net/weixin_39634034/article/details/102524239

    Mann_Kendall <- function(timeserial){#Mann Kendall 突变检验，传递参数
      Mann_Kendall_sub <- function(timeserial){#需要进行两次秩的分析，因此在函数中嵌套了一个函数
        r <- c()#分析的三个变量，具体含义可以参照魏凤英老师的书
        s <- c()#秩序列。
        U <- c()
       
        for(i in 2:length(timeserial))#进行大小比较，从第二个开始与以前的数据进行大小比较
        {r <- 0
         
         for(j in 1:i)
         {
           if(timeserial>timeserial[j]){r=r+1}#如果后面的数大于前面的数，则秩加1
         }
         
         
         s <- 0
         for (ii in 2:i){
           s <- s+r[ii]#秩序列。Sk是第i时刻数值大于ii时刻数值个数的累计数
         }
         
         
         U <- 0
         U <- (s- (i*(i-1)/4))/sqrt(i*(i-1)*(2*i+5)/72)
         
        }
       
        r[1] <- 0
        s[1] <- 0
        U[1] <- 0
       
        LST <- list(r = r, s = s, U = U)
       
        return (LST)
      }
      timeserial_rev <- rev(timeserial)#生成逆序列
      
      y1 <- Mann_Kendall_sub(timeserial)#计算正序列
      y2 <- Mann_Kendall_sub(timeserial_rev)#计算逆序列
      
      y2$U <- -(rev(y2$U))#转换符号与顺序
      
      LST <- list(UF=y1,UB=y2)
      return(LST)  
    }



    #这里是你要修改的地方
    
WD <- getwd() #获取工作目录
if (!is.null(WD))
  setwd(WD)
# choose.dir(WD) #选择目录
library(rio)
od <- import(choose.files(WD)) # 导入数据
# Variable <- c("Jan","Feb","Mar","Apr","May","Jun")
Variable <-od[ ,1]
#修改上面代码

#可以自己定义，或者根据数据自动生成
rows <- length(od$Year)
startyear <- as.numeric(od[1,1])
years <- od$Year
#######MK检验#######
d<-Mann_Kendall(od[,2]) #进行突变检验
yUF <- as.data.frame(d$UF[3])$U
yUB <- as.data.frame(d$UB[3])$U

#####绘图#####
plot(x=c(1:length(od$Year)),y=yUF, type="l", ylim=c(min(yUF,yUB,-1.96),max(yUF,yUB,1.96)),lwd=1, lty=5, ylab="", cex=0.5,xaxt="n",mgp=c(1,0.1,0),tck=-0.02)
points(x=c(1:length(od$Year)),y=yUB,type="l",lty=3,col=6,lwd=1)
abline(h=1.969,lty=4,lwd=0.5) # 1.969是a=0.05的显著性水平
abline(h=-1.96,lty=4,lwd=0.5)
abline(h=0,col="gray",lwd=0.5)
    
axis(side=1, at=years, tck=-0.04, hadj=0.4, labels=years,mgp=c(1,0.4,0), cex.axis=1) # add x-axis to the last figure
axis(side=1, at=1:length(od$Year), tck=-0.01, hadj=0.4, labels=NA, cex.axis=1) # add month labels to the x-axis
mtext("年份",side=1,line=1.5)

 #####输出突变年份#####
year_mk<-od[,1]
year_point<-year_mk[which((as.numeric(yUF)-(yUB))>0)[1]-1]
print(year_point)

如果有多个交叉点 突变点，怎么办？

提供思路：

1 计算yUF-yUB的值du,在突变前后，du会从正值变为负值，或从负值变为正值。

2 将du中正值赋值为1，负值负值为0。计算du1=du[i+1]-du

3 筛选du1中的非零值，这些数值对应的时间是突变点。