1、0“工大出版社杯”第十六届西北工业大学数学建模竞赛暨全国大学生数学建模竞赛选拔赛题目A 题密封号 2015 年 5 月 4 日剪 切 线密封号 2015 年 5 月 4 日自动化 学院 第 八十八 队队员 1 队员 2 队员 3姓名 宋亚龙 李浩男 张建红班级 09021403 09021403 100214011目 录摘要 2模型分析 3模型假设 3模型参数的假设 3模型建立 3模型的求解 82电阻率数据插值加密及成像问题摘 要实际中一个形状不规则，质量和密度不均匀物体内部和表面的每一个点的电阻率很难或不能实际测量，这对我们研究物理问题形成了阻碍，为了解决显示的物理问题，我们必须知道或估测出
2、很接近的某一点的电阻率数值，因此，本题意欲通过数学建模求出特定点的电阻率数值，并且证明插值后的极值在同一个位置求得，进而求得加密后的每一点电阻率数据和原网格及其加密后网格的平均值和标准差并对两种算法进行评估，然后对加密后的进行颜色图示表示，观察与原图的对比，最后对两种方法的效果进行定量表示。第一问建立三次线性插值模型和反加权插值模型计算出定点的电阻率数值，然后通过数学证明得到极值点未发生移动。第二问运用 matlab 软件进行插值拟合，将步距由 10 调成 1，借助计算机求得与原数据每一个点对应的电阻率数值，并进而求得插值前后的平均值和标准差，以及进行评估。第三问运用 matlab 进行绘图，
3、并分别令 z=0 和 50,得到平面二维的颜色图示，直观地看出电阻率在整个物体以及一个界面上变化情况。第四问建立结果和定值的关系，定量的分析这两个加密方法的优缺。关键字：MATLAB 软件 数据插值 曲线插值拟合三维模型 插值法 颜色图示3一 模型分析1问题背景物体的电阻是一个很常用的物理量，有很大的运用价值，尤其是在物理学中，它涉及到很多方面，因此知道物体的电阻率很重要，但是实际物体通常呈现不规则形装和不均匀的质量和密度，导致我们很难知道整个物体各个地方的电阻率，因此迫切地需要一种办法来计算出物体各个地方的电阻率，此题便是基于这个问题所提出来的。2. 问题分析对于问题一，用附件中给出的数据，
4、用 matlab 插值法建立三维模型，对于问题二，基于问题一给出的两种方法，在 matlab 里计算出网格大小为1*1*1 时的电阻率数据，再用均值法计算出加密前后的平均值和标准差。对于问题三，在 matlab 中画出 z 为定值时的三维立体图像并着色。而对于问题四，定量的比较出两种方法的效果，并做评估。二 模型假设（1）物体的外部形状不随时间的变化而变化。（2）电阻率不随时间和温度等外界因素变化而变化。（3）取样点的数据较好地反映了该物体的电阻率。（4）测量的个别数据对整体没有影响。（5）物体内部的电阻率连续变化。三 模型参数的假设x 表示物体内部待求点的横坐标y 表示物体内部待求点的纵坐标
5、z 表示物体内部待求点的竖坐标 表示物体内部待求点处的电阻率四 模型建立在此对一维插值法做如下简介：插值：求过已知有限个数据点的近似函数。拟合：已知有限个数据点，求近似函数，不要求过已知数据点，只要求在某种意义下它在这些点上的总偏差最小。 （1）4插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。而面对一个实际问题，究竟应该用插值还是拟合，有时容易确定，有时则并不明显。插值方法：下面介绍一种基本的、常用的插值：拉格朗日多项式插值拉格朗日多项式插值插值多项式用多项式作为研究插值的工具，称为代数插值。其基本问题是：已知函数f(x)在区间 a,b上
6、n +1 个不同点 x0,x1,xn处的函数值 yi=f(xi)(i=0,1,n)，求一个至多 n 次多项式n(x)=a 0+a1x+anxn (1)使其在给定点处与 f(x) 同值，即满足插值条件n(x i)=f(xi)=yi(i=0,1,n) (2)称为插值多项式，x i(i=0,1,n) 称为插值节点，简称节点， a,b称为插值区间。从几何上看， n 次多项式插值就是过 n+1 个点(x i,f(xi)(i=0,1,n), 作一条多项式曲线 y=n(x)近似曲线 y=f(x).n 次多项式（1）有 n +1 个待定系数，由插值条件（2）恰好给出 n +1 个方程a0+a1x0+a2x02
7、+anx0n=y0a0+a1x1+a2x12+anx1n=y1 (3) a0+a1xn+a2xn2+anxnn=yn记此方程组的系数矩阵为 A ，则1 x0 x02 x0n1 x1 x12 x1ndet(A)= 1 x2 x22 x2n1 xn xn2 xnn是范德蒙特(Vandermonde)行列式。当 x0,x1,xn互不相同时，此行列式值不为零。因此方程组（3）有唯一解。这表明，只要 n+1 个节点互不相同，满足插值要求（2）的插值多项式（1）是唯一的。插值多项式与被插函数之间的差Rn(x)=f(x)- n (x)称为截断误差，又称为插值余项。当 f (x)充分光滑时Rn(x)=f(x)
8、-Ln(x)=fn+1()/(n+1)! n+1(x), (a,b)其中 n+1（x）= =0(x-xi)。现在建立模型如下：1. 基于断层图像分割的三维匹配插值如下图所示5由已知断层图像 (x i,yj,zk)和 （x i,yj,zk+1）插值出位于它们之间的新的一层图像 （x i,yj,zk） 。不失一般性，假设 k 层图像距离新插值图像较远，则对于差值图像上的每个像素点 （x i,yj,z） ，在 k 层图像上，以 (x i,yj,zk)为中心取宽度 WW 的窗口，窗口中的每一个点(x n,ym,zk)和点(x i,yj,z)的连线于 k+1 层图像有交点 （x i,yj,zk+1） ，
9、这些 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+int (xW-m-1(k)-xi)21ym(k+1)=yi+int (
10、yW-m-1(k)-yi)21m,n=0,1,2,W-1由上述的四个公式，就可以建立(x m(k),yn(k)和(x m(k+1),yn(k+1)匹配点对。其中 d1=z-zk 为点（x i,yj,z）到 (x i,yj,zk)的距离，d 2=zk+1-z 为点到(x i,yj,zk+1)的距离。（2） 确定最佳匹配点匹配插值的关键是从众多的匹配点对中确定最佳匹配点对。此文提出最佳6匹配点对应该满足以下条件： 匹配点对连线短的距离较小； 匹配点对具有同一个旋转方向； 匹配点对之间的电阻率值相同或者相近； 匹配点对的电阻率梯度值相近；根据这四个条件，构造了如下的向量函数表示一对匹配点对(x,y,
11、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+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,y,z)和(x,y,z)在水平面上的投影的(-)+(-)距离。这里，电阻率梯度 D(x,y,z)是取D(x,y,z)的模。D(x,y,z)= (x+1,y,
12、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)k12由（5）式可得：|R(xm(k),yn(k),zk,xm(k+1),yn(k+1),zk+1|=(x m(k),yn(k),zk)- (x m(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根据以上的最佳匹配点对需要满足的四个条件，如
13、果点（x p(k),yp(k),zk）和（x p(k+1),yp(k+1),zk+1）构成最佳匹配点对，则此两点必须满足：|R(xp(k),yp(k),zk,(xp(k+1),yp(k+1),z(k+1)|=R(xp(k),yp(k),zk,(xp(k+1),yp(k+1),z(k+1)|minm,n=0,1,-1|有上述步骤求出窗口中的最佳匹配点对后，可由这两点的电阻率值作线性插值，得到目标点(x i,yi,z)的电阻率值：(x i,yj,z)= (x p(k),yq(k),zk)d2+(x p(k+1),yq(k+1),zk+!)d1(3)11+2（3） 给予分割的三维匹配插值 给出要插值
14、的断层图像到上下断层图像的距离分别是 d1和 d2； 求出最佳匹配点对(x p(k),yq(k),zk)和(x p(k+1),yq(k+1),zk+1)； 当上层区域大于下层区域时，将 作为缩放因子对于上层图像的区域缩小， 作为插值数据；当上层区域小于下层区域时，则将 作为缩放因子对11+27上层图像的区域放大，作为差值数据。上数算法关键之处就是确定最家匹配点对的位置，如果在相同的密度物质内的区域里，就利用匹配插值，在不同的密度物质区域里时，就是对相应的区域进行缩放处理，然后把缩放区域作为插值数据，得到新的断层图像数据。这样就可以避免单纯的匹配插值所带来的误差，提高了插值图像的质量。为验证基于
15、断层图像分割的匹配插值算法的有效性，采用一组 z 为定值时的切面进行实验。一般用标准差来衡量插值结果的好坏，其和分别是估值与真值为总点数。两种算法的标准差见下表：2 利用反距离权重法测定三维物体某点电阻率 (8)理论依据：由于我们在实际条件下测定空间中某一点的电阻率，因此我们利用反距离权重法进行测算，即对所要测定的某一点，其周边离它越近的点的电阻率与其电阻率相差越小（或者说二者距离越近，两者的电阻率关联性越强） ，离它越远的点的电阻率的值与其电阻率相差越大。那么，我们希望，在进行计算时，距离测量点越远的点对该点的的测算影响越小，于是我们利用= ini=1 i *wi=（1/d ip）/（ ip
16、）=11/其中， i 为测量点周围已知的第 i 个整数点的电阻率，w i为权重，di= (-)+(-)+(-)p 为设定的影响指数，基于上述论述，在相同条件下，p 的取值越大，精确度越高，考虑实际计算量，我们取 p=4。当测量点周围的已知点取得越多时，测量值越精确，因此对于 n，我们希望越大越好。数据检验：已知量 待求量x 10 0 10 10 10y 10 0 0 10 10z 40 40 30 30 41 199.301 0.0056 0.0021 0.02 199.3308x 40 40 40 40 45.8y -40 -30 -30 -30 -32.7z 60 60 70 80 68.
17、2 18.803 28.22 109.7 16.6183 199.44198五 模型求解问题一：（1）计算： 第一种方法是基于断层图像分割的三维匹配插值法。对于此方法的模型建立及其推倒和运用上述已经给出此处不再赘述，现在给出公式：(x i,yj,z)= (x p(k),yq(k),zk)d2+(x p(k+1),yq(k+1),zk+!)d111+2计算在（45.8，-32.7,68.2）处的电阻率：第二种方法是反距离权重法，上述亦已经给出证明和公式，此处只简要说明这种方法的主要思想是一个点的电阻率和它周围的点有关联，并且周围的点对其的影响随着距离的增大而减小。因此可以通过周围的点来求出特定的
18、点。现在给出公式：= ini=1 i *计算在（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,
19、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标准差：
20、14.0464（2） 算法二得出的 1*1*1 的网格的数据：平均值：196.6598标准差：13.2348（3） 对两种计算方法的评估：问题三：已知 z=40 时的原图像和加密后的图像分别如下：记 处电阻率的最小值，中间值和最大值分别记为 。40z 321,r9对每一幅需要对比显示效果的图，请将最小值置为纯蓝色（RGB 为(0,0,255)） ，中间值置为纯绿色（RGB 为(0,255,0)） ，最大值置为纯红色（RGB 为(255,0,0)） ，中间数值采用过渡的颜色，可自行设计。我们设计的电阻率与 R，G，B 像素值之间的对应关系为.),(25,0)( 3231rrrR.),(25,)(
21、 32311rrrG.,0,),()(32211rBclc, 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=40 对应的 x,y 坐标r40=r0(ind) %提出 z=40 的电阻率r1=min(r40), r2=median(r40), r3=max(r40) %z=40 时电阻率的最小，中间，最大值b(n
22、,n,1)=0; b(n,n,2)=0; b(n,n,3)=0; %电阻率对应的 RGB 矩阵的初始化Rr=(r)255/(r3-r2)*(r-r2).*(rr2 Rr0=Rr(r40) %红色像素的取值Gr0=Gr(r40) %绿色像素的取值Br0=Br(r40) %蓝色像素的取值for k=1:length(ind)i=find(xb=x0(ind(k); j=find(yb=y0(ind(k);b(i,j,1)=Rr0(k); b(i,j,2)=Gr0(k); b(i,j,3)=Br0(k);endb=uint8(b); imshow(b) %显示题目中的第一幅图片imwrite(b,xibei.jpg) %把 3 维矩阵，即彩色图片写到 jpg 文件。利用 matlab 进行绘图可以得到z=0 时的加密前后的图像如下图：10z=50 时的加密前后图像如下所示：11问题四：由于本队技术力量不够，无法给出定量指标。本文用到的参考文献：插值与拟合三次样条插值在工程中的运用基于断层图像分割的三维匹配插值方法研究一元 n 次多项式根的展开公式及其求根算法12三维放射治疗计划系统的研究数码相机定位-2008 高教社杯全国大学生数学建模竞赛2011 年数学建模 A 题城市表层土壤重金属污染反距离权重法的运用
