八十八队数模A题.doc

“工大出版社杯”第十六届西北工业大学数学建模竞赛暨全国大学生数学建模竞赛选拔赛题目A题密封号2015年5月4日剪 切 线密封号2015年5月4日 自动化 学院 第 八十八 队队员1队员2队员3姓名宋亚龙李浩男张建红班级090214030902140310021401目 录 摘要 2 模型分析 3 模型假设 3 模型参数的假设 3 模型建立 3 模型的求解 8 电阻率数据插值加密及成像问题摘 要 实际中一个形状不规则，质量和密度不均匀物体内部和表面的每一个点的电阻率很难或不能实际测量，这对我们研究物理问题形成了阻碍，为了解决显示的物理问题，我们必须知道或估测出很接近的某一点的电阻率数值，因此，本题意欲通过数学建模求出特定点的电阻率数值，并且证明插值后的极值在同一个位置求得，进而求得加密后的每一点电阻率数据和原网格及其加密后网格的平均值和标准差并对两种算法进行评估，然后对加密后的进行颜色图示表示，观察与原图的对比，最后对两种方法的效果进行定量表示。第一问建立三次线性插值模型和反加权插值模型计算出定点的电阻率数值，然后通过数学证明得到极值点未发生移动。第二问运用matlab软件进行插值拟合，将步距由10调成1，借助计算机求得与原数据每一个点对应的电阻率数值，并进而求得插值前后的平均值和标准差，以及进行评估。第三问运用matlab进行绘图，并分别令z=0和50,得到平面二维的颜色图示，直观地看出电阻率在整个物体以及一个界面上变化情况。第四问建立结果和定值的关系，定量的分析这两个加密方法的优缺。关键字：MATLAB软件 数据插值 曲线插值拟合三维模型 插值法 颜色图示一 模型分析1问题背景 物体的电阻是一个很常用的物理量，有很大的运用价值，尤其是在物理学中，它涉及到很多方面，因此知道物体的电阻率很重要，但是实际物体通常呈现不规则形装和不均匀的质量和密度，导致我们很难知道整个物体各个地方的电阻率，因此迫切地需要一种办法来计算出物体各个地方的电阻率，此题便是基于这个问题所提出来的。2. 问题分析 对于问题一，用附件中给出的数据，用matlab插值法建立三维模型， 对于问题二，基于问题一给出的两种方法，在matlab里计算出网格大小为1*1*1时的电阻率数据，再用均值法计算出加密前后的平均值和标准差。 对于问题三，在matlab中画出z为定值时的三维立体图像并着色。 而对于问题四，定量的比较出两种方法的效果，并做评估。二 模型假设（1）物体的外部形状不随时间的变化而变化。（2）电阻率不随时间和温度等外界因素变化而变化。（3）取样点的数据较好地反映了该物体的电阻率。（4）测量的个别数据对整体没有影响。（5）物体内部的电阻率连续变化。三 模型参数的假设x表示物体内部待求点的横坐标y表示物体内部待求点的纵坐标z表示物体内部待求点的竖坐标表示物体内部待求点处的电阻率四 模型建立在此对一维插值法做如下简介：插值：求过已知有限个数据点的近似函数。拟合：已知有限个数据点，求近似函数，不要求过已知数据点，只要求在某种意义下它在这些点上的总偏差最小。（1）插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。而面对一个实际问题，究竟应该用插值还是拟合，有时容易确定，有时则并不明显。 插值方法：下面介绍一种基本的、常用的插值：拉格朗日多项式插值拉格朗日多项式插值插值多项式 用多项式作为研究插值的工具，称为代数插值。其基本问题是：已知函数f(x)在区间a,b上n +1个不同点x0,x1,xn处的函数值yi=f(xi)(i=0,1,n)，求一个至多n 次多项式n(x)=a0+a1x+anxn (1)使其在给定点处与f(x) 同值，即满足插值条件n(xi)=f(xi)=yi(i=0,1,n) (2)称为插值多项式，xi(i=0,1,n) 称为插值节点，简称节点，a,b称为插值区间。从几何上看，n次多项式插值就是过n+1个点(xi,f(xi)(i=0,1,n), 作一条多项式曲线 y=n(x)近似曲线 y=f(x). n次多项式（1）有n +1个待定系数，由插值条件（2）恰好给出n +1个方程 a0+a1x0+a2x02+anx0n=y0a0+a1x1+a2x12+anx1n=y1. (3) a0+a1xn+a2xn2+anxnn=yn 记此方程组的系数矩阵为 A ，则1 x0 x02 x0n 1 x1 x12 x1ndet(A)=1 x2 x22 x2n 1 xn xn2 xnn是范德蒙特(Vandermonde)行列式。当x0,x1,xn互不相同时，此行列式值不为零。因此方程组（3）有唯一解。这表明，只要n+1个节点互不相同，满足插值要求（2）的插值多项式（1）是唯一的。插值多项式与被插函数之间的差 Rn(x)=f(x)-n (x) 称为截断误差，又称为插值余项。当 f (x)充分光滑时 Rn(x)=f(x)-Ln(x)=fn+1()/(n+1)!n+1(x), (a,b) 其中n+1（x）=j=0n(x-xi)。现在建立模型如下：1. 基于断层图像分割的三维匹配插值如下图所示 由已知断层图像(xi,yj,zk)和（xi,yj,zk+1）插值出位于它们之间的新的一层图像（xi,yj,zk）。不失一般性，假设k层图像距离新插值图像较远，则对于差值图像上的每个像素点（xi,yj,z），在k层图像上，以(xi,yj,zk)为中心取宽度WW的窗口，窗口中的每一个点(xn,ym,zk)和点(xi,yj,z)的连线于k+1层图像有交点（xi,yj,zk+1），这些k层和k+1层上的对应点对构成一组初始匹配点，共WW对。（1） 计算候选匹配点位置 匹配点队都在选择的窗口内，因此匹配窗口的大小决定匹配点对的数目。窗口宽度W是一个固定值，为了保证对称，一般选择奇数宽窗，如3x3,5x5等。实验证明，窗宽选择比较大时，不到计算量大，算法效率低，而且插值出的效果也差强人意。 在K层断层图像中，匹配点的坐标为：xm(k)=xi+ 2m+1-W2yn(k)=yi+ 2n+1-W2 相对于K层的匹配点，在 K+1层断层图像中，匹配点对的对应坐标为：xm(k+1)=xi+intd2d1(xW-m-1(k)-xi)ym(k+1)=yi+intd2d1(yW-m-1(k)-yi)m,n=0,1,2,W-1 由上述的四个公式，就可以建立(xm(k),yn(k)和(xm(k+1),yn(k+1)匹配点对。其中d1=z-zk 为点（xi,yj,z）到(xi,yj,zk)的距离，d2=zk+1-z 为点到(xi,yj,zk+1)的距离。（2） 确定最佳匹配点 匹配插值的关键是从众多的匹配点对中确定最佳匹配点对。此文提出最佳匹配点对应该满足以下条件： 匹配点对连线短的距离较小； 匹配点对具有同一个旋转方向； 匹配点对之间的电阻率值相同或者相近； 匹配点对的电阻率梯度值相近；根据这四个条件，构造了如下的向量函数表示一对匹配点对(x,y,z)和(x,y,z)的匹配程度： R(x,y,z, x,y,z)=(x,y,z)- (x,y,z)i+D(x,y,z)- D(x,y,z)j +(x,y,z)- (x,y,z)k+x-x+y-y l其中(x,y,z)，D(x,y,z)和(x,y,z)分别表示点(x,y,z)的电阻率值，电阻率梯度大小和梯度方向，(x,y,z)，D(x,y,z)和(x,y,z)分别表示点(x,y,z)的电阻率值，电阻率梯度大小和梯度方向。x-x+y-y是匹配点对(x,y,z)和(x,y,z)在水平面上的投影的距离。这里，电阻率梯度D(x,y,z)是取D(x,y,z)的模。 D(x,y,z)=12(x+1,y,z)- (x-1,y,z)i+12(x,y+1,z)- (x,y-1,z)j+12(x,y,z+1)- (x,y,z-1)k 由（5）式可得：|R(xm(k),yn(k),zk,xm(k+1),yn(k+1),zk+1|=(xm(k),yn(k),zk)- (xm(k+1),yn(k+1),zk+1)2+(D(xm(k),yn(k),zk)- D(xm(k+1),yn(k+1),zk+1)2+( (xm(k),yn(k),zk)- (xm(k+1),yn(k+1),zk+1)2+(xm(k)- xm(k+1)2+(yn(k) -yn(k+1)2根据以上的最佳匹配点对需要满足的四个条件，如果点（xp(k),yp(k),zk）和（xp(k+1),yp(k+1),zk+1）构成最佳匹配点对，则此两点必须满足：|R(xp(k),yp(k),zk,(xp(k+1),yp(k+1),z(k+1)|=minm,n=0,1,W-1| R(xp(k),yp(k),zk,(xp(k+1),yp(k+1),z(k+1)| 有上述步骤求出窗口中的最佳匹配点对后，可由这两点的电阻率值作线性插值，得到目标点(xi,yi,z)的电阻率值：(xi,yj,z)=1d1+d2(xp(k),yq(k),zk)d2+(xp(k+1),yq(k+1),zk+!)d1(3)（3） 给予分割的三维匹配插值 给出要插值的断层图像到上下断层图像的距离分别是d1和d2； 求出最佳匹配点对(xp(k),yq(k),zk)和(xp(k+1),yq(k+1),zk+1)； 当上层区域大于下层区域时，将 作为缩放因子对于上层图像的区域缩小， 作为插值数据；当上层区域小于下层区域时，则将d1d1+d2作为缩放因子对上层图像的区域放大，作为差值数据。上数算法关键之处就是确定最家匹配点对的位置，如果在相同的密度物质内的区域里，就利用匹配插值，在不同的密度物质区域里时，就是对相应的区域进行缩放处理，然后把缩放区域作为插值数据，得到新的断层图像数据。这样就可以避免单纯的匹配插值所带来的误差，提高了插值图像的质量。 为验证基于断层图像分割的匹配插值算法的有效性，采用一组z为定值时的切面进行实验。一般用标准差来衡量插值结果的好坏，其和分别是估值与真值为总点数。两种算法的标准差见下表：2 利用反距离权重法测定三维物体某点电阻率(8)理论依据：由于我们在实际条件下测定空间中某一点的电阻率，因此我们利用反距离权重法进行测算，即对所要测定的某一点，其周边离它越近的点的电阻率与其电阻率相差越小（或者说二者距离越近，两者的电阻率关联性越强），离它越远的点的电阻率的值与其电阻率相差越大。那么，我们希望，在进行计算时，距离测量点越远的点对该点的的测算影响越小，于是我们利用 =i=1ni *wiwi=（1/dip）/（i=1n1/dip）其中，i 为测量点周围已知的第i个整数点的电阻率，wi为权重，di=xi-x+yi-y+zi-zp为设定的影响指数，基于上述论述，在相同条件下，p的取值越大，精确度越高，考虑实际计算量，我们取p=4。当测量点周围的已知点取得越多时，测量值越精确，因此对于n，我们希望越大越好。数据检验：已知量待求量x100101010y10001010z4040303041199.3010.00560.00210.02199.3308x4040404045.8y-40-30-30-30-32.7z6060708068.218.80328.22109.716.6183199.4419五 模型求解问题一： （1）计算： 第一种方法是基于断层图像分割的三维匹配插值法。对于此方法的模型建立及其推倒和运用上述已经给出此处不再赘述，现在给出公式： (xi,yj,z)=1d1+d2(xp(k),yq(k),zk)d2+(xp(k+1),yq(k+1),zk+!)d1计算在（45.8，-32.7,68.2）处的电阻率： 第二种方法是反距离权重法，上述亦已经给出证明和公式，此处只简要说明这种方法的主要思想是一个点的电阻率和它周围的点有关联，并且周围的点对其的影响随着距离的增大而减小。因此可以通过周围的点来求出特定的点。现在给出公式：=i=1ni *wi计算在（45.8，-32.7,68.2）处的电阻率：199.4419（2）证明：根据问题，我们采用Matlab，用程序运行的结果证明了极值与位置不变。clc, cleara=load(data3D_.txt);x0=a(:,1); y0=a(:,2); z0=a(:,3); r0=a(:,4);r00=griddata(x0,y0,z0,r0,45.8,-32.7,68.2) %计算定点的插值结果x=-100:100; y=-100:100; z=0:100;x,y,z=meshgrid(x,y,z); %数据网络化r1=griddata(x0,y0,z0,r0,x,y,z); %111网格线性插值r1=r1(:); %把插值结果展开成列向量mu1=mean(r1), s1=std(r1) %计算均值和标准差r2=griddata(x0,y0,z0,r0,x,y,z,natural); %111网格“Triangulation-based natural neighbor interpolation”插值r2=r2(:); mu2=mean(r2), s2=std(r2) %计算均值和标准差问题二：原网格数据的平均值：197.2689 标准差：13.73325（1） 算法一得出的1*1*1的网格的数据：平均值：194.4135标准差：14.0464（2） 算法二得出的1*1*1的网格的数据：平均值：196.6598标准差：13.2348（3） 对两种计算方法的评估： 问题三： 已知z=40时的原图像和加密后的图像分别如下：记处电阻率的最小值，中间值和最大值分别记为。对每一幅需要对比显示效果的图，请将最小值置为纯蓝色（RGB为(0,0,255)），中间值置为纯绿色（RGB为(0,255,0)），最大值置为纯红色（RGB为(255,0,0)），中间数值采用过渡的颜色，可自行设计。我们设计的电阻率与R，G，B像素值之间的对应关系为clc, cleara=load(data3D_.txt);x0=a(:,1); y0=a(:,2); z0=a(:,3); r0=a(:,4);xb=-100:10:100; yb=-100:10:100; %初始的网格划分n=length(xb)ind=find(z0=40);x40=x0(ind), y40=y0(ind) %提出z=4