插值算法(六):趋势面插值(Trend Surface)

标签:
trend趋势面插值实现c |
分类: GIS |
http://s8/mw690/6316e2afgdce6fa09e6e7&690Surface)" TITLE="插值算法(六):趋势面插值(Trend Surface)" />何时使用趋势插值法
使用趋势插值法可获得表示感兴趣区域表面渐进趋势的平滑表面,此种产值法适用于一下几种情况
1、感兴趣区域的表面在各位置间出现渐变时,可以将该表面与采样点拟合,例如,工业区的污染情况
2、检查或排除长期趋势或全局趋势的影响,此类情况下,次啊用方法通常为趋势面分析。
此外,计算得出的表面对于异常值(极高值和极低值)非常敏感,尤其边缘处。
趋势插值法的类型
趋势插值法共有两种基本类型:线性和逻辑型
线性趋势
线性趋势面插值法用于创建浮点型栅格。它将通过多项式回归将最小二乘表面与各输入点进行拟合。使用线性选项可以控制用于拟合表面的多项式阶数。要理解趋势面法的线性选项,请考虑一阶多项式。一阶多项式趋势面插值法将对平面与一组输入点进行最小二乘拟合。
利用趋势插值法可创建平滑表面。生成的表面几乎不能穿过各原始数据点,因为对整个表面执行的是最佳拟合。如果所用多项式阶数高于一阶,插值器所生成的栅格的最大值和最小值可能会超过输入要素数据输入文件中的最大值和最小值。
逻辑型趋势
可生成趋势面的逻辑型选项适用于预测空间中给定的一组位置(x,y)处某种现象存在与否(以概率的形式)。z值是仅会产生两种可能结果的分类随机变量,例如,濒临灭绝的物种存在与否。生成的两种z值分别编码为1和0。逻辑型选项可以根据值为0和1的歌像元值创建连续的概率网格。
可以使用最大可能性估计直接计算出非线性概率表面模型,而无需现将该模型转换成线性形式。
- 求解:http://s10/bmiddle/6316e2afgdce70227ff99&690Surface)" TITLE="插值算法(六):趋势面插值(Trend
Surface)" /> - //////////////////////////////////////////////////////////////////////////
- ////趋势面插值
- ////dpx,pdy,dpz为采样点数据
- ////ndp为采样点数
- ////ipx,ipy,ipz为插值点数据
- ////nip为插值点数
- ///////////最小二乘法/////////////////////////////////////////////////////
- ////
A B C - ////|--------------| |----| |---|
- ////| n |
X | Y | | b0 | | Z | - ////| X | X*X | XY |*| b1 |=| XZ|
- ////| Y | X*Y | Y*Y| | b2 | | YZ|
- ////|--------------| |----| |---|
- //// A的逆矩阵与C矩阵相乘可得到B的解
- //////////////////////////////////////////////////////////////////////////
- void CInterpolationTrend::InterpolateByTrend(double *dpx,double *dpy,double *dpz,int ndp,double *ipx,double *ipy,double *ipz,int nip)
- {
- int n=i_eTrendPower;
- int ML=GetParamerCount(n);
- ///初始矩阵
- double **A =new double*[ML];
- double *B=new double[ML];
- double *C=new double[ML];
- double **temp=new double*[ML];
- double dl,dr,dd;
- //////////////////////////////////////////////////////////////////////////
- ///初始化矩阵A
- for(int i=0; i
- A[i] = new double[ML];
- }
- for (int i=0;i
- {
- for (int j=0;j
- {
- A[i][j]=0.0;
- }
- }
- //////////////////////////////////////////////////////////////////////////
- ///第一个矩阵A
- ///第三个矩阵C
- //////////////////////////////////////////////////////////////////////////
- for (int i=0;i
- {
- int xr=GetXPower(n,i);
- int yr=GetYPower(n,i);
- for (int j=0;j
- {
- int xc=GetXPower(n,j);
- int yc=GetYPower(n,j);
- dl=0.0;
- for (int ps=0;ps
- {
- dl+=pow(dpx[ps],xr)*pow(dpy[ps],yr)*pow(dpx[ps],xc)*pow(dpy[ps],yc);
- }
- A[i][j]=A[j][i]=dl;
- }
- dr=0.0;
- for (int ps=0;ps
- {
- dr+=pow(dpx[ps],xr)*pow(dpy[ps],yr)*dpz[ps];
- }
- C[i]=dr;
- }
- //////////////
- //对A矩阵求逆
- temp=InverseMatrix(A,ML);
- //求出系数矩阵B
- MultiplyMatrix(temp,C,ML,B);
- //根据插值点和系数矩阵求出预测值
- for (int i=0;i
- {
- dd=0.0;
- for (int j=0;j
- {
- int xc=GetXPower(n,j);
- int yc=GetYPower(n,j);
- dd+=pow(ipx[i],xc)*pow(ipy[i],yc)*B[j];
- }
- ipz[i]=dd;
- }
- delete []A;delete []B;delete[]C;
- }