论文2-通过三维数据插值进行电阻率的拟合及估计

.z.摘要在实际问题中，我们无法对一个物体进行密集测量，一般等间隔的选取部分点进行测量，这时候就需要采用数学方法，根据已知位置的数据计算未知位置的数据。本文根据题目要求建立了两种合理的差值加密模型——反距离加权插值模型、克里金插值模型，对文件中给出的坐标网格大小为10m*10m*10m的三维电阻率数据进行插值，得到坐标网格大小为1m*1m*1m的数据。我们借助多种数学软件，处理大量数据并对其进行可视化处理，绘制出颜色图反映插值效果，根据插值拟合特点对插值效果给出定量指标。针对问题一，首先对所给出的三维的电阻率数据进行分析，根据不同空间插值方法的优势，结合已知数据的网格化特点，确定插值模型——反距离加权插值模型、克里金插值模型。反距离加权插值模型算法简单易于实现，但对权重函数的选择十分敏感，而网格化的数据恰好规避了这一劣势。克里金插值是在变异函数理论及结构分析基础上，进行无偏、最优估计的一种方法，不仅考虑了观测点与待估计点的相对位置，而且考虑了各观测点之间的相对关系，插值效果较好。利用两种方法分别计算出空间*点(45.8,-32.7,68.2)处的电阻率数值。反距离加权插值法结果为Z1=196.18；克里金插值结果为Z2=194.13。针对问题二：利用问题一中两种插值方法，分别计算出网格大小为1m*1m*1m的三维电阻率数据（由于数据量过大，具体数据见电子版1、2），分析两种模型的计算流程可知，克里金算法计算过程较为复杂，数据量较大，计算时间与占用内存大，而反距离加权插值计算过程较为简单，计算时间明显短于克里金方法，但计算精度较低。将原网格数据及两种方法加密网格后数据利用E*cel计算出平均值与标准差，见下表。平均值E方差S2原始数据197.26893713.73324911反距离加权插值数据196.54990515.70823291克里金差值数据197.64752115.47241092针对问题三：对插值加密网格后的数据，利用颜色图直观展示。将电阻率最小值置为纯蓝色，中间值置为纯绿色，最大值置为纯红色，中间数值采用过渡的颜色，利用matlab对切片Z=0,50分别给出原数据，两种方法加密网格后数据的颜色对比图。另外分别对切片*=82，Y=47,Z=88，给出两种加密方法得到数据的颜色图。通过同一切面颜色图的比较，直观反映出原数据、反距离加权插值加密、与克里金插值加密的插值效果。针对问题四：利用插值加密前后电阻率均值的偏离程度以及电阻率沿*，y，z三个轴向变化趋势的偏离程度给出定量指标，定量反应两种插值方法的差值效果。关键词：反距离加权插值变差函数加权系数克里金插值一、问题重述在实际问题中，由于技术与测量成本等原因，只能等间隔的选取部分点进行测量。实际中需要采用数学方法获得更多位置数据，根据已知位置的数据计算未知位置的数据。根据坐标网格大小为10m*10m*10m的空间三维体电阻率数据，通过插值加密后获得网格大小为1m*1m*1m的三维电阻率数据，必须保证插值后的电阻率数据极值不变，且坐标位置相同。另外，必须保证插值后的三维成像结果与插值前的成像结果形态基本一致，只是前者像素更高。（1）请给出符合条件的两种计算方法，并给出相应数学公式，证明该方法插值后的数据极值与位置不变。分别计算出空间*点(45.8,-32.7,68.2)处的电阻率数值。(2)利用（1）中两种方法分别计算出网格大小为1m*1m*1m的三维电阻率数据，给出计算流程，并对两种方法的计算量或计算复杂性进行评估。同时给出原网格数据及两种方法加密网格后数据的平均值与标准差。(3)对加密网格后的直观效果，可采用颜色图展示。对每一幅需要对比显示效果的图，请将最小值置为纯蓝色（RGB为(0,0,255)），中间值置为纯绿色（RGB为(0,255,0)），最大值置为纯红色（RGB为(255,0,0)），中间数值采用过渡的颜色，可自行设计。采用这种方法，请对切片Z=0,50分别给出原数据，两种方法加密网格后数据的颜色对比图。另外分别对切片*=82，Y=47,Z=88，请给出两种加密方法得到数据的颜色图。（4）对不同的插值加密方法的效果，给出定量的指标，分别计算出两种不同加密方法得到数据的指标值，并给出评价。二、模型假设1、基本假设：1)假设data3D_.t*t所给电阻率数据背景值真实；2)被测电阻各向同性；3)不考虑环境因素对电阻率的影响；2、针对反距离加权插值法：未知值的点受较近控制点的影响比较远控制点的影响更大；3、针对克里金插值法：本征假设：设Z（*）为空间点*处的属性值即电阻率，Z(*)与Z(*+h)之间的相关性不依赖于在电阻内的特定位置。实际中只需假设其1,2阶矩存在且平稳。当区域变化量的协方差函数不存在时，只需满足式（1）、（2），即满足本征假设。①区域化变量Z(*)的增量[Z(*)-Z(*+h)]的数学期望为0,即:E[Z(*)-Z(*+h)]=0(*,h)(1)②增量[Z(*)-Z(*+h)]的方差存在且平稳,即:Var[Z(*)-Z(*+h)]=E[Z(*)-Z(*+h)]-{E[Z(*)-Z(*+h)]}2=E[Z(*)-Z(*+h)]2,*,h(2)三、符号说明空间中待插值点待插值点与其邻域内第i个点之间的距离待插值点P处的电阻率空间点*处的电阻率滞后距实验变差函数理论变差函数块金常数基台值变程原始数据与加密后数据期望的偏离程度沿*轴向加密数据与原始数据变化率的偏离程度沿y轴向加密数据与原始数据变化率的偏离程度沿z轴向加密数据与原始数据变化率的偏离程度四、模型的分析、建立与求解4.1问题一问题分析由于技术和测量成本等问题，我们只能获得有限个电阻率数据，而为了获得更多位置的数据，则需要利用空间插值加密的方法建立电阻率空间分布模型，并保证插值前后电阻率数据极值与位置不变。反距离加权插值法设电阻中待插点为P(*,y,z),P点邻域内有已知散乱点Qi(*i,yi,zi),i=1,2,…,n，利用反距离加权插值法对P点的属性值（即电阻率）进行插值，其插值原理是以插值点与样本点间的距离为权重进行加权平均，离插值点越近的样本点赋予的权值越大，权的大小是距离k（k一般取2）次方的倒数。即：Zp=（3）其中：di为待插点与其邻域内第i个点之间的距离。即：（4）4.1.3本题中克里金的研究对象(电阻率)是区域化变量,而表征区域化变量空间相关性的工具就叫做变差函数。使用已知点的属性值所计算的变差函数叫做实验变差函数。数学表达式为：（5）其中,表示在*轴向上的实验变差函数值;表示在点*处的属性值;表示在点处的属性值。（1）计算实验变差函数在本征假设的前提下，区域化变量Z(*)的增量[Z(*)-Z(*-h)]只依赖于向量h，并不依赖于位置*，被向量h分割的每对数据可以看成一对随机变量的一个不同实现。这时，对于每个滞后距h，我们可以算出的值，其计算公式为：（6）而本题需要计算不同轴向方向的多个实验变差函数，所以将公式变形为：（7）（2）理论变差函数球状模型理论变差函数模型是未知的，所以需要从有限的空间已知样本中去估计，通过不同的滞后距计算不同的，再用理论变差函数球状模型来拟合。球状模型的基本公式为：（8）式中：为块金常数，为基台值，a为变程。接近原点时，变差函数呈线性形状，在变程处达到基台值。其图形如下：图4-1-3-1（3）拟合实验变差函数假设对于不同的滞后距已经计算出相应的实验变差函数值，又对每个滞后距参与计算的数据对的数目为Ni。球状模型表达式：（9）当时，令，得到：，则有：（10）目标函数取最小值时，可求得其中参数。目标函数分别对求偏导，得到：（11）将上述方程组整理并写成矩阵形式，有：（12）解得：（13）（4）求解待插值点函数值（14）问题求解：（1）由于以上两种插值加密方法的插值条件保证了插值函数在插值节点处函数值与原数据相同，故插值后的电阻率数据极值与位置不变。（2）由反距离插值方法得：Z（45.8,-32.7,68.2）=196.18由克里金插值方法得：Z（45.8,-32.7,68.2）=194.134.2问题二：三维电阻率数据分别采用问题一中所建立的反距离插值模型与克里金插值模型计算出网格大小为1m*1m*1m的三维电阻率数据见1、2。计算流程将反距离插值法与克里金插值法的计算流程分别用图（4-2-1）、图（4-2-2）表示。待插值节点P（*，y，z）待插值节点P（*，y，z）P点邻域内Qi(*i,yi,zi)图4-2-1反距离加权插值法计算实验变差函数计算已知数据点*i与*0（0,0,0）之间距离di由di找出不同滞后距h，及其个数N（h）理论变差函数基本模型加权最小二乘拟合理论变差函数利用计算图4-2-2克里金插值法计算量与计算复杂性评估（1）在反距离插值算法中：需要根据目标点到已知数据点间的距离的平方的倒数来进行加权求出加权系数，进而求出待测点的属性值，随着数据量及拟和精度的提高的增加，计算量以及存储数据所需要的空间会成倍增加，计算数据所需时间较长，所需内存空间较大。在计算复杂度方面主要是进行乘法和加法运算，算法的时间复杂度和空间复杂度较大。且差值效果一般。（2）在克里金插值算法中：主要是对变差函数的计算以及求解加权系数。在计算变差函数时，由于所给区域化变量为三维空间内的函数，在进行计算时，主要由三个参数：搜索方向，最大搜索半径，基本滞后距。本文采用等距搜索，由于数据量较大搜索次数较多，采用简化方法后，以原点为基本点，减少了搜索次数。在运算中主要解一个矛盾方程组（克里金方程组），求解加权系数过程与反距离插值相同，因此克里金插值法德计算量比反距离插值法大得多，在时间和空间复杂度上均是反距离插值法的数倍，但插值效果较好。平均值与标准差平均值E方差S2原始数据197.26893713.73324911反距离加权插值数据196.54990515.70823291克里金差值数据197.64752115.47241092表4-2-34.3问题三问题分析对原始数据进行插值加密后，可以得到更加密集的数据，利用matlab绘制颜色图，使加密后绘制图的像素明显提高。绘制颜色图（1）Z=0原数据：图4-3-2-1反距离加权插值后数据：图4-3-2-2克里金插值后数据：图4-3-2-3（2）Z=50原数据：图4-3-2-4反距离加权插值后数据：图4-3-2-5克里金插值后数据：图4-3-2-6（3）*=82反距离插值后数据：图4-3-2-7克里金插值后数据：图4-3-2-8（4）Y=47反距离加权插值后数据：图4-3-2-9克里金插值后数据：图4-3-2-10（5）Z=88反距离加权插值后数据：图4-3-2-11克里金插值后数据：图4-3-2-12问题四：4.1问题分析：数据插值加密效果主要从数据均值和数据变化趋势两方面分析，故给出定量衡量指标与。4.2衡量指标：(1)指标：表征原始数据与差值加密后数据期望的偏离程度设Z（*,y,z）为三维空间点（*,y,z）的电阻率值。其中表示原始数据的期望，表示差值加密后的数据期望。(2)指标S：可分别估计出沿*，y，z三个方向上插值加密数据与原数据变化率的偏离程度。设三维空间点（*,y,z）取z为插值*围内任意常数z0，则在z=z0平面内：设沿*轴的插值节点为，两点间的一阶差商为，利用差值方法计算点处的函数值。表示在z=z0平面内，沿y轴方向，插值加密数据与原数据变化率的偏离程度。设沿y轴的插值节点为，两点间的一阶差商为，利用差值方法计算点处的函数值。表示在z=z0平面内，沿*轴方向，插值加密数据与原数据变化率的偏离程度。设三维空间点（*,y,z）取*为插值*围内任意常数*0，则在*=*0平面内：设沿z轴的插值节点为，两点间的一阶差商为，利用差值方法计算点处的函数值。表示在*=*0平面内，沿z轴方向，插值加密数据与原数据变化率的偏离程度。4.3数据分析：平均值分析均值E平面原数据反距离加权插值克里金插值*=0204.4205.18204.4724109266y=0204.5205.78204.2246618829z=0199.0197.22199.0055699984表4-31：整体差值效果分析：反距离加权插值：2.07克里金插值：=0.79由平均值计算结果可知克里金插值结果较好，与原数据相差较小，说明整体差值结果与原数据吻合度较高，而反距离插值效果相对较差。2：数据整体变化趋势分析：为了更好地衡量差值效果本模型提出用相邻点的一阶差分与原数据的一阶差分进行对比，通过原数据点相邻区间内差值点的差商与原数据的差商的差值作为偏离程度，并绘出各轴向偏离程度图线进行对比如下：由图线可以看出在各轴向上都有偏离程度均有少许震荡，但都稳定在零偏附近，说明克里金插值效果与电阻率的实际分布规律较为吻合，进一步说明差值效果较好。五、模型评估5.1模型优点本文采用的差值模型均为精确差值方法，较全面的运用了已经观测点数据，能够较好的与实际模型相吻合，相对精确地再现了模型的真实值，在整体效果与变化趋势方面均能与已有数据相吻合，通过matlab绘图所得结果与原数据绘图结果基本吻合，具有较高精确度，属于较精确的差值模型。5.2缺点：在计算量和计算复杂度方面，需要根据目标点到已知数据点间的距离的平方的倒数来进行加权求出加权系数，进而求出待测点的属性值，随着数据量及拟和精度的提高的增加，计算量以及存储数据所需要的空间会成倍增加，计算数据所需时间较长，所需内存空间较大。在计算复杂度方面主要是进行乘法和加法运算，算法的时间复杂度和空间复杂度较大。5.3模型改进：（1）在计算权重时，我们采用的是简化模型，计算精度较低，若采用其他克里金插值如指数克里金插值，泛克里金插值法，高斯模型，差值效果可能会有所提高。有待进一步测试。（2）在计算变差函数进行搜索时，考虑到算法复杂度，采用模型过于简化，可能会影响计算结果的精度，下一步可才能考有关论文，对搜索方法进行进一步改善提高插值的精确度和缩小计算步骤，降低算法复杂度。（3）在进行算法的推导方面，由于时间仓促，没有进行大量修改，因此程序运行时间较长，后期可对相关算法程序进行改进。六、参考文献[1]翟进乾.克里金插值方法在煤层分布检测中的应用研究.[学位论文].**：**理工大学.2008[2]蔡元菲.快速克里金算法的研究与实现.[学位论文].**：电子科技大学.2013[3]靳国栋，*衍聪，牛文杰等.距离加权反比插值与克里金插值法的比较.**大学学学报.2003[4]阮晓青，周义仓等数学建模引论..高等教育，2005.7[5]清源计算机工作室编.MATLAB基础与应用..机械工业3：相关算法程序1、Z切面程序formatlongl=1;for*2=-100:1:100;m=1;fory2=-100:1:100forz1=0:10:100for*1=-100:10:100fory1=-100:10:100D((*1+110)/10,(y1+110)/10,(z1+10)/10)=1/((*2-*1)^2+(y2-y1)^2+(88-z1)^2)^3;endendendfori=1:11;forj=1:21fork=1:21C(j,k,i)=B(k+(j-1)*21+(i-1)*21*21);endendendf=0;fori=1:4851;f=f+D(i);endE=C.*D;e=0;fori=1:4851;e=e+E(i);endr(l,m)=e/f;m=m+1;endl=l+1;end>>V=r;pcolor(V);figure(gcf);gca=pcolor(*,y,V);set(gca,'LineStyle','none')title('Z=0数据成像')title('Z=88数据成像')2、*切面程序l=1;forz2=0:100m=1;fory2=-100:1:100forz1=0:10:100for*1=-100:10:100fory1=-100:10:100D((*1+110)/10,(y1+110)/10,(z1+10)/10)=1/(sqrt((82-*1)^2+(y2-y1)^2+(z2-z1)^2))^3;endendendfori=1:11;forj=1:21fork=1:21C(j,k,i)=B(k+(j-1)*21+(i-1)*21*21);endendendf=0;fori=1:4851;f=f+D(i);endE=C.*D;e=0;fori=1:4851;e=e+E(i);endr(l,m)=e/f;m=m+1;endl=l+1;end>>V=r;[*,y]=meshgrid(1:201,1:101);pcolor(*,y,r);gca=pcolor(*,y,r);set(gca,'linestyle','none');>>title('*=82数据成像')3、Y切面程序l=1;forz2=0:100m=1;for*2=-100:1:100forz1=0:10:100for*1=-100:10:100fory1=-100:10:100D((*1+110)/10,(y1+110)/10,(z1+10)/10)=1/((*2-*1)^2+(47-y1)^2+(z2-z1)^2);endendendfori=1:11;forj=1:21fork=1:21C(j,k,i)=B(k+(j-1)*21+(i-1)*21*21);endendendf=0;fori=1:4851;f=f+D(i);endE=C.*D;e=0;fori=1:4851;e=e+E(i);endr(l,m)=e/f;m=m+1;endl=l+1;end4、平均值程序：B=data3D_(:,4);formatlongl=1;forz2=0:1:100;m=1;for*2=-100:1:100forz1=0:10:100for*1=-100:10:100fory1=-100:10:100D((*1+110)/10,(y1+110)/10,(z1+10)/10)=1/((*2-*1)^4+(1-y1)^4+(z2-z1)^4);endendendfori=1:11;forj=1:21fork=1:21C(j,k,i)=B(k+(j-1)*21+(i-1)*21*21);endendendf=0;fori=1:4851;f=f+D(i);endE=C.*D;e=0;fori=1:4851;e=e+E(i);endr(l,m)=e/f;m=m+1;endl=l+1;ends1=sum(sum(r))/203015、变化趋势偏离度：Z轴向B=data3D_(:,4);fori=1:11forj=1:21A(i)=B((10+(i-1)*21)*21+11);endendformatlongm=1;forz2=0:1:100forz1=0:10:100for*1=-100:10:100fory1=-100:10:100D((*1+110)/10,(y1+110)/10,(z1+10)/10)=1/((1-*1)^2+(1-y1)^2+(z2-z1)^2)^3;endendendfori=1:11;forj=1:21fork=1:21C(j,k,i)=B(k+(j-1)*21+(i-1)*21*21);endendendf=0;fori=1:4851;f=f+D(i);endE=C.*D;e=0;fori=1:4851;e=e+E(i);endr(m)=e/f;m=m+1;endfori=2:11C1(i-1)=(A(i)-A(i-1))/10;endfori=2:101C2(i-1)=(r(i)-r(i-1));endfori=1:10forj=1:10C3((i-1)*10+j)=(C2((i-1)*10+j)-C1(i));endend6、*轴向fori=1:21forj=1:21D1(i)=B((i-1)*21+11)endendformatlongm=1;for*2=-100:1:100forz1=0:10:100for*1=-100:10:100fory1=-100:10:100D((*1+110)/