Python實現(xiàn)克里金插值法的過程詳解

更新時間：2022年11月21日 08:51:04 作者：指尖聽?wèi)?

克里金算法提供的半變異函數(shù)模型有高斯、線形、球形、阻尼正弦和指數(shù)模型等，在對氣象要素場插值時球形模擬比較好。本文將用Python實現(xiàn)克里金插值法，感興趣的可以了解一下

一、克里金插值法介紹

克里金算法提供的半變異函數(shù)模型有高斯、線形、球形、阻尼正弦和指數(shù)模型等，在對氣象要素場插值時球形模擬比較好。既考慮了儲層參數(shù)的隨機性，有考慮了儲層參數(shù)的相關(guān)性，在滿足插值方差最小的條件下，給出最佳線性無偏插值，同時還給出了插值方差。

與傳統(tǒng)的插值方法（如最小二乘法、三角剖分法、距離加權(quán)平均法）相比，克里金法的優(yōu)勢：

1、在數(shù)據(jù)網(wǎng)格化的過程中考慮了描述對象的空間相關(guān)性質(zhì)，使插值結(jié)果更科學(xué)、更接近于實際情況；

2、能給出插值的誤差（克里金方差），使插值的可靠程度一目了然

插值方差：就是指實際參數(shù)值 zv 與估計值 zv* 兩者偏差平方的數(shù)學(xué)期望：

而插值點的 zv*，通過N個離散點獲得；

其中λ與N個離散點指的是加權(quán)系數(shù); *變差函數(shù)的理論模型*

變差函數(shù)與隨機變量的距離h存在一定的關(guān)系，這種關(guān)系可以用理論模型表示。常用的變差函數(shù)理論模型包括球狀模型、高斯模型與指數(shù)模型（還包括：具基臺值線性模型、冪函數(shù)模型、無基臺值線性模型）；

1、球狀模型公式：

2、高斯模型公式：

3、指數(shù)模型公式：

4、具基臺值線性模型：

5、冪函數(shù)模型：

式中：為冪指數(shù)；不存在基臺值。兩邊取對數(shù)得ln(γ(h))=αlnh，線性化為γ(hi)=b1X1,i

6、無基臺值線性模型：

式中：k為直線斜率；不存在基臺值和變程，當(dāng)h>0時, γ(hi)=b0+b1X1,i

普通克里格方法的基本步驟如下：

二、算法實現(xiàn)

代碼實現(xiàn):

from gma.algorithm.spmis.interpolate import *

class Kriging(IPolate):
    '''克里金法插值'''
 # 繼承 gma 的標(biāo)準(zhǔn)插值類 IPolate。本處不再做詳細說明。
    def __init__(self, Points, Values, Boundary = None, Resolution = None, 
                 InProjection = 'WGS84', 
                 VariogramModel = 'Linear',
                 VariogramParameters = None,
                 **kwargs):
        
        IPolate.__init__(self, Points, Values, Boundary, Resolution, InProjection)
        
        self.eps = eps
        
        # 初始化半變異函數(shù)及參數(shù)
        self.VariogramModel, self.VParametersList = GetVariogramParameters(VariogramModel, VariogramParameters)
        self.VariogramFUN = getattr(variogram, self.VariogramModel)
        if self.VParametersList is None:
            self.VParametersList = self._INITVariogramModel(**kwargs)
        
        # 調(diào)整輸入數(shù)據(jù)
        if self.GType == 'PROJCS':
            self.Center = (self.Points.min(axis = 0) + self.Points.max(axis = 0)) * 0.5
            self.AnisotropyScaling = AnisotropyScaling
            self.AnisotropyAngle = AnisotropyAngle
            self.DistanceMethod = cdist
        else:
            # 方便后期優(yōu)化
            self.Center = np.array([0,0])
            self.AnisotropyScaling = 1.0
            self.AnisotropyAngle = 0.0
            self.DistanceMethod = GreatCircleDistance
        
        self.AdjustPoints = AdjustAnisotropy(self.Points, self.Center, 
                                             [self.AnisotropyScaling], 
                                             [self.AnisotropyAngle])
        self.XYs = AdjustAnisotropy(self.XYs, self.Center,
                                    [self.AnisotropyScaling], 
                                    [self.AnisotropyAngle])
        
    def _INITVariogramModel(self, **kwargs):
        '''初始化參數(shù)'''
        
        if 'NLags' in kwargs:
            NLags = kwargs['NLags']
            initialize.ValueType(NLags, 'pint')
        else:
            NLags = 6
            
        if 'Weight' in kwargs:
            Weight = ToNumericArray(kwargs['Weight']).flatten().astype(bool)[0]
        else:
            Weight = False

        Lags, SEMI = INITVariogramModel(self.Points, self.Values, NLags, self.GType)
        
        # 為求解自動參數(shù)準(zhǔn)備
        if self.VariogramModel == "Linear":
            X0 = [np.ptp(SEMI) / np.ptp(Lags), np.min(SEMI)]
            BNDS = ([0.0, 0.0], [np.inf, np.max(SEMI)])
        elif self.VariogramModel == "Power":
            X0 = [np.ptp(SEMI) / np.ptp(Lags), 1.1, np.min(SEMI)]
            BNDS = ([0.0, 0.001, 0.0], [np.inf, 1.999, np.max(SEMI)])
        else:
            X0 = [np.ptp(SEMI), 0.25 * np.max(Lags),  np.min(SEMI)]
            BNDS = ([0.0, 0.0, 0.0], [10.0 * np.max(SEMI), np.max(Lags), np.max(SEMI)])
        
        # 最小二乘法求解默認(rèn)參數(shù)
        def _VariogramResiduals(Params, X, Y, Weight):
            if Weight:
                Weight = 1.0 / (1.0 + np.exp(-2.1972 / (0.1 * np.ptp(X)) * (0.7 * np.ptp(X) + np.min(X) - x))) + 1 
            else:
                Weight = 1
            return (self.VariogramFUN(X, *Params) - Y) * Weight

        RES = least_squares(_VariogramResiduals, X0, bounds = BNDS, loss = "soft_l1",
                            args = (Lags, SEMI, Weight))
                            
        return RES.x
    
    def _GetKrigingMatrix(self):
        """獲取克里金矩陣"""
            
        LDs = self.DistanceMethod(self.AdjustPoints, self.AdjustPoints)
        
        A = -self.VariogramFUN(LDs, *self.VParametersList)
        A = np.pad(A, (0, 1), constant_values = 1)
        # 填充主對角線
        np.fill_diagonal(A, 0.0)
 
        return  A
    
    def _UKExec(self, A, LDs, SearchRadius):
        """泛克里金求解"""
        Args = LDs.argsort(axis = 1)[:,:SearchRadius]
        Values = self.Values[Args.T].T
        
        # A 的逆矩陣
        AInv = inv(A)
        B = -self.VariogramFUN(LDs, *self.VParametersList)
        B[np.abs(LDs) <= self.eps] = 0.0
        B = np.pad(B, ((0,0),(0,1)), constant_values = 1)

        X = np.dot(B, AInv)
        
        B = B[np.ogrid[:len(B)], Args.T].T
        X = X[np.ogrid[:len(X)], Args.T].T
        X = X / X.sum(axis = 1, keepdims = True)

        UKResults = np.sum(X * Values, axis = 1), np.sum((X * -B), axis = 1)

        return UKResults
    
    def _OKExec(self, A, LDs, SearchRadius):
        """普通克里金求解"""
        Args = LDs.argsort(axis = 1)[:,:SearchRadius]
        LDs = LDs[np.ogrid[:len(LDs)], Args.T].T

        B = -self.VariogramFUN(LDs, *self.VParametersList)
        B[np.abs(LDs) <= self.eps] = 0.0
        B = np.pad(B, ((0,0),(0,1)), constant_values = 1)
 
        OKResults = np.zeros([2, len(LDs)])

        for i, b in enumerate(B):
            
            BSelector = Args[i] 
            ASelector = np.append(BSelector, len(self.AdjustPoints))
            a = A[ASelector[:, None], ASelector]  
            x = solve(a, b)
            
            OKResults[:, i] = x[:SearchRadius].dot(self.Values[BSelector]), -x.dot(b)

        return OKResults
    

    def Execute(self, SearchRadius = 12, KMethod = 'Ordinary'):
        '''克里金插值'''
        initialize.ValueType(SearchRadius, 'pint')
        SearchRadius = np.min([SearchRadius, len(self.AdjustPoints)])
        
        A = self._GetKrigingMatrix()

        LDs = self.DistanceMethod(self.XYs, self.AdjustPoints)
        
        if KMethod not in ['Universal', 'Ordinary']:
            raise ValueError("Undefined Kriging method. Please select 'Universal' or 'Ordinary'!")
        elif KMethod == 'Universal':
            KResults = self._UKExec(A, LDs, SearchRadius)
        else:
            KResults = self._OKExec(A, LDs, SearchRadius)
            
        NT = namedtuple('Kriging', ['Data', 'SigmaSQ', 'Transform'])
    
        return NT(KResults[0].reshape(self.YLAT, self.XLON),
                  KResults[1].reshape(self.YLAT, self.XLON), self.Transform)

三、差值應(yīng)用

示例數(shù)據(jù)可從：https://gma.luosgeo.com/ 獲取

在 gma 1.0.13.15 之后的版本可以直接引用。這里基于 1.0.13.15之后的版本引用做示例。

import gma
import pandas as pd

Data = pd.read_excel("Interpolate.xlsx")
Points = Data.loc[:, ['經(jīng)度','緯度']].values
Values = Data.loc[:, ['值']].values

# 普通克里金（球面函數(shù)模型）插值
KD = gma.smc.Interpolate.Kriging(Points, Values, Resolution = 0.05, 
                                 VariogramModel = 'Spherical', 
                                 VariogramParameters = None,
                                 KMethod = 'Ordinary',
                                 InProjection = 'EPSG:4326')

# 泛克里金類似，這里不做示例

gma.rasp.WriteRaster(r'.\gma_OKriging.tif',
                     KD.Data,
                     Projection = 'WGS84',
                     Transform = KD.Transform, 
                     DataType = 'Float32')