IDL之通用版克里金Kriging插值程序

15195775117

一、缘起

克里金法是我工作之后才知道的插值法，其插值效果极其优美，是我平时用的最多的插值法。

二、特性

克里金法（Kriging）是依据协方差函数对随机过程/随机场进行空间建模和预测（插值）的回归算法。在特定的随机过程，例如固有平稳过程中，克里金法能够给出最优线性无偏估计，因此在地统计学中也被称为空间最优无偏估计器，被应用于地理科学、环境科学、大气科学等领域。
克里金法是典型的地统计学算法，以鄙人浅显的经验，它的插值结果与原始散点的值并不是绝对对应，而是整体一致但可能有局部差别。

三、IDL中的克里金

IDL中有二维克里金插值函数Krig2D，但并不是全自动，克里金插值法有几种常用插值模型，需要人为提供模型参数。
我也买了2本专业书来研究克里金，但是推导公式太麻烦了，于是我习惯于使用help中提供的参数。
没有了配置模型参数的烦恼，克里金法的参数只有经度、纬度、值3者，根据经验，一个位置不能有多个值，于是我写了去除重复位置的函数，对于同一个位置的值，做取平均处理。

四、去除重复位置的函数

;通用Krig的辅助函数
;剔除重复[x,y]位置的程序,lat,lon,ext是向量
;对于同一的位置的值，取平均
function lat_lon_ext,lat,lon,ext
  ;如果向量大小不一致就报错：
  n1=n_elements(ext)
  n2=n_elements(lon)
  n3=n_elements(lat)
  if(n1 ne n2 or n2 ne n3)then begin
    warning=dialog_message('数量不一致!')
    stop
  endif

  n=n1;数量一致用n表示
  book=[-123,-123,-123,-123]
  for i=0,n-1 do begin
    p=where(lon eq lon and lat eq lat,count)
    book=[[book],[lat,lon,count,mean(ext[p])]]
  endfor
  book=book[*,1:-1]
  sizebook=size(book)
  book2=[-123,-123,-123,-123]
  for i=0,sizebook[2]-1 do begin
    p=where(book2[0,*] eq book[0,i] and book2[1,*] eq book[1,i],count)
    if(count eq 0)then begin
      book2=[[book2],[book[*,i]]]
    endif
    if count eq 1 then continue
    if(count gt 1 or count lt 0)then begin
      warning=DIALOG_MESSAGE('出错!')
      stop
    endif
  endfor
  book2=book2[*,1:-1]
  lon2=transpose(book2[1,*])
  lat2=transpose(book2[0,*])
  value2=transpose(book2[3,*])
  re={lon:lon2,lat:lat2,value:value2}
  return,re
end

五、通用克里金函数

;返回值结构体中增设了fail用于判断插值是否失败
function GeneralKrig,lon1,lat1,value1
  ;插值成败指示变量fail：
  fail=0
  ;去除重复地点值：
  k=lat_lon_ext(lat1,lon1,value1)
  ;krig插值，使用了IDL的krig2d的help中的模型参数：
  value2=KRIG2D(k.value,k.lon,k.lat, EXPONENTIAL=[0.5,0.2,1], $
    NX=200, NY=300, XOUT=xout, YOUT=yout);NX和NY表示返回值数组的大小
  ;如果value2中有nan值，则插值失败：
  p=where(finite(value2) eq 0,count)
  if(count gt 0)then begin
    warning=dialog_message('克里金插值失败!')
    fail=1
  endif
  re={lon:xout,lat:yout,value:value2,fail:fail}
  return,re
end

六、调用方法

去除重复点的lat_lon_ext函数不用管，它是给GeneralKrig使用的，
re=GeneralKrig(经度序列,纬度序列,值序列)
对于返回值先看re.fail是否是0,0是成功,1是插值失败
插值后的经度数组是re.lon
插值后的纬度数组是re.lat
插值后的值数组是re.value

赵家清欢

您真是好厉害！这也是我要学的！

15195775117

赵家清欢 发表于 2021-10-9 16:53
您真是好厉害！这也是我要学的！

过誉了，共同进步！

richie727

共同进步！

金今

EXPONENTIAL如何调整，有没有推荐呢？

15195775117

金今 发表于 2022-7-17 08:36
EXPONENTIAL如何调整，有没有推荐呢？

这里的模型参数我也不太懂