公司网站模板源码,win7和WordPress的关系,青岛胶州网站建设,网站建设制作确认单有一篇大神的文章写得非常的具体#xff0c;
https://xg1990.com/blog/archives/222
下面写下一些学习笔记#xff1a;
1、关于克里金插值的基本原理
克里金插值来源于地理学#xff0c;它的前提是地理学第一定律#xff1a;所有事物都与其他事务相关#xff0c;但是近…有一篇大神的文章写得非常的具体
https://xg1990.com/blog/archives/222
下面写下一些学习笔记
1、关于克里金插值的基本原理
克里金插值来源于地理学它的前提是地理学第一定律所有事物都与其他事务相关但是近处的事物比远处的事物更相关。这里就引出克里金插值的
第一个核心思想在一定范围内两点属性值的差异性不相关性与二者间距离在一定距离范围内成正相关两点属性值相关性与二者之间距离正相关。
可以举一个简单的例子苏州和上海的气温差异要小于苏州和北京的差异但是也有距离限制比如苏州和伦敦的差异就不一定比苏州和北京的差距大
这里有一个关键结论一点的属性值可以由其周围点的属性值推导出
2、克里金插值的数学形式
做一些基础的假设已知点坐标及属性值已知待求点的坐标已知点间距离矩阵已知点-待求点之间的距离矩阵。根据克里金插值的第一核心思想可以得到空间内某未知点x0,y0估计值Z0是部分已知点值的加权和。可被表示为 其中n代表点的数量zi代表第 i 个已知点的属性值 wi 代表第 i 个已知点的权重Zi 都是已知点但是 wi 目前还不是一个已知量。具体需要一个什么样的 wi 才能求得准确的 Z0 呢为了得到更加准确的估测结果权重系数 wi 应使得点x0,y0处估计值与真实值之差最小 我们需要做如下期望
2.1 估计值 Z0 与真实值 Z0 之间方差尽可能小
估计值 Z0 与真实值 Z0 之间的差值为期望为0 估计值 Z0 与真实值 Z0 之间方差最小 即该假设使得结果尽可能准确 2.2 空间偏差的方差不变
为了维持空间是平稳的空间任意一点处的值 z z(x,y)由区域平均值 c 和随机偏差 R(x,y) 组成其中偏差的方差均为常数。 写成期望的形式 这个假设就保证了可以通过已知点值估算未知点值空间平稳、具备一定的趋势。对于 z(x0,y0)估计值 Z0、真实值 Z0、使用n个已知属性值的领域点求其值每个点有一直的属性值zi权重 wi 经过文章所示的推导过程克里金推导然后稍加求逆可以得到权重系数的表达矩阵 2.3 半方差的矩阵
对 Zi Z(xi,yj)Zj Z(xi,yj)其半方差函数可被表示为 这个是非常关键的一个点。他表示了第 i , j 个领域点之间的半方差第0个领域点表示待求点本身。Z1....Zn已知y11 和 ynn 可计算得到。但真实值 Z0是未知的y10yn0是未知的。如何求y10...yn0是需要解决的问题
2.4 解算的思路建立拟合函数
从上面的分析可知我们已知每个点之间的距离有了每个点之间距离组成的n*n的矩阵并且还已知维持空间稳定的半方差矩阵是不是可以理解为我既知道了已知点和未知点之间的相对距离还知道了一部分这个空间不同点之间的一个约束关系那么是否拟合出一个非常接近真实情况的函数形式来表达不同点在不同距离的变化趋势呢
建立一个拟合函数 该拟合函数的 d 即为待求点的位置假如将其代入到上面的权重系数矩阵可得 如果这个拟合函数已知那么 w1- wn 就可以被求得。最后通过方程即可完成 z0 的估计。 2.5 为什么可以通过函数拟合
在进行拟合函数构建之前需要讨论的是为什么可以通过函数拟合这个时候回到地理学第一定律在一定范围内两点差异随着两点间距离增大而增大并且空间是稳定的也就是说他必然是存在趋势的因此这个拟合是可以接受的。在已经有了半方差 矩阵之后我们可以通过这种潜在的约束和已知点的集合去对真实的空间变化趋势进行拟合。
2.6 如何构建拟合函数呢
根据克里金差值的原理两点差异随着两点间距离增大而增大因此这个函数一定是一个增函数。然后这个差异的变化趋势是在一定范围内满足的因此这个函数一定存在某个极值或者边界。这个函数自变量增大到一定的值后因变量不再增大。已有的模型可以用来借鉴主其中球型模型、指数模型、高斯模型都比较符合。 为了得到更加精确的结果采用先分组再拟合的方式已知n*n组距离-半方差关系 按照距离d大小排序这n*n组关系再分组得到m个距离-半方差关系组g1、g2...gm。计算每组平均距离和平均半方差 根据每组的平均距离和平均半方差构建拟合函数。
如果我们整体回顾一下整个克里金差值模型的思路为了完成 Z0真实值的估计我们进行了如下步骤
第一步建立期望——差值最小、方差最小得到估计值表达式要得到估计值就必须要求得权重系数矩阵 第二步建立权重系数表达矩阵——依据优化目标构建权重系数表达矩阵。y11-ynn均已知。要求解权重系数矩阵必须先求解y10...yn0矩阵。 第三步求解y10...yn0矩阵——构建拟合函数 y f(d)要构建拟合函数必须进行模型选择及数据处理
第四步选择模型、数据分组分段函数拟合
所以最关键的就是这些数据的处理、模型的选择、函数的拟合下面就针对这个进行学习。
3 函数的拟合
首先有这个么数据其二维平面的分布如下图所示 第一步求取距离与半方差求每个点之间的距离和半方差半方差和距离都是75*75的矩阵最后将其变为一个二维的矩阵即代码 d_r[d_r;[dij(i,j),rij(i,j)]]; 其中dij是距离rij是半方差。
S[0.700000000000000,59.6000000000000;2.10000000000000,82.7000000000000;4.70000000000000,75.1000000000000;4.80000000000000,52.8000000000000;5.90000000000000,67.1000000000000;6,35.7000000000000;6.40000000000000,33.7000000000000;7,46.7000000000000;8.20000000000000,40.1000000000000;13.3000000000000,0.600000000000000;13.3000000000000,68.2000000000000;13.4000000000000,31.3000000000000;17.8000000000000,6.90000000000000;20.1000000000000,66.3000000000000;22.7000000000000,87.6000000000000;23,93.9000000000000;24.3000000000000,73;24.8000000000000,15.1000000000000;24.8000000000000,26.3000000000000;26.4000000000000,58;26.9000000000000,65;27.7000000000000,83.3000000000000;27.9000000000000,90.8000000000000;29.1000000000000,47.9000000000000;29.5000000000000,89.4000000000000;30.1000000000000,6.10000000000000;30.8000000000000,12.1000000000000;32.7000000000000,40.2000000000000;34.8000000000000,8.10000000000000;35.3000000000000,32;37,70.3000000000000;38.2000000000000,77.9000000000000;38.9000000000000,23.3000000000000;39.4000000000000,82.5000000000000;43,4.70000000000000;43.7000000000000,7.60000000000000;46.4000000000000,84.1000000000000;46.7000000000000,10.6000000000000;49.9000000000000,22.1000000000000;51,88.8000000000000;52.8000000000000,68.9000000000000;52.9000000000000,32.7000000000000;55.5000000000000,92.9000000000000;56,1.60000000000000;60.6000000000000,75.2000000000000;62.1000000000000,26.6000000000000;63,12.7000000000000;69,75.6000000000000;70.5000000000000,83.7000000000000;70.9000000000000,11;71.5000000000000,29.5000000000000;78.1000000000000,45.5000000000000;78.2000000000000,9.10000000000000;78.4000000000000,20;80.5000000000000,55.9000000000000;81.1000000000000,51;83.8000000000000,7.90000000000000;84.5000000000000,11;85.2000000000000,67.3000000000000;85.5000000000000,73;86.7000000000000,70.4000000000000;87.2000000000000,55.7000000000000;88.1000000000000,0;88.4000000000000,12.1000000000000;88.4000000000000,99.6000000000000;88.8000000000000,82.9000000000000;88.9000000000000,6.20000000000000;90.6000000000000,7;90.7000000000000,49.6000000000000;91.5000000000000,55.4000000000000;92.9000000000000,46.8000000000000;93.4000000000000,70.9000000000000;94.8000000000000,71.5000000000000;96.2000000000000,84.3000000000000;98.2000000000000,58.2000000000000];
Y[34.1000000000000;42.2000000000000;39.5000000000000;34.3000000000000;37;35.9000000000000;36.4000000000000;34.6000000000000;35.4000000000000;44.7000000000000;37.8000000000000;37.8000000000000;43.9000000000000;37.7000000000000;42.8000000000000;43.6000000000000;39.3000000000000;42.3000000000000;39.7000000000000;36.9000000000000;37.8000000000000;41.8000000000000;43.3000000000000;36.7000000000000;43;43.6000000000000;42.8000000000000;37.5000000000000;43.3000000000000;38.8000000000000;39.2000000000000;40.7000000000000;40.5000000000000;41.4000000000000;43.3000000000000;43.1000000000000;41.5000000000000;42.6000000000000;40.7000000000000;42;39.3000000000000;39.2000000000000;42.2000000000000;42.7000000000000;40.1000000000000;40.1000000000000;41.8000000000000;40.1000000000000;40.9000000000000;41.7000000000000;40.8000000000000;38.7000000000000;41.7000000000000;40.8000000000000;38.7000000000000;38.6000000000000;41.6000000000000;41.5000000000000;39.4000000000000;39.8000000000000;39.6000000000000;38.8000000000000;41.6000000000000;41.3000000000000;41.2000000000000;40.5000000000000;41.5000000000000;41.5000000000000;38.9000000000000;39;39.1000000000000;39.7000000000000;39.7000000000000;40.3000000000000;39.5000000000000];nsize(S,1);
rijzeros(n,n);
dijzeros(n,n);
d_r[];for i1:nfor j1:nrij(i,j)0.5*(Y(i,1)-Y(j,1)).^2; % 半方差dij(i,j)sqrt(sum((S(i,:)-S(j,:)).^2));%欧式距离 绝对距离
% dij(i,j)sum(abs(S(i,:)-S(j,:)));%曼哈顿距离(在优化问题中易于线性化)d_r[d_r;[dij(i,j),rij(i,j)]]; % dij是i x j y是的距离rij是半方差两个dr就是持续往下叠叠乐下end
end
将其画为点状图的形式可得
绘图代码
S[0.700000000000000,59.6000000000000;2.10000000000000,82.7000000000000;4.70000000000000,75.1000000000000;4.80000000000000,52.8000000000000;5.90000000000000,67.1000000000000;6,35.7000000000000;6.40000000000000,33.7000000000000;7,46.7000000000000;8.20000000000000,40.1000000000000;13.3000000000000,0.600000000000000;13.3000000000000,68.2000000000000;13.4000000000000,31.3000000000000;17.8000000000000,6.90000000000000;20.1000000000000,66.3000000000000;22.7000000000000,87.6000000000000;23,93.9000000000000;24.3000000000000,73;24.8000000000000,15.1000000000000;24.8000000000000,26.3000000000000;26.4000000000000,58;26.9000000000000,65;27.7000000000000,83.3000000000000;27.9000000000000,90.8000000000000;29.1000000000000,47.9000000000000;29.5000000000000,89.4000000000000;30.1000000000000,6.10000000000000;30.8000000000000,12.1000000000000;32.7000000000000,40.2000000000000;34.8000000000000,8.10000000000000;35.3000000000000,32;37,70.3000000000000;38.2000000000000,77.9000000000000;38.9000000000000,23.3000000000000;39.4000000000000,82.5000000000000;43,4.70000000000000;43.7000000000000,7.60000000000000;46.4000000000000,84.1000000000000;46.7000000000000,10.6000000000000;49.9000000000000,22.1000000000000;51,88.8000000000000;52.8000000000000,68.9000000000000;52.9000000000000,32.7000000000000;55.5000000000000,92.9000000000000;56,1.60000000000000;60.6000000000000,75.2000000000000;62.1000000000000,26.6000000000000;63,12.7000000000000;69,75.6000000000000;70.5000000000000,83.7000000000000;70.9000000000000,11;71.5000000000000,29.5000000000000;78.1000000000000,45.5000000000000;78.2000000000000,9.10000000000000;78.4000000000000,20;80.5000000000000,55.9000000000000;81.1000000000000,51;83.8000000000000,7.90000000000000;84.5000000000000,11;85.2000000000000,67.3000000000000;85.5000000000000,73;86.7000000000000,70.4000000000000;87.2000000000000,55.7000000000000;88.1000000000000,0;88.4000000000000,12.1000000000000;88.4000000000000,99.6000000000000;88.8000000000000,82.9000000000000;88.9000000000000,6.20000000000000;90.6000000000000,7;90.7000000000000,49.6000000000000;91.5000000000000,55.4000000000000;92.9000000000000,46.8000000000000;93.4000000000000,70.9000000000000;94.8000000000000,71.5000000000000;96.2000000000000,84.3000000000000;98.2000000000000,58.2000000000000];
Y[34.1000000000000;42.2000000000000;39.5000000000000;34.3000000000000;37;35.9000000000000;36.4000000000000;34.6000000000000;35.4000000000000;44.7000000000000;37.8000000000000;37.8000000000000;43.9000000000000;37.7000000000000;42.8000000000000;43.6000000000000;39.3000000000000;42.3000000000000;39.7000000000000;36.9000000000000;37.8000000000000;41.8000000000000;43.3000000000000;36.7000000000000;43;43.6000000000000;42.8000000000000;37.5000000000000;43.3000000000000;38.8000000000000;39.2000000000000;40.7000000000000;40.5000000000000;41.4000000000000;43.3000000000000;43.1000000000000;41.5000000000000;42.6000000000000;40.7000000000000;42;39.3000000000000;39.2000000000000;42.2000000000000;42.7000000000000;40.1000000000000;40.1000000000000;41.8000000000000;40.1000000000000;40.9000000000000;41.7000000000000;40.8000000000000;38.7000000000000;41.7000000000000;40.8000000000000;38.7000000000000;38.6000000000000;41.6000000000000;41.5000000000000;39.4000000000000;39.8000000000000;39.6000000000000;38.8000000000000;41.6000000000000;41.3000000000000;41.2000000000000;40.5000000000000;41.5000000000000;41.5000000000000;38.9000000000000;39;39.1000000000000;39.7000000000000;39.7000000000000;40.3000000000000;39.5000000000000];figure(4)
scatter(S(:,1),S(:,2)) 第二步数据的预处理。从上面的图可以看到 距离-方差 散点分布图非常的混乱没有任何规律常规的是根据距离进行点的分组但是目前这种情况也没办法直接根据不同的距离进行分组。因此需要对数据进行一个基本的处理查找了一下大佬们所用的方法用k-means聚类对距离进行聚类聚到一类的距离和半方差均取平均值再进行拟合
k-means聚类是对基础的数据进行分簇处理英文直译是这个词
K-means 算法是一种迭代聚类算法它将数据集划分为 k 个不同的簇每个簇包含具有相似特征的数据点。该算法的基本原理如下 初始化 随机选择 k 个数据点作为初始质心簇的中心。 分配 将每个数据点分配到离其最近的质心所在的簇。 更新质心 计算每个簇的新质心即该簇中所有数据点的平均值。 重复 重复步骤 2 和 3直到簇不再发生变化或者达到预定的迭代次数。
K-means 的目标是最小化每个数据点与其所属簇的质心之间的平方欧氏距离的总和。这被称为簇内平方和inertia。整个过程的目标是最小化以下损失函数 其中
n 是数据点的数量。k 是簇的数量。xi 是第 i 个数据点。cj 是第 j 个簇的质心。
K-means 算法的迭代过程不保证全局最优解因此对于不同的初始质心可能会得到不同的聚类结果。通常采用多次运行算法即设置 Replicates 参数来找到最佳的聚类。
在matlab中。这个函数的一般语法如下
[idx, C] kmeans(X, k)其中
X 是包含观测值的矩阵每行代表一个观测值每列代表一个特征。k 是要分为的簇的数量。idx 是一个列向量包含了每个观测值所属的簇的索引。C 是一个矩阵包含了每个簇的质心的坐标。
以下是一个简单的示例
% 创建一个示例数据集
X [randn(50,2)2; randn(50,2)-2];% 使用 kmeans 进行聚类将数据分为 k 2 个簇
k 2;
[idx, C] kmeans(X, k);% 绘制聚类结果
scatter(X(:,1), X(:,2), 20, idx, filled);
hold on;
plot(C(:,1), C(:,2), rx, MarkerSize, 10, LineWidth, 2);
hold off;
title(K-Means Clustering Result);在这个例子中X 是一个包含两个聚类的示例数据集。kmeans 函数将数据分为两个簇并返回每个点的簇索引和质心坐标。最后使用散点图可视化聚类结果每个簇用不同的颜色表示质心用红色叉号标记。 对目标数据函数进行处理的结果为绘图代码与结果如下所示
figure(8)
scatter(d_r(:,1),d_r(:,2),.)
hold on
scatter(d_r_divided(:,1),d_r_divided(:,2));
hold on;
plot(d_r_divided(:,1),d_r_divided(:,2), rx, MarkerSize, 10, LineWidth, 2);
hold on 此时的K值选择的是20也可以选择为40或者更多
第三步选择拟合模型进行曲线拟合。选择高斯模型为目标模型
高斯半变异函数是半变异函数的一种类型通常用于描述数据在空间上的平滑性和相关性。
高斯半变异函数的数学表达式通常如下 其中
γ(h) 是半变异函数表示两个点之间的半变异。σ2 是变异的总体方差。ℎ是两个点之间的空间 lag间隔或距离。a 是半变异函数的范围参数控制了半变异函数的横向展布。
高斯半变异函数的图形呈钟形bell-shaped或类似正态分布的形状。该函数在 lag 达到一定距离 a 时趋于平稳表示空间相关性在这个距离上基本消失这和地理学相符。在地统计学和地理信息系统 (GIS) 中通过半变异函数的拟合和分析可以对空间数据的相关性结构进行建模有助于空间插值、空间预测和地理空间数据分析。
这个时候关键一步调用matlab拟合工具箱
cftool(d_r_divided(:,1),d_r_divided(:,2))
cftool 是 MATLAB 中的 Curve Fitting 工具是一个交互式工具用于拟合曲线和表面拟合。这个工具使用户能够使用不同的拟合模型如线性、多项式、指数、对数等来拟合数据并通过可视化工具来调整拟合参数以达到最佳拟合效果。操作过程如下图所示1选择数据源2选择高斯模型3选择高斯阶数4观察曲线拟合效果 符合的话就按照公式与系数撰写一个单独的函数
function [y] Guassian_func_3(x)a1 7.81 ;b1 50.88 ;c1 29.73 ;a2 6.685 ;b2 99.01 ;c2 17.29 ;[m,n]size(x);
yzeros(m,n);
for i1:mfor j1:ny(i,j)a1*exp(-((x(i,j)-b1)/c1)^2) a2*exp(-((x(i,j)-b2)/c2)^2);end
end
end
拟合结果如下所示当然是阶数越高这个效果越好但是计算量也越大这个按照合适的定。 第四步根据拟合曲线进行预测。
在完成了 距离-半方差 曲线的拟合后需要返回来求取具体的权重系数方程然后对Z值进行估计因此将距离值带回到高斯过程。
求取权重系数方程之前先要通过拟合的模型求解y10...yn0矩阵这个要用到拟合出来的f函数d就是待预测点与已知点之间的距离把待预测点与所有已知点之间的半方差就可以求出来了这个半方差矩阵就是y10...yn0矩阵
K_line 75;
X gridsamp([0 0;100 100], K_line ); % gridsamp就是线性的生成0-100的40个点一共是几维度数组就是40^n个数目然后排列组合。%这个排列组合的形式是线性生成 0-10040个点%线性第二列生成0-10040个点然后这第一列的40个点与第二列40个点之间排列组合。[m,~]size(X);
YXzeros(m,1);xX(i,:); % 第i行所有的值
nsize(S,1);
dixsqrt(sum((repmat(x,n,1)-S).^2,2));
rixGuassian_func_3(dix); 然后依据如下公式求解权重系数的矩阵先构建这个rij最右侧列11110最下一行11110组合代码 temp[[rij,ones(size(rij,1),1)];[ones(1,size(rij,1)),0]];
最后除一下就可以得到权重系数矩阵lamba利用最开始的那个预测方程即可完成预测 预测结果如下图所示
