zhuangwei的个人博客分享 http://blog.sciencenet.cn/u/zhuangwei

博文

USE R: SCEUA优化方法进行Shuttleworth-Wallace 模型的参数优化

已有 1454 次阅读 2019-7-9 17:02 |个人分类:R|系统分类:科研笔记

sceua_SW.zip


分享一个SCEUA进行参数优化的例子,模型是Shuttleworth-Wallace 模型,其中的冠层导度子模型采用自Leuning et al. (2008) WRR. 用到了Ozflux一个站点的数据,数据为白天的日均值,优化的参数是最大气孔导度gmax,和土壤表面阻力有关的参数b3。感兴趣的老师同学可以看看。这里写的R代码是在jupyter lab运行的,把HTML中的代码,粘贴到R或者RStudio中运行结果是相同的。


Shuffled Complex Evolution (SCE-UA) Method: Parameter Optimization for SW model

This file will get the optimized parameter value, which will be used in 'conduct_simulation.ipynb' to get prediction results over the five selected savanas sites reference: https://cran.r-project.org/web/packages/rtop/rtop.pdf https://www.mathworks.com/matlabcentral/fileexchange/7671-shuffled-complex-evolution-sce-ua-method


In [1]:

# load library
library(rtop)
Warning message:
"package 'rtop' was built under R version 3.5.3"Loading required package: sp
Warning message:
"package 'sp' was built under R version 3.5.2"

In [2]:
# read input data for the modele
data <- read.csv(file='ade_flux_lai_daily.csv',head=T)
edata$date <- as.Date(edata$date) 
edata <- subset(edata, fe_qcflag<0.2 & fe >0)

In [3]:
# constant
kr    =  0.41  # Karman constance
gamma  = 0.067   #  0.067 kPa/K, Goudriaan 1977
rho    = 1.225    # air density 1.293 kg/m3
cp    = 1005     # J / (kg * K)

# Parameters:
# Input |value |Description
lwidth =  0.03  # characteristic leaf width (m)
z0gs  =  0.01  # roughness parameter of soil surface (m) # (0.01 m) (Van Bavel and Hillel,1976). zhang 2008
km     =  2.5    # wind/diffusivity extinction coefficient
rhotp  =  2      # ratio of total leaf area to projected leaf area 

# parameter in surface conductance
kq     =  0.6  # extinction coefficient of visible energy 
q50    =  30   # W/m2; the visible radiation flux at half maximum of stomatal conductance
vpd50  =  0.7  # kPa ; vapor pressure deficit at half maximum of stomatal conductance

# load site specific pararmeters
hc = 12.5; zr = 15 ; swc_min = 0.116 ; swc_max = 0.234

# parameter in soil surface resistance
b1 = 3; b2 = 4; # b3 = 600

  
In [4]:
# load parameter range for parameter needing optimization
parm_range = data.frame(parname = c('gmax', 'b3'),
             lower   = c(0.004 , 100), 
             best    = c(0.008 , 500), 
             upper   = c(0.015 , 1000))
                                          
# load reference parameter definition (best, upper, lower)
  refPars <- parm_range[, 2:4]
  rownames(refPars) = c('gmax', 'b3')
  refPars


lowerbestupper
gmax4e-038e-031.5e-02
b31e+025e+021.0e+03


In [5]:

# Objective function, summing up squared differencesOFUN = function(pars, x, yobs) {
  yvals = sw_model(x, pars)
  sum((yvals-yobs)^2)}

  

In [6]:

# load modelsource('sw_model.R')# conduct SCE opitimization
  driver = na.omit(edata)
  driver$et_pre <- sw_model(driver, refPars$best)
    
  sceuares = sceua(OFUN, pars = refPars$best, 
                   lower = refPars$lower,
               upper = refPars$upper, 
               x = driver, 
               yobs = driver$et_ob)

52 best 115 function convergence 200 parameter convergence 6103.69 

83 best 115 function convergence 200 parameter convergence 3056.697 
113 best 115 function convergence 200 parameter convergence 1834.914 
141 best 115 function convergence 200 parameter convergence 835.5961 
171 best 115 function convergence 200 parameter convergence 431.5789 
201 best 115 function convergence 0.0727 parameter convergence 203.3519 
229 best 115 function convergence 0.0757 parameter convergence 104.6551 
259 best 115 function convergence 0.0198 parameter convergence 53.47075 
289 best 115 function convergence 0.02 parameter convergence 26.28713 
321 best 115 function convergence 0.0152 parameter convergence 11.07059 
349 best 115 function convergence 0.00312 parameter convergence 5.43802 
379 best 115 function convergence 0.000146 parameter convergence 3.059501 
407 best 115 function convergence 0.000146 parameter convergence 1.520457 
435 best 115 function convergence 2.16e-05 parameter convergence 0.8048972

In [7]:

sceuares$par

0.00920174471959075

  1.  279.09394161195


In [8]:

# conduct simulation with optimized parameters
edata$et_sw <- sw_model(edata, sceuares$par)


In [9]:

# model evaluationd_lm = with(edata,lm(et_sw ~ et_ob))
  mevals = c(r2 = summary(d_lm)$r.squared,   # R2
             itcp = d_lm$coefficients[1],      # intercept
             slp  = d_lm$coefficients[2],      # slope
             num  = dim(edata)[1],               # N
             mae  = sum(abs(edata$et_ob - edata$et_sw))/dim(edata)[1],           # MAE
             rmse = (sum((edata$et_ob - edata$et_sw)^2)/dim(edata)[1])^0.5)      # RMSE


In [10]:

mevals


  • r2

  • 0.904015521163113

  • itcp.(Intercept)

  • 0.436897410134342

  • slp.et_ob

  • 0.862991255900808

  • num

  • 480

  • mae

  • 0.374697611822356

  • rmse

  • 0.488530739708623

In [11]:

with(edata, plot(et_ob, et_sw, 
                 xlim=c(0,7), ylim=c(0,7),
                 xlab='observed ET mm/day',
                 ylab='predicted ET mm/day'))





http://blog.sciencenet.cn/blog-526092-1188776.html

上一篇:气孔导度单位换算 更正
下一篇:Century 模型 R代码分享

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-11-24 09:43

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部