2023年高教社杯全国大学生数学建模竞赛题目-CT系统参数标定及成像

2023西交数模第一次模拟赛

数学建模论文首页

选题A

队伍编号696

班级学号姓名

队长钱学森73班2173714071孔令辉

队员1钱学森73班2171311098杨峥

队员2电气642160400092李佳航

2023年6月30日

摘要

本题针对一种二维CT获得样品内部结构信息的工作方式及成像原理，意在

通过借助已知结构样品进行参数标定，消退系统误差，而后对未知结构的样品进

行成像，得出该未知介质的相关信息。

问题一中，要求通过标定模板的相关，确定CT系统的旋转中心、探测器单

元之间的距离以及该CT系统运用的X射线的180个方向。本文利用标定模板几

何参数，通过条数、以及探测器单元间距相等等信息，首先计算出探测器单元之

间的距离为0.2778mm。依据X射线与标定模板的几何特性，以椭圆短轴所在直

线为x轴，长轴所在直线为y轴建立直角坐标系，得旋转中心在所建立的坐标系

中的坐标为(-9.0396,5.7060)。最终通过近似确定X射线180个方向的旋转角度

具有高度线性相关性，依据部分确定数据拟合出整体旋转角度。前十五组旋转角

度为：

29.649331.002731.557232.647633.679534.6478

35.648336.649037.647538.648639.648840.6491

41.648742.649743.6492

问题二中，要求通过已求得的标定参数，确定未知介质在正方形托盘中的位

置、几何形态和吸取率等信息，并具体给出所要求十个位置的吸取率。本文依据

CT层析成像原理，利用逆拉东变换作出重构图像，并利用Excel中数据分布计算

出原位置介质的相关性质。所求十个位置的吸取率为：

序号12345

吸取率0.01850.9701-0.00011.17031.0287

序号678910

吸取率1.44251.2815-0.0074-0.00250.0194

问题三中，与问题二计算方法相像得十个位置的吸取率为:

序号12345

吸取率0.07362.53057.13270.03320.8015

序号678910

吸取率3.14875.57140.05308.13620.0272

问题四中，要求对问题一中的标定模型进行改进以减小误差并增加稳定性，

本文利用在问题一求解过程中遇到的问题进行思索，首先适当增大模板减小偶然

误差，其次做出投射图像后应简洁找到极值，并且图像应有确定的对称性；图像

扫描后应尽量少地得出重复数据。

关键词：CT层析成像Radon变换与逆变换吸取率MATLAB算法

1.问题重述

1.1问题背景

CT（ComputedTomography）可以在不破坏样品的状况下，利用样品对射线能量的吸取特

性对生物组织和工程材料的样品进行断层成像，由此获得样品内部的结构信息。

1.2相关信息

本题介绍了一种二维CT系统，平行入射的X射线垂直于探测器平面，每个探测器单元

看成一个接收点，且等距排列。X射线的放射器和探测器相对位置固定不变，整个放射-接

收系统绕某固定的旋转中心逆时针旋转180次。对每一个X射线方向，在具有512个等距

单元的探测器上测量经位置固定不动的二维待检测介质吸取衰减后的射线能量，并经过增益

等处理后得到180组接收信息。并且，为消退安装误差，须要对安装好的CT系统进行参数

标定，即借助于已知结构的样品（称为模板）标定CT系统的参数，并据此对未知结构的样

品进行成像。题目附件中供应了标定模板的几何信息，接受信息，待测介质的接收信息以及

图3所给位置的相应数据。

1.3须要解决的问题

1.3.1问题一

在正方形托盘上放置两个匀整固体介质组成的标定模板，模板的几何信息如图2所示，

请依据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元

之间的距离以及该CT系统运用的X射线的180个方向。

1.3.2问题二

附件3是利用上述CT系统得到的某未知介质的接收信息。利用第一问中得到的标定参

数，确定该未知介质在正方形托盘中的位置、几何形态和吸取率等信息。另外，请具体给出

图3所给的10个位置处的吸取率。

1.3.3问题三

利用给出上述CT系统得到的另一个未知介质的接收信息。利用第一问中得到的标定参

数，给出该未知介质的相关信息。另外，请具体给出图3所给的10个位置处的吸取率。

1.3.4问题四

分析第一问中参数标定的精度和稳定性。在此基础上自行设计新模板、建立对应的标定

模型，以改进标定精度和稳定性，并说明理由。

图1.CT系统小意图图2.模板小意图（单位：mm）

2.问题分析

2.1问题一分析

结合图2和附件1表中数据，可以首先计算出CT系统探测器个数和模板长

度度量的比值，运用程序1.1可以得出模板的几何形态如图2.1.1,可大致认为它

是对称的。对于附件2,由于对-180个方向尚无清晰地相识，首先用同样的方

式做出数据分布图2.1.2,视察到图像比较平滑，因此认为按表格的依次180个

方向是相邻较密、不错位的。（图中有色区域表示该点有吸取率，蓝色部分表示

吸取率大于100）

Figure1—0X.Figure1一口X

文件（D蒙球1）ffS（Y）瓶入（1）IHCD如IQ«n（w）WK（H）、文件(F)«Si(E)fifi(V)JLM)ia(T)£30(0)«C(W)HMn(H)»

d）k•、-、考公*/•息口国♦・口□dd琳♦、二学@里4•③□3■□

250

200

150

100

50

0

020406080100120140160180

图2.1.1附件1的数据分布图图2.1.2附件2的数据分布图

在图2.1.2中易观查到有一红色条形图案，这是在不同的方向扫描到圆形时

留下的，可以此为突破口首先求出探测器单元之间的距离。求出探测器单元之间

的距离与增益比率之后，可以依据几何关系，自己设立坐标系并通过数学运算计

算出旋转中心。最终在建立的坐标系内，将待求解的180个旋转方向转换成X

射线的斜率进行数学运算。

2.2问题二分析

题目给出了未知介质的接收信息，要求出介质的相关信息，可以搜寻相应的

数学模型，将附件中给出的依据射线条数与旋转角度列成的表格，一一对应为相

应坐标点的吸取状况，从而依据各坐标点的不同性质，还原回该未知介质的几何

信息与吸取率等信息。

2.3问题三分析

问题三与问题二类似，可以大致看出数据分布更具有一般性，不简洁描述出

未知介质的相关信息，可以通过图形大致描绘出介质相关信息。

2.4问题四分析

问题四要求分析题目所给的二维CT系统，设计新的标定模板以提高原系统

的精确度与稳定性。可以搜寻相关资料，依据第一问的求解思路与求解过程，以

规避求解过程中因标定模板自身性质而出现的误差为原则进行思路拓展，以设计

高精度与稳定性的标定模板。

3.模型假设

♦假设在射线经过介质时能量只损失在介质中，及不考虑衍射等现象;

♦假设附件中所给出数据是正确的、可以干脆利用的；

♦假设旋转中心在相邻两条射线的中间直线上；

4.模型的建立与求解

4.1问题一模型的建立与求解

4.1.1建立模型

从题目中可以得出，由于x射线之间得间距相等，不管x射线怎么旋转，穿

过托盘上圆的X射线条数应当是大致相等的。可将穿过圆形标定模板的X射线与

模板建立模型示意图如下：

图4.1.1.IX射线穿过圆形标定模板示意图

如图4.1.1.1所示可以设穿过圆的在圆心两旁的射线被圆所截的距离分别为

X1，X2,X3……和丫1，丫2，丫3……。

求解旋转中心时，以椭圆形标定模板短轴所在直线为X轴，长轴所在直线为

Y轴建立直角坐标系，将托盘进行划分。由附件二的数据分布图的不对称性可知,

旋转中心应在对称轴某一侧，大致确定旋转中心方位后，依据旋转中心两侧探测

器个数不变且同一探测器接收的与距旋转中心的距离不变具体确定旋转中心的

位置。

求解X射线旋转角度时,设穿过椭圆的最边缘的射线到椭圆中心的距离为R,

取椭圆中心为原点，模板的对称中心为x轴建立平面直角坐标系，X射线所在直

线的斜率为ko

4.1.2符号说明

d相邻探测器之间的距离

N对应于圆形模板的射线条数或探测器个数

U标定模板的吸取率

X处理数据时的增益率

4)圆形模板的直径

4.1.3模型求解

求探测器单元间距离：

通过MATLAB编程求出计算出的条形带平均涉及探测器个数为

万=28.8023,沿直径方向上的平均吸取率〃=14.1766；由图二圆形模板直径

6=8mm,计算出探测器之间的平均距离：

Q

------mm«0.2778mm

N28.8023

探测系统平均增益率:

〃14.1766

1.7721

O--8-

求旋转中心：

分析图2.1.2可知，由于蓝色区域仅分布在后部分角度范围内，因此估计旋

转中心的位置在对称轴的某一侧，红色区域分为两部分时表示该角度下由两部分

射线分别照射经过椭圆形和圆形；为了探讨便利，现对正方形托盘做出如下划分:

■1—□X

MBM£iwooaxu）工AID闻agvaoaomon，

Q♦、-/・3□0■□

2SO

200

150

100

50

0

801001902002SO300

图4.1.3.1圆盘划分

简洁分析可知，假如旋转点在1区域，则沿180个方向照射后不会出现射线

分成两部分的状况；在n或in或IV区域，当吸取率出现最大值的时候射线也被

分成两部分，而不是像图1.2那样成为一部分，也解除；综合各因素可推断旋转

中心应当在V区域。因为放射-接收系统逆时针旋转，且图2,1.2中两部分红色区

域，间距缩小，说明旋转之后在垂直于放射-接收方向上二者的距离是缩短的，

从而确定旋转中心在短轴的上半侧，先运用Excel对附件二第1列数据作图（即

画出第一个方向上的扫描图像）如图4.1.3.2,发觉此时得到两部分图像，说明两

模板之间有一部分射线干脆被探测器接收，直到第14列数据两部分图像结合在

了一起，如图4.1.3.3。

结合托盘的几何特征，在垂直于对称轴方向上应当会出现最大的吸取率，利

用MATLAB求得出现最大吸取率的方向为第151个方向，在此方向结合增益率得

出的模板长度为79.9884,80mm进一步验证了结果；平行于对称轴方向上最大的

吸取率出现在沿对称轴的直线上，计算得出为第61个方向，此方向的模板长度

37.9988®30+8mm,也验证了结果。依据极近似水平方向为第61方向、最大吸

取率出现在235号探测器，极近似竖直方向为151方向、最大吸取率出现在223

号探测器，设旋转中心到竖直轴、水平轴的距离分别为x,y。

为求解还须要另一个方向的等量关系，选取椭圆和圆的一条外公切线的方向

,经计算得出经椭圆与圆公切线所在直线的X射线为第372号射线。则在所建立

坐标系内，切线方程为：

y=-0.8148^+41.8255

进而列出二元方程组:

y=%+(235—223)d

41.8225+0.8148行

>+--------.=——-=(372-235)a

V1+0.81482

解之得:

%=9.0396

y=5.7060

因此旋转中心在所建立的坐标系中的坐标为(-9.0396,5.7060)。

求180个旋转方向：

由于x射线的发生装置是连续旋转的，所以在512条射线中穿过模板的射线

条数应当是连续变更的，所以用MATLAB编写程序，计算每次旋穿过模板的x射

线条数，并画出图像如图4.1.3.4。

图4.1.3.4穿过模板射线条数随旋转次数变更曲线

图中横坐标为旋转次数，纵坐标为穿过模板的射线条数。从图中可以看出附

表2中的数据是依据射线放射装置旋转的依次依次给出的，而且可以看出，180

次旋转后装置共旋转了180度，每次旋转的角度近乎相等。

考虑到在前50次旋转中，穿过整个装置的X射线条数与穿过椭圆的条数相

等。设穿过椭圆的最边缘的射线到椭圆中心的距离为R,取椭圆中心为原点，模

板的对称中心为x轴建立平面直角坐标系，由于射线可以近似看成与椭圆和圆

x2+y2=R2都相切，可以得到以下方程组：

(x2+y2=R2

[y=kx+m

2z2

x乙y

---1---=

152402

y=kx+m

化简以上方程组可得:

利用MATLAB编程得到图像如图4.1.3.5所示

旋转角度与旋转次数的关系

90

8o

7o

.S6o

S5O

40-.•,

30

20---------1---------1---------'---------'---------1---------1---------'---------'---------1---------1

05101520253035404550

旋转次数

图4.1.3.5前50次旋转角度变更曲线

用同样方式得出最终面25组数据的图像为:

旋转角度与旋转次数的关系

30

155160165170175180

旋转次数

图4.1.3.6后25次旋转角度变更曲线

从图像中可以看出每次旋转角度近似为1度。

为了更加精确地计算旋转角度以便据此得到估算其他角度地依据，本文设计

了另外一种算法。从附表二中可以看出，前15组数据中穿过椭圆的射线与穿过

圆的射线没有交叉，所以数据中换算出来的最大吸取距离就是近似穿过椭圆中心

的射线被椭圆所截的距离，弦长公式为：

(x2y2

J152十402

(y=kx+m

且这个距离由直线的斜率k唯一确定，已知k、m时可解得:

1200V/C2+1

1=7152/c2-m2+402

152/c2+402

所以用MATLAB编程计算得出了较为精确的前15次旋转角度的结果如下:

29.649331.002731.557232.647633.679534.6478

35.648336.649037.647538.648639.648840.6491

41.648742.649743.6492

可以看出，这些结果的线性相关性特殊好，据此利用Exce1进行拟合，得到

如下图表:

旋转角度与旋转次数的关系

50

40

30

20

10

51015

旋转次数

图4.1.3.7旋转角度与旋转次数拟合公式

从图表中可以得出拟合公式为：y=0.9938+28.718,相关系数为：0.9996

依据拟合公式利用MATLAB编写程序计算得出全部方向角度。

4.2问题二模型的建立与求解

4.2.1建立模型

通过分析问题与附件数据可以发觉，题目所给数据与介质相关性质的二维分

布具有对应关系，X射线将介质一条线上的性质投影为一点。据此，本文利用

radon变换与radon逆变换进行运算，以通过投影后的讯号重建原始未知介质相

关性质的二维分布。该变换的定义为:令密度函数f(X)=f(x,y)是一个的定义域为心

的紧致台(compactsupport)o令R为radon变换的运算子，则Rf(x,y)是一个定义

在R2空间中的直线Lo

R/(n=J/(x)k/xi

其基本思想为：radon变换可以理解为图像在空间的投影，空间的每一点对应一

条直线，而radon变换是图像像素点在每一条直线上的积分。因此，图像中高灰

度值的直线会在空间形成亮点，而低灰度值的线段在空间形成暗点。对直线的检

测转化为在变换区域对亮点、暗点的检测。

Radon变换是一幅图像在一个特定的角度下的径向线方向的投影，一幅图像

的radon变换是每一个像素radon变换的集合。对于MATLAB中语句R=radon(l,

theta),假如theta是一个标量，R则是一个包含在theta的列向量。假如theta

是一个向量，R则是一个矩阵，矩阵的每一列是对应其中一个theta的radon变

换。

而radon变换的逆运算，就可以将CT系统对于每始终线上的X射线吸取率数

据还原回未知介质的物理性质。radon反变换的公式是:

f(x,y)=/g(xcos0+ysin0,0)d0

Jo

该反变换操作比较简洁，思路清晰，可以借助数学运算软件计算。

4.2.2模型求解

在MATLAB上将附件2和附件3的数据导入，利用radon逆变换将附件中数

据重建未知介质信息，得到图像如图4.2.2.1。

50

100

150

200

250

50100150200250100200300100200300

100

200

300

100200300100200300

图4.2.2.1radon逆变换重建介质图像(1)

为使逆拉东变换后得到的图像与原图像大小相等，在逆radon变换公式中取

362条射线，每两条射线相距0.2770,所以得到的图形长宽也为100,由第一问

求解得知，X射线从约29度位置起先旋转，为消退重建模型与真实模型的旋转

角度差异，将从29度逆radon变换得到的数据导入Excel表格可得变换后图像

坐标水平方向平移了31X0.2770=8.587,竖直方向平移了22X0.2770=

6.094,所以将原图中的坐标按上述数值平移就得到了变换后的坐标。

从模板的逆radon变换产生的矩阵中可以发觉，模板中全部点对应的灰度都

近似为0.5,又由于模板的吸取率为1,所以相对比例近似为2o据此利用MATLAB

编程可以算出图中对应十个点的吸取率如下表所示：

表4.2.2.1所求十点的吸取率(1)

序号12345

吸取率0.01850.9701-0.00011.17031.0287

序号678910

吸取率1.44251.2815-0.0074-0.00250.0194

然后依据逆radon变换的结果，将得的数据导入Excel表格，在不同范围内的数

值填充成不同颜色对不同吸取率的部分进行色块填充，得到结果如图4.2.2.2

所示。再利用Excel表格中找寻各椭圆定点的坐标，由于拉东变换中x射线之间

的距离都是0.2770,所以表格中两组数据在实际物体上的距离也是0.2770。由

100+0.2770《362得知表格取362组数据时，总长度与原图基本相同，此时可

以长度为基础，计算坐标变换公式。

依据前文坐标变换的逆变换，用MATLAB编写程序运算各色块（由运算可得

各色块均为椭圆形）的坐标及长短轴相应数据。将问题二中数据进行radon逆变

换后的图像最低点在Excel行数和列数，将行数和列数乘以倍率0.2770即距离

图像边界的距离，由于第一个图距离两边界的距离已知，可以得到平移的方向和

距离，具体结果如表4.2.2.2所示。

图422.2不同区域吸取率关系图（D

以每个椭圆中心的吸取率代表整个椭圆的吸取率则A,B,C,D,E,F的吸取率

分别为：0、1.1870、1.2914、0、0.9877、1.0632o

表422.2不同吸取率介质在托盘中的位置、几何形态与吸取率

椭圆编号ABCDEF

x063.017546.259054.707540.580551.245048.3365

yo33.101575.067071.189030.331553.184055.5385

halfl11.19205.25577.08529.645422.68871.8478

half26.11913.465312.88726.384539.53471.9439

吸取率01.18701.291400.98771.0632

其中X0、yO、halfl、half2分别表示椭圆中心横纵坐标和两个半轴。数据均以

左下角的点为坐标原点建立坐标系求得。

4.3问题三模型的建立与求解

4.3.1建立模型

与其次问类似，仍利用radon变换的思想建立模型，将空间每一条射线所投

影的点还原回一条直线，将数据合并后重建未知介质的几何性质与吸取率。

4.3.2模型求解

将附件5的数据导入MATLAB,利用radon逆变换将附件中数据重建未知介质

信息，得到图像如图4.3.2.1。

100

200

300

100200300100200300

100100

200200

300300

400400

500500

100200300400500100200300400500

图4.3.2.1radon逆变换重建介质图像（2）

其中图三图四显示完全，图二图四经角度修正后得出。

为具体算出所要求十个点的吸取率，将数据带入MATLAB中运算后结果如下

表所示：

表4.3.2.1所求十点的吸取率（2）

序号12345

吸取率0.07362.53057.13270.03320.8015

序号678910

吸取率3.14875.57140.05308.13620.0272

依据radon逆变换的结果，将得的数据导入Excel表格，在不同范围内的数

值填充成不同颜色对不同吸取率的部分进行色块填充，得到结果如图4.3.2.2

所示。

图4.3.2.2不同区域吸取率关系图（2）

其中无色处吸取率为0,黄色吸取率为0到2,红色为2到4,浅蓝色为4

到6,绿色为6到8,紫色为8到9。

从表格中大致取出图案中心，调用程序得到图案中心在(51.2450,47.9210)

(若以椭圆中心为原点，则图案中心在(1.2450,-2.0790))旁边。

4.4问题四的分析与求解：

为了便于求相邻探测器之间的距离，考虑最好仍选择圆形模板，但是要适当

增大圆形模板的半径，从而减小偶然误差；为了利用180次旋转得到的投射图像

求出旋转中心的坐标，鉴于在原来的标定模板中椭圆和圆的内公切线的选择有较

大误差，新的模板中应尽量使得切线简洁取得,且做出投射图像后简洁找到极值,

并且图像应有确定的对称性；图像扫描后应尽量少地得出重复数据，因此两个图

像差别应较大；考虑到便利地识别投射位置以确定旋转中心，应至少设置两个模

板且相互隔开。基于以上分析，设计出如下新模板：

图4.虫1新设计标定模板

标定方法：首先借助于圆形模板很简洁求得相邻探测器间距，利用对称性,

更简洁求得射线水平、竖直和一条倾斜方向的位置，因此用和第1问相同的思路,

此模板相对来说更能精确地确定CT系统的参数。

5.结果的分析与检验

5.1问题一结果分析与检验

第一问求得的单元间距，旋转中心与旋转角度与所搜集资料的实际值相差不

大且符合现实认知，在运算旋转角度时，由于线性关系良好，可以印证旋转中心

与单元间距计算误差不大。

5.2问题二结果分析

其次问十个位置的吸取率计算，由于是干脆由题目所给数据计算得出，结果

较精确，所得结果与逆radon变换所得图形相对应。在求解具体坐标与几何关系

时，通过Excel表格数据干脆计算得出，可能存在误差，但由具体结果运用radon

变换检验后可以看出，基本符合题目所给数据。

6.模型的优缺点分析

第一问中求单元间距与解旋转中心时利用了切线的特殊性质，但切线的选择

不愿定精确，因为射线宽度远小于探测器宽度且所用射线不愿定恰好为标定模板

切线位置。从题目中可以得出，由于x射线之间得间距相等，不管x射线怎么旋

转，穿过托盘上圆的X射线条数应当是大致相等的。依据附件二可以发觉，前

13组穿过圆的x射线条数都是29条，由此可以粗略的计算出x射线之间的距离

为空=02758。

29

但是这种做法并不精确，因为在29条射线中最两边的x射线并不是与圆相切的,

为尽量减小误差，本文利用图4.1.1.1进行误差检验，则如下方程成立：

JR2_堤_JR2_吗_]=d

设模板的吸取率为P，第n个数据为汇,上述方程可以转化为:

JR2—02s渡—JR2_p2sMi—d

即

2,R2_p2sl-JR2-p2sMi-J-2_P2sMi=°

利用MATLAB编程带入多组数据求其平均值可得k=1.7725,进而求得d=

0.2770.从这里可以看出在圆两边的射线不管是否与圆完全相切，对结果的影响

不是很大，所以该建模方式可以运用。

7.模型的改进

简洁地说第2,3问依据逆radon矩阵求出坐标值存在确定的偶然误差，计

算量大。而且本文的模型没有考虑到噪声等其他因素的影响，因此输出图像模糊

有光晕。为了得到清晰的图像，可以进行频域滤波。首先二维傅里叶变换对为：

F(u,v)=③/(x,y)e圻(皿^dxdy

J00J8

OO/•00

/Fe⑼/M卬叫如小;

/00J—OO

引入傅里叶切片定理，其中3是频率重量：

OO

g(p,0)e，2m「dp

/00

00/800

y)6(xcos0+ysin0—p)ej2vwpdxdydp

00

00■

3(xcos0+t/sin9—p)ei21fwpdpdxdy

00J-00J-00

/•0000

f(x,y)emmg861Vsioa)壹的

00

=F(3COS9,3sin9)

这说明一个投影的一维傅里叶变换，是二维投影矩阵的二维傅里叶变换的一个切

片，执行换元操作后，引入窗函数计算积分计算式并滤波，从而得到一个相对较

好的结果。

参考文献：工业CT技术刘丰林

工业CT系统旋转中心定位方法探讨刘明进

附录

%程序1.1作附件1的数据分布图

axisequal

fori=l:256

forj=1:256

if(A(izj)>0)

plot(j,257-iz,r**)

holdon

end

end

end

%程序1.2作附件2的数据分布图

fori=l:512

forj=l:180

if(AS(izj)>0)

if(AS(izj)>100)

plot(j,513-iz

holdon

else

plot(jz513z'r*-')

holdon

end

end

end

end

holdoff

%程序1.3计算探测系统的增益率及探测器的平均距离

tic

d=ones(86,1);%统计圆模板对应探测器的平均个数

m=zeros(86,1);%统计圆模板对应探测器的最大吸取率平均值

fori=l:14

k=0;x=0;

forj=374:430

if(AS(jzi)>0)

k=k+l;

end

if(AS(j,i)>x)

x=AS(jzi);

end

end

d(i)=k;m(i)=x;

end

fori=109:180

k=0;x=0;

forj=45:110

if(AS(jzi)>0)

k=k+l;

end

if(AS(jzi)>x)

x=AS(j,i);

end

end

d(i-94)=k;m(i-94)=x;

end

a=mean(d)

b=mean(m)

toe

%程序1.4求出竖直射线近似方向

d=zeros(180,1);k=0

forj=l:180

if(k<max(AS(:,j)))

d(j)=max(AS(:,j));

end

end

d

fori=l:180

if(k<d(i))

k=d(i)

i

end

end

宅程序1.5推出水平射线近似方向

d=zeros(180,1);u=l.7721;

forj=14:109

d(j)=max(AS(:zj));

d(j)

j

end

%程序1.6依据你和公式用matlab编写程序:

th=zeros(1,180);

fori=l:1:180

th(i)=0.9938*i+28.718;

end

th

计算得出全部的方向为:

1至15列

29.711830.705631.699432.693233.687034.680835.6746

36.668437.662238.656039.649840.643641.637442.631243.6250

16至30列

44.618845.612646.606447.600248.594049.587850.5816

51.575452.569253.563054.556855.550656.544457.538258.5320

31至45列

59.525860.519661.513462.507263.501064.494865.4886

66.482467.476268.470069.463870.457671.451472.445273.4390

46至60列

74.432875.426676.420477.414278.408079.4018803956

81.389482.383283.377084.370885.364686.358487.352288.3460

61至75列

89.339890.333691.327492.321293.315094.308895.3026

96.296497.290298.284099.2778100.2716101.2654102.2592103.2530

76至90列

104.2468105.2406106.2344107.2282108.2220109.2158110.2096

111.2034112.1972113.1910114.1848115.1786116.1724117.1662118.1600

91至105列

119.1538120.1476121.1414122.1352123.1290124.1228125.1166

126.1104127.1042128.0980129.0918130.0856131.0794132.0732133.0670

106至120列

134.0608135.0546136.0484137.0422138.0360139.0298140.0236

141.0174142.0112143.0050143.9988144.9926145.9864146.9802147.9740

121至135列

148.9678149.9616150.9554151.9492152.9430153.9368154.9306

155.9244156.9182157.9120158.9058159.8996160.8934161.8872162.8810

136至150列

163.87