- 积分
- 3638
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-10-21
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 15195775117 于 2020-12-24 09:39 编辑
一、缘起
克里金法是我工作之后才知道的插值法,其插值效果极其优美,是我平时用的最多的插值法。
二、特性
克里金法(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 |
|