Cokriging代理模型理论相关推导

Cokriging模型是20世纪70年代发展起来的一种地质统计学插值模型。在地质统计学领域，为了提高对某个抽样比较困难的量的预测精度，提出了采用更容易抽样的量进行辅助预测的Kriging模型，称为Cokriging模型，本文对Cokriging模型相关理论进行了推导。

2.1 代理模型问题的基本描述

对于一个有$m$个设计变量的优化问题，在设计空间中同时进行高可信度分析和低可信度抽样，以建立所谓变可信度模型，变可信代理模型在达到相同近似精度的条件下，可以显著提高建立代理模型的效率。设高、低可信度分析程序的抽样位置分别为：
$$\{\begin{array}{l}{\mathbf{S}}_1=[{\mathbf{x}}_1^{(1)}{\mathbf{x}}_1^{(2)}\dots {\mathbf{x}}_1^{(n_1)}{]}^T\in {\mathbb{R}}^{n_1\times m}\\ {\mathbf{S}}_2=[{\mathbf{x}}_2^{(1)}{\mathbf{x}}_2^{(2)}\dots {\mathbf{x}}_2^{(n_2)}{]}^T\in {\mathbb{R}}^{n_2\times m}\end{array}\text{\hspace{1em}(1)}$$
上式中，下标 “1” 和 “2” 分别代表高、低可信度，例如 $n_1$ 和 $n_2$ 分别代表高、低可信度样本点数（可合理假设 $n_2\gg n_1$ ）。相应的目标函数或约束函数响应值为：
$$\{\begin{array}{l}{\mathbf{y}}_1=[y_1^{(1)}{y}_1^{(2)}\dots y_1^{(n_1)}{]}^T\in {\mathbb{R}}^{n_1}\\ {\mathbf{y}}_2=[y_2^{(1)}{y}_2^{(2)}\dots y_2^{(n_2)}{]}^T\in {\mathbb{R}}^{n_2}\end{array}\text{\hspace{1em}(2)}$$
结合高可信度数据集 $({\mathbf{S}}_1,{\mathbf{y}}_1)$ 和低可信度数据 $({\mathbf{S}}_1,{\mathbf{y}}_2)$ 可以建立 Cokriging 模型。

2.2 Cokriging 模型及其预估值

Cokriging 模型预估值定义为：
$${\hat{y}}_1(\mathbf{x})={\mathbf{\lambda }}^T{\mathbf{y}}_s={\mathbf{\lambda }}_1^T{\mathbf{y}}_1+{\mathbf{\lambda }}_2^T{\mathbf{y}}_2\text{\hspace{1em}(3)}$$
上式中：${\mathbf{\lambda }}_1$ 、${\mathbf{\lambda }}_2$ 分别为对高低可信度响应值的加权系数，只要给出加权系数 ${\mathbf{\lambda }}_1=[{\lambda }_1^{(1)}{\lambda }_1^{(2)}\dots {\lambda }_1^{n_1}]$ 和 ${\mathbf{\lambda }}_2=[{\lambda }_2^{(1)}{\lambda }_2^{(2)}\dots {\lambda }_2^{(n_2)}]$ ，就可以得到设计空间中任意点的响应值。为了计算这两个加权系数，引入统计学假设，假设存在分别与 $y_1$ 、$y_2$ 对应的两个高斯静态随机过程：
$$\{\begin{array}{l}{Y}_1(\mathbf{x})={\beta }_1+Z_1(\mathbf{x})\\ Y_2(\mathbf{x})={\beta }_2+Z_2(\mathbf{x})\end{array}\text{\hspace{1em}(4)}$$
上式中，${\beta }_1$ 、${\beta }_2$ 均为未知常数，也称为全局趋势模型，它们分别代表 $Y_1(\mathbf{x})$ 和 $Y_2(\mathbf{x})$ 的数学期望值；$Z_1(\cdot )$ 、$Z_2(\cdot )$ 是均值为0，方差分别为 ${\sigma }_1^2$ 和 ${\sigma }_2^2$ 的静态随机过程。在设计空间不同位置处，随机变量之间的协方差和交叉协方差定义为：
$$\{\begin{array}{l}Cov(Z({\mathbf{x}}_1^{(i)}),Z({\mathbf{x}}_1^{(j)}))={\sigma }_1^2{R}^{(11)}({\mathbf{x}}_1^{(i)},{\mathbf{x}}_1^{(j)})\\ Cov(Z({\mathbf{x}}_2^{(i)}),Z({\mathbf{x}}_2^{(j)}))={\sigma }_2^2{R}^{(22)}({\mathbf{x}}_2^{(i)},{\mathbf{x}}_2^{(j)})\\ Cov(Z({\mathbf{x}}_1^{(i)}),Z({\mathbf{x}}_2^{(j)}))={\sigma }_1{\sigma }_2{R}^{(12)}({\mathbf{x}}_1^{(i)},{\mathbf{x}}_2^{(j)})\end{array}\text{\hspace{1em}(5)}$$
上式中，${\sigma }_1^2$ 和 ${\sigma }_2^2$ 分别为随机过程 $Y_1(\mathbf{x})$ 和 $Y_2(\mathbf{x})$ 的过程方差。

基于上述假设，Cokriging 模型寻找最优加权系数 ${\mathbf{\lambda }}_1$ 和 ${\mathbf{\lambda }}_2$ 使得均方差：
$$MSE({\hat{y}}_1(\mathbf{x}))=E[({\mathbf{\lambda }}_1^T{\mathbf{Y}}_1+{\mathbf{\lambda }}_2^T{\mathbf{Y}}_2-Y_1(\mathbf{x}){)}^2]\text{\hspace{1em}(6)}$$
最小，并引入无偏估计：
$$E[{\mathbf{\lambda }}_1^T{\mathbf{Y}}_1+{\mathbf{\lambda }}_2^T{\mathbf{Y}}_2]=E[Y_1]\text{\hspace{1em}(7)}$$
无偏估计化简得：
$${\beta }_1(\sum _{i=1}^{n_1}{\lambda }_1^{(i)}-1)+{\beta }_2\sum _{j=1}^{n_2}{\lambda }_2^{(j)}=0\text{\hspace{1em}(8)}$$
为了让无偏估计条件不依赖于 ${\beta }_1$ 和 ${\beta }_2$ 的取值，取如下更强一些的无偏估计条件：
$$\sum _{i=1}^{n_1}{\lambda }_1^{(i)}-1=0,\sum _{i=1}^{n_2}{\lambda }_2^{(i)}=0\text{\hspace{1em}(9)}$$

采用拉格朗日乘数法，令：
$$L({\mathbf{\lambda }}_1,{\mathbf{\lambda }}_2,{\mu }_1,{\mu }_2)=E[({\mathbf{\lambda }}_1^T{\mathbf{Y}}_1+{\mathbf{\lambda }}_2^T{\mathbf{Y}}_2-Y_1(\mathbf{x}){)}^2]+{\mu }_1(\sum _{i=1}^{n_1}{\lambda }_1^{(i)}-1)+{\mu }_2\sum _{j=1}^{n_2}{\lambda }_2^{(j)}\text{\hspace{1em}(10)}$$
令：
$$\{\begin{array}{l}\frac{\partial L({\mathbf{\lambda }}_1,{\mathbf{\lambda }}_2,{\mu }_1,{\mu }_2)}{\partial {\lambda }_1^{(i)}}=E[2({\mathbf{\lambda }}_1^T{\mathbf{Y}}_1+{\mathbf{\lambda }}_2^T{\mathbf{Y}}_2-Y_1(\mathbf{x}))\ast Y_1({\mathbf{x}}_1^{(i)})]+{\mu }_1=0\\ \frac{\partial L({\mathbf{\lambda }}_1,{\mathbf{\lambda }}_2,{\mu }_1,{\mu }_2)}{\partial {\lambda }_2^{(i)}}=E[2({\mathbf{\lambda }}_1^T{\mathbf{Y}}_1+{\mathbf{\lambda }}_2^T{\mathbf{Y}}_2-Y_1(\mathbf{x}))\ast Y_2({\mathbf{x}}_2^{(i)})]+{\mu }_2=0\\ \frac{\partial L({\mathbf{\lambda }}_1,{\mathbf{\lambda }}_2,{\mu }_1,{\mu }_2)}{\partial {\mu }_1}={\sum }_{i=1}^{n_1}{\lambda }_1^{(i)}-1=0\\ \frac{\partial L({\mathbf{\lambda }}_1,{\mathbf{\lambda }}_2,{\mu }_1,{\mu }_2)}{\partial {\mu }_2}={\sum }_{i=1}^{n_2}{\lambda }_2^{(i)}=0\end{array}\text{\hspace{1em}(11)}$$
以式 (9) 中第一行式子为例演示展开过程：
$$\begin{array}{rl}\frac{\partial L({\mathbf{\lambda }}_1,{\mathbf{\lambda }}_2,{\mu }_1,{\mu }_2)}{\partial {\lambda }_1^{(i)}}& =2E[(\sum _{j=1}^{n_1}{\lambda }_1^{(j)}{Y}_1({\mathbf{x}}_1^{(j)}))\cdot Y_1({\mathbf{x}}_1^{(i)})]+2E[(\sum _{j=1}^{n_2}{\lambda }_2^{(j)}{Y}_2({\mathbf{x}}_2^{(j)}))\cdot Y_1({\mathbf{x}}_1^{(i)})]-2E[Y_1({\mathbf{x}}_1)\cdot Y_1({\mathbf{x}}_1^{(i)})]+{\mu }_1\\ & =2\sum _{j=1}^{n_1}{\lambda }_1^{(j)}({\beta }_1^2+{\sigma }_1^2{R}^{11}({\mathbf{x}}_1^{(i)},{\mathbf{x}}_1^{(j)}))+2\sum _{j=1}^{n_2}{\lambda }_2^{(j)}({\beta }_1{\beta }_2+{\sigma }_1{\sigma }_2{R}^{12}({\mathbf{x}}_1^{(i)},{\mathbf{x}}_2^{(j)}))-2({\beta }_1^2+{\sigma }_1^{2}R({\mathbf{x}}_1^{(i)},{\mathbf{x}}_1))+{\mu }_1\\ & =2(\sum _{j=1}^{n_1}{\lambda }_1^{(j)}-1){\beta }_1^2+2{\sigma }_1^2(\sum _{j=1}^{n_1}{\lambda }_1^{(j)}R({\mathbf{x}}_1^{(i)},{\mathbf{x}}_1^{(j)}))+2{\beta }_1{\beta }_2\sum _{j=1}^{n_2}{\lambda }_2^{(j)}+2{\sigma }_1{\sigma }_2(\sum _{j=1}^{n_2}{\lambda }_2^{(j)}{R}^{12}({\mathbf{x}}_1^{(i)},{\mathbf{x}}_2^{(j)}))-2{\sigma }_1^2{R}^{11}({\mathbf{x}}_1,{\mathbf{x}}_1^{(i)})+{\mu }_1\\ & =-2{\beta }_1(\sum _{j=1}^{n_2}{\lambda }_2^{(j)}{\beta }_2)+2{\sigma }_1^2(\sum _{j=1}^{n_1}{\lambda }_1^{(j)}R({\mathbf{x}}_1^{(i)},{\mathbf{x}}_1^{(j)}))+2{\beta }_1{\beta }_2\sum _{j=1}^{n_2}{\lambda }_2^{(j)}+2{\sigma }_1{\sigma }_2(\sum _{j=1}^{n_2}{\lambda }_2^{(j)}{R}^{12}({\mathbf{x}}_1^{(i)},{\mathbf{x}}_2^{(j)}))-2{\sigma }_1^2{R}^{11}({\mathbf{x}}_1,{\mathbf{x}}_1^{(i)})+{\mu }_1\\ & =2{\sigma }_1^2(\sum _{j=1}^{n_1}{\lambda }_1^{(j)}R({\mathbf{x}}_1^{(i)},{\mathbf{x}}_1^{(j)})+2{\sigma }_1{\sigma }_2(\sum _{j=1}^{n_2}{\lambda }_2^{(j)}{R}^{12}({\mathbf{x}}_1^{(i)},{\mathbf{x}}_2^{(j)}))-2{\sigma }_1^2{R}^{11}({\mathbf{x}}_1,{\mathbf{x}}_1^{(i)})+{\mu }_1\end{array}\text{\hspace{1em}(12)}$$
同样的过程可以展开式(9)中第二行式子得：
$$\frac{\partial L({\mathbf{\lambda }}_1,{\mathbf{\lambda }}_2,{\mu }_1,{\mu }_2)}{\partial {\lambda }_2^{(i)}}=2{\sigma }_1{\sigma }_2(\sum _{j=1}^{n_1}{\lambda }_1^{(j)}{R}^{12}({\mathbf{x}}_1^{(j)},{\mathbf{x}}_2^{(i)}))+2{\sigma }_2^2(\sum _{j=1}^{n_2}{\lambda }_2^{(j)}{R}^{22}({\mathbf{x}}_2^{(i)},{\mathbf{x}}_2^{(j)}))-2{\sigma }_1{\sigma }_2{R}^{12}({\mathbf{x}}_1,{\mathbf{x}}_2^{(i)})+{\mu }_2\text{\hspace{1em}(13)}$$
则式 (9) 简化如下：
$$\{\begin{array}{l}2{\sigma }_1^2({\sum }_{j=1}^{n_1}{\lambda }_1^{(j)}R({\mathbf{x}}_1^{(i)},{\mathbf{x}}_1^{(j)})+2{\sigma }_1{\sigma }_2({\sum }_{j=1}^{n_2}{\lambda }_2^{(j)}{R}^{12}({\mathbf{x}}_1^{(i)},{\mathbf{x}}_2^{(j)}))-2{\sigma }_1^2{R}^{11}({\mathbf{x}}_1,{\mathbf{x}}_1^{(i)})+{\mu }_1=0\\ 2{\sigma }_1{\sigma }_2({\sum }_{j=1}^{n_1}{\lambda }_1^{(j)}{R}^{12}({\mathbf{x}}_1^{(j)},{\mathbf{x}}_2^{(i)}))+2{\sigma }_2^2({\sum }_{j=1}^{n_2}{\lambda }_2^{(j)}{R}^{22}({\mathbf{x}}_2^{(i)},{\mathbf{x}}_2^{(j)}))-2{\sigma }_1{\sigma }_2{R}^{12}({\mathbf{x}}_1,{\mathbf{x}}_2^{(i)})+{\mu }_2=0\\ {\sum }_{i=1}^{n_1}{\lambda }_1^{(i)}-1=0\\ {\sum }_{i=1}^{n_2}{\lambda }_2^{(i)}=0\end{array}\text{\hspace{1em}(14)}$$
改写为矩阵形式有：
$$\left[\begin{array}{cccc}{\mathbf{C}}^{(11)}& {\mathbf{C}}^{(12)}& {\mathbf{F}}_1& 0\\ {\mathbf{C}}^{(21)}& {\mathbf{C}}^{(22)}& 0& {\mathbf{F}}_2\\ {\mathbf{F}}_1^T& 0^T& 0& 0\\ 0^T& {\mathbf{F}}_2^T& 0& 0\end{array}\right]\left[\begin{array}{c}{\mathbf{\lambda }}_1\\ {\mathbf{\lambda }}_2\\ {\mu }_1/2\\ {\mu }_2/2\end{array}\right]=\left[\begin{array}{c}{\mathbf{c}}_1\\ {\mathbf{c}}_2\\ 1\\ 0\end{array}\right]\text{\hspace{1em}(13)}$$
上式中：
$$\begin{array}{rl}& {\mathbf{C}}^{(11)}=({\sigma }_1^2{R}^{11}({\mathbf{x}}_1^{(i)},{\mathbf{x}}_1^{(j)}){)}_{i,j}\in {\mathbb{R}}^{n_1\times n_1}\\ & {\mathbf{C}}^{(12)}=({\sigma }_1{\sigma }_2{R}^{12}({\mathbf{x}}_1^{(i)},{\mathbf{x}}_2^{(j)}){)}_{i,j}=({\mathbf{C}}^{(21)}{)}^T\in {\mathbb{R}}^{n_1\times n_2}\\ & {\mathbf{C}}^{(22)}=({\sigma }_2^2{R}^{22}({\mathbf{x}}_2^{(i)},{\mathbf{x}}_2^{(j)}){)}_{i,j}\in {\mathbb{R}}^{n_2\times n_2}\\ & {\mathbf{c}}_1=({\sigma }_1^2{R}^{11}({\mathbf{x}}_1^{(i)},{\mathbf{x}}_1){)}_i\in {\mathbb{R}}^{n_1}\\ & {\mathbf{c}}_2=({\sigma }_1{\sigma }_2{R}^{12}({\mathbf{x}}_2^{(i)},{\mathbf{x}}_1){)}_i\in {\mathbb{R}}^{n_2}\\ & {\mathbf{F}}_1=[11\dots 1{]}_i^T,i\in [1,n_1]\\ & {\mathbf{F}}_2=[11\dots 1{]}_i^T,i\in [1,n_2]\end{array}\text{\hspace{1em}(14)}$$
为了在让 ${\sigma }_1$ 和 ${\sigma }_2$ 在式（13）中消去，设：
$$\begin{array}{rl}& {\tilde{\mathbf{\lambda }}}_1={\mathbf{\lambda }}_1,{\tilde{\mathbf{\lambda }}}_2=\frac{{\sigma }_2}{{\sigma }_1}{\mathbf{\lambda }}_2\\ & {\tilde{\mu }}_1={\mu }_1/(2{\sigma }_1^2),{\tilde{\mu }}_2={\mu }_2/(2{\sigma }_1{\sigma }_2)\end{array}\text{\hspace{1em}(15)}$$
则式（13）化为：
$$\left[\begin{array}{cccc}{\mathbf{R}}^{(11)}& {\mathbf{R}}^{(12)}& {\mathbf{F}}_1& 0\\ {\mathbf{R}}^{(21)}& {\mathbf{R}}^{(22)}& 0& {\mathbf{F}}_2\\ {\mathbf{F}}_1^T& 0^T& 0& 0\\ 0^T& {\mathbf{F}}_2^T& 0& 0\end{array}\right]\left[\begin{array}{c}{\tilde{\mathbf{\lambda }}}_1\\ {\tilde{\mathbf{\lambda }}}_2\\ {\tilde{\mu }}_1\\ {\tilde{\mu }}_2\end{array}\right]=\left[\begin{array}{c}{\mathbf{r}}_1(\mathbf{x})\\ {\mathbf{r}}_2(\mathbf{x})\\ 1\\ 0\end{array}\right]\text{\hspace{1em}(16)}$$
上式中：
$$\begin{array}{rl}& {\mathbf{R}}^{(11)}=(R^{11}({\mathbf{x}}_1^{(i)},{\mathbf{x}}_1^{(j)}){)}_{i,j}\in {\mathbb{R}}^{n_1\times n_1}\\ & {\mathbf{R}}^{(12)}=(R^{12}({\mathbf{x}}_1^{(i)},{\mathbf{x}}_2^{(j)}){)}_{i,j}={\mathbf{R}}^{(21)}{)}^T\in {\mathbb{R}}^{n_1\times n_2}\\ & {\mathbf{R}}^{(22)}=(R^{22}({\mathbf{x}}_2^{(i)},{\mathbf{x}}_2^{(j)}){)}_{i,j}\in {\mathbb{R}}^{n_2\times n_2}\\ & {\mathbf{r}}_1(\mathbf{x})=(R^{11}({\mathbf{x}}_1,{\mathbf{x}}_1^{(i)}){)}_i\in {\mathbb{R}}^{n_1}\\ & {\mathbf{r}}_2(\mathbf{x})=(R^{12}({\mathbf{x}}_1,{\mathbf{x}}_2^{(i)}){)}_i\in {\mathbb{R}}^{n_2}\end{array}\text{\hspace{1em}(17)}$$
则 Cokriging 模型表达式：
$$\hat{y}(\mathbf{x})={\tilde{\mathbf{\lambda }}}_1^T{\mathbf{Y}}_{1,s}+\frac{{\sigma }_1}{{\sigma }_2}{\tilde{\mathbf{\lambda }}}_2{\mathbf{Y}}_{2,s}=\left[\begin{array}{cc}{\tilde{\mathbf{\lambda }}}_1& {\tilde{\mathbf{\lambda }}}_2\end{array}\right]\left[\begin{array}{c}{\mathbf{y}}_{1,s}\\ \frac{{\sigma }_1}{{\sigma }_2}{\mathbf{y}}_{2,s}\end{array}\right]\text{\hspace{1em}(18)}$$
结合式（16），则其矩阵形式为：
$${\hat{y}}_1(\mathbf{x})={\left[\begin{array}{c}{\mathbf{r}}_1(\mathbf{x})\\ {\mathbf{r}}_2(\mathbf{x})\\ 1\\ 0\end{array}\right]}^T{\left[\begin{array}{cccc}{\mathbf{R}}^{(11)}& {\mathbf{R}}^{(12)}& {\mathbf{F}}_1& 0\\ {\mathbf{R}}^{(21)}& {\mathbf{R}}^{(22)}& 0& {\mathbf{F}}_2\\ {\mathbf{F}}_1^T& 0^T& 0& 0\\ 0^T& {\mathbf{F}}_2^T& 0& 0\end{array}\right]}^{-1}\left[\begin{array}{c}{\mathbf{y}}_{1,s}\\ \frac{{\sigma }_1}{{\sigma }_2}{\mathbf{y}}_{2,s}\\ 0\\ 0\end{array}\right]\text{\hspace{1em}(19)}$$
参考 Kriging 模型推导过程，可得 Cokriging 模型表达式如下：
$${\hat{y}}_1(\mathbf{x})=\mathbf{\phi }\tilde{\mathbf{\beta }}+\mathbf{r}(\mathbf{x}){\mathbf{R}}^{-1}({\tilde{\mathbf{y}}}_s-\mathbf{F}\tilde{\mathbf{\beta }})\text{\hspace{1em}(20)}$$
上式中：
$$\begin{array}{rl}& \mathbf{\phi }=\left[\begin{array}{c}1\\ 0\end{array}\right],\tilde{\mathbf{\beta }}=\left[\begin{array}{c}{\tilde{\beta }}_1\\ {\tilde{\beta }}_2\end{array}\right]=({\mathbf{F}}^T{\mathbf{R}}^{-1}\mathbf{F}{)}^{-1}{\mathbf{F}}^T{\mathbf{R}}^{-1}{\tilde{\mathbf{y}}}_s,\\ & \mathbf{r}=\left[\begin{array}{c}{\mathbf{r}}_1(\mathbf{x})\\ {\mathbf{r}}_2(\mathbf{x})\end{array}\right],\mathbf{R}=\left[\begin{array}{cc}{\mathbf{R}}^{(11)}& {\mathbf{R}}^{(12)}\\ {\mathbf{R}}^{(21)}& {\mathbf{R}}^{(22)}\end{array}\right],\\ & {\tilde{\mathbf{y}}}_s=\left[\begin{array}{c}{\mathbf{y}}_1\\ \frac{{\sigma }_1}{{\sigma }_2}{\mathbf{y}}_2\end{array}\right],\tilde{\mathbf{F}}=\left[\begin{array}{cc}{\mathbf{F}}_1& 0\\ 0& {\mathbf{F}}_2\end{array}\right]\end{array}\text{\hspace{1em}(21)}$$
在式（20）中，令 $v_{cok}={\mathbf{R}}^{-1}({\tilde{\mathbf{y}}}_s-\mathbf{F}\tilde{\mathbf{\beta }})$ ，那么 $\mathbf{\phi }\tilde{\beta }$ 和 $v_{cok}$ 只与已知样本点有关，已未知点无关，可以在模型训练结束后一次性计算并存储。之后，预测任意 $\mathbf{x}$ 处的响应值只需要计算 $\mathbf{r}(\mathbf{x})$ 与 ${\mathbf{V}}_{krig}$ 间的点乘。

同理得到均方差为：
$$MSE[{\hat{y}}_1(\mathbf{x})]={\sigma }_1^2[1.0-{\mathbf{r}}^T{\mathbf{R}}^{-1}\mathbf{r}+({\mathbf{r}}^T{\mathbf{R}}^{-1}\tilde{\mathbf{F}}-\mathbf{\phi })({\tilde{\mathbf{F}}}^T{\mathbf{R}}^{-1}\mathbf{F}{)}^{-1}({\mathbf{r}}^T{\mathbf{R}}^{-1}\tilde{\mathbf{F}}-\mathbf{\phi }{)}^T]\text{\hspace{1em}(22)}$$

2.3 相关函数

在式 CoKriging模型表达式中，相关矩阵 $\mathbf{R}$ 和相关矢量 $\mathbf r $ 构造涉及相关函数的选择和计算，这一过程与 Kriging 模型建立过程是类似的。下表为常用的相关函数：

相关函数	$\mathbf{R}(\mathbf{x}-{\mathbf{x}}^{\prime })$
Exponential function	R ( x , x ′ ) = e x p { − ∑ i = 1 m θ i / ∣ x i − x i ′ ∣ } \mathbf R(\mathbf x , \mathbf x')=exp\{-\sum_{i=1}^m\theta_i/ \lvert {x_i -x_i'} \rvert \} R(x,x′)=exp{−∑i=1m​θi​/∣xi​−xi′​∣}
Power function	R ( x , x ′ ) = e x p { − ∑ i = 1 m θ i / ( x i − x i ′ ) 2 } \mathbf R(\mathbf x , \mathbf x')=exp\{-\sum_{i=1}^m\theta_i/(x_i-x_i')^2\} R(x,x′)=exp{−∑i=1m​θi​/(xi​−xi′​)2}
Gaussian function	R ( x , x ′ ) = e x p { − ∑ i = 1 m θ i ( x i − x i ′ ) p i } \mathbf R(\mathbf x , \mathbf x')=exp\{-\sum_{i=1}^m\theta_i(x_i-x_i')^{p_i}\} R(x,x′)=exp{−∑i=1m​θi​(xi​−xi′​)pi​}
Linear function	R ( x , x ′ ) = m a x { 0 , 1 − ∑ i = 1 m θ i / ∣ x i − x i ′ ∣ } \mathbf R(\mathbf x , \mathbf x')=max\{0,1-\sum_{i=1}^m\theta_i/ \lvert x_i-x_i' \rvert \} R(x,x′)=max{0,1−∑i=1m​θi​/∣xi​−xi′​∣}
Ball function	R ( x , x ′ ) = 1 − 1.5 ζ + 0.5 ζ 3 , ζ = m i n { 1 , ∑ i = 1 m θ i ∣ x i − x i ′ ∣ } \mathbf R (\mathbf x,\mathbf x')=1-1.5\zeta+0.5 \zeta ^3, \\ \zeta=min\{ 1,\sum_{i=1}^m\theta_i \lvert x_i-x_i' \rvert \} R(x,x′)=1−1.5ζ+0.5ζ3,ζ=min{1,∑i=1m​θi​∣xi​−xi′​∣}
Spline function	$\mathbf{R}(\mathbf{x},{\mathbf{x}}^{\prime })=\{\begin{array}{l}1-15({\sum }_{i=1}^m{\theta }_i|x_i-x_i^{\prime }|{)}^2+30({\sum }_{i=1}^m{\theta }_i|x_i-x_i^{\prime }|{)}^{3}0\le {\sum }_{i=1}^m{\theta }_i|x_i-x_i^{\prime }|\le 0.2\\ 1.25(1-({\sum }_{i=1}^m{\theta }_i|x_i-x_i^{\prime }|{)}^3)0.2\le {\sum }_{i=1}^m{\theta }_i|x_i-x_i^{\prime }|\le 1\\ 0{\sum }_{i=1}^m{\theta }_i|x_i-x_i^{\prime }|\ge 1\end{array}$

2.4 模型参数训练

采用“最大似然估计”对 Cokriging 模型超参数进行训练，由前述引入的高斯静态随机过程的假设，在任意一点 $\mathbf{x}$ 处的响应值 $\hat{y}(\mathbf{x})$ 服从正态分布 $N({\beta }_0,{\sigma }^2)$ ，则每个样本点处的概率密度函数为：

$$P[Y(\mathbf{x})]=\frac{1}{\sqrt{2\pi }\sigma }exp\left(-\frac{(y^{(i)}-{\beta }_0{)}^2}{2{\sigma }^2}\right),i=1,\dots ,n_1+n_2\text{\hspace{1em}(23)}$$
则所有样本点的联合分布密度（即似然函数）为：
$$L({\tilde{\mathbf{\beta }}}_0,{\sigma }_1/{\sigma }_2,{\mathbf{\theta }}^{(11)},{\mathbf{\theta }}^{(12)},{\mathbf{\theta }}^{(22)})=\frac{1}{\sqrt{(2\pi {\sigma }_1^2{)}^{(n_1+n_2)}|\mathbf{R}|}}\cdot exp(-\frac{1}{2}\frac{({\tilde{\mathbf{y}}}_s-\mathbf{F}\tilde{\mathbf{\beta }}{)}^T{\mathbf{R}}^{-1}({\tilde{\mathbf{y}}}_s-\mathbf{F}\tilde{\mathbf{\beta }})}{{\sigma }^2})\text{\hspace{1em}(24)}$$

与Kriging模型推导过程相似，让似然函数最大，可以得到：
$$\begin{array}{rl}& \tilde{\mathbf{\beta }}=({\mathbf{F}}^T{\mathbf{R}}^{-1}\mathbf{F}{)}^{-1}{\mathbf{F}}^T{\mathbf{R}}^{-1}{\tilde{\mathbf{y}}}_s\\ & \frac{{\sigma }_1}{{\sigma }_2}={\left(\begin{array}{c}{\left[\begin{array}{c}0\\ {\mathbf{y}}_2\end{array}\right]}^T{\mathbf{R}}^{-1}\left[\begin{array}{c}0\\ {\mathbf{y}}_2\end{array}\right]\end{array}\right)}^{-1}{\left[\begin{array}{c}0\\ {\mathbf{y}}_2\end{array}\right]}^T{\mathbf{R}}^{-1}\left[\begin{array}{c}-({\mathbf{y}}_1-{\mathbf{F}}_1{\tilde{\beta }}_1)\\ -{\mathbf{F}}_2{\tilde{\beta }}_2\end{array}\right]\\ & {\sigma }_1^2=\frac{({\tilde{\mathbf{y}}}_s-\mathbf{F}\tilde{\mathbf{\beta }}{)}^T{\mathbf{R}}^{-1}({\tilde{\mathbf{y}}}_s-\mathbf{F}\tilde{\mathbf{\beta }}{)}^T}{n_1+n_2}\end{array}\text{\hspace{1em}(25)}$$
将上式代入式（22）并对等式两边取对数有：
$$ln[L({\mathbf{\theta }}^{(11)},{\mathbf{\theta }}^{(12)},{\mathbf{\theta }}^{22})]=-\frac{1}{2}[(n_1+n_2)\cdot ln({\sigma }_1^2)+ln(|\mathbf{R}|)]\text{\hspace{1em}(26)}$$
无法解析地求出 ${\mathbf{\theta }}^{(11)}$ 、${\mathbf{\theta }}^{(12)}$和 ${\mathbf{\theta }}^{(22)}$ 的最优值，需要采用数值优化算法寻优，如梯度优化算法等。

参考文献

[1] Han Z , Zimmerman R , GoRtz S . Alternative Cokriging Method for Variable-Fidelity Surrogate Modeling[J]. Aiaa Journal, 2012, 50(5):1205-1210.