基于直接线性变换算法的普通数码相机检校的应用研究
孔 建 黄建魏 沈 周
(西南交通大学 四川成都 610031 中铁十局 山东济南 520000) 摘要：本文采用直接线性变换（DLT ）算法，完成了普通数码相机检校的应用研究。通过编程实验，解算普通数码相机在不同焦距情况下内方位元素（00,x y ,f ）以及畸变参数（径向畸变系数1k ，2k 、偏心畸变系数1p ，2p ），同时对直接线性变换方法中l 初值的问题给出解决方案。提出了解决控制点布设在一个近似平面上解算l 系数初始值的方法，并且依据实验数据分析了在不同焦距下，相机内方位元素和光学畸变参数的变化情况。
关键字：直接线性变换；相机检校；径向畸变；偏心畸变
Abstract
In this paper, to complete a common application of digital camera calibration by using the direct linear transformation algorithm. This paper have solved different elements of interior orientation （00,x y ,f ）and distortion parameters （Radinal Distortion 1k ，
2k ,Decentering Distortion 1p ，2p ）of ordinary digital camera focal length by the programming experiments and meanwhile, put forward the solutions of the initial value problem in the direct linear transformation method. Proposed a solution in an approximate control points for solving plane initial value coefficient method, and analyzed the changes of the camera orientation elements and optical distortion parameters in the base of experimental data at different focal lengths.
1 概述
在数字摄影测量中，数字影像的获取，通常采用的是专业的摄影设备。这些专业设备的价格昂贵，对非专业部门是无法应用的。随着数码相机技术的发展与进步，普通数码相机在数字摄影测量领域中得到了广泛的应用，尤其是在近景数字摄影测量、无人机低空摄影测量的应用中，表现出了巨大的优势。普通数码相机不仅价格便宜，且操作方便，是专业摄影机不能比拟的。随着数码相机技术的
不断进步，其像幅、分辨率不断提高，数码相机相对于专业摄影机的优势更加凸显。
数码相机是非量测相机，应用到数字摄影测量中有其自身的缺陷，例如其像主点位置是未知的；相机镜头存在较大的光学畸变等。因此，为了保证影像的精度，获得可靠的测量结果，必须对数码相机的内方位元素及相关参数进行精确的检校确定。
非量测数码相机的内方位元素及相关参数检校的算法有很多，其中一种常用的方法是直接线性变换（DLT ）法。该方法解算过程中无需内、外方位元素的初始值，故特别适用于非量测相机的检校。其基本思想是在空间设置足够数量（大于等于6个）的控制点，用相机对控制点所在的目标进行摄影，最后将像点坐标和空间坐标代人直接线性变换的方程，引入光学畸变误差进行整体评差，即可解算出该相机的内外方位元素及相关参数。
2 非量测相机检校的内容和原理
2.1 相机检校内容
非量测像机检校的主要内容有：
（1） 像主点位置（00,x y ）与主距f 的测定；
（2） 光学畸变系数的测定，包括径向畸变系数（1k ，2k ）和偏心畸变系数
（1p ，2p ）的测定；
（3） 变焦后畸变差变化的测定； （4） 变焦后主距变化的测定； （5） 摄影机偏心常数的测定。
本文对摄影机检校研究的内容主要是内方位元素的检校和光学畸变参数的解算。
2.2非量测相机检校原理
2.2.1 内方位元素的确定
像片的内方位元素是恢复摄影时光束状态的要素；内方位元素确定了摄影中心S 与所摄像片P 相对位置关系，依据相对位置关系可以恢复摄影时光束的形状。借助内方位元素可以确定摄影中心与所摄像片间的位置关系，即恢复光束（光
线a S ，b S ，c S ）在摄影时候的形状。由此可见，确定内方位元素对摄影测量来说是重要的一步。理论上，像主点在底片（模拟影像）或者CCD （数字影像）的中心，但由于相机设计制造时的偏差，因此无论是量测相机还是非量测相机，像主点通常不与框标的连线或者像片中心重合，因此需要解算出像主点的偏差值（00,x y ）。如图2—1。若像主点在框标坐标系中的原点坐标为00,x y ，量测出来的像点坐标化算到以像主点为原点的像平面坐标系的坐标为（00,x x y y --）。
图2—1 内方位元素
2.2.2 非量测相机的光学畸变差 （1） 径向畸变差（Radial Distortion ）
径向畸变差使得构像点a 沿向径方向偏离其准确理想位置0a 。设在主距为f 的标准像片0P 上，物方点A 的标准位置为0a ，实际构像于点
a ，而且像方构像角'α与物方角α不等，如图2—3。由图可知，径向畸变差与物方点的入射角α有关，影像上不同点位的畸变差，因其α角度的不一样而不同。
图2—2 径向畸变
像点的径向畸变可以用如下奇次多项式来表达：
240122
4
012......
x k x k xr k xr y k y k yr k yr =+++=+++ （2—1）
其中，x 和y 是以像主点为原点的像片面坐标。
本文进行径向畸变改正时候只考虑了1k ，2k ，因为0k 与焦距具有相同的影响，因此不能与焦距校准同时发生[1]
。1k ，2k 的精度取决于控制点的精度。控制点像素坐标的量测精度应该高于0.2个像素才能用于相机校准[1]
。 （2） 偏心畸变差（Decentering Distortion ）
物镜系统的各单元透镜，因装配和震动偏离了轴线或者歪斜，从而引起像点偏离其理想位置。如图2—4，当轴线偏离其设计轴线12o o 或者旋转了一个角度，影响将会产生变形。
图2—3 引起偏心误差的原因
在近景数字摄影测量应用中，偏心畸变会随着焦距D 的变化而变化；而且不在焦距D 上的物点也存在偏心畸变差。一般情况下，偏心畸变比径向畸变小，但是这样的畸变仍然不能忽视。偏心畸变差的表达式如下：
22102002
2
20100[2()]2()()[2()]2()()
D D x p r x x p x x y y y p r y y p x x y y ?=+-+--?=+-+-- （2—2）
式中：
,D D x y ??——焦距为D 时的偏心畸变差分量；
f ——焦距为D 时的摄影时的主距；
12,p p ——偏心畸变系数；
r ——像点向径
]；
00,x y ——是像主点坐标；
图2—4即为径向畸变差r ?和偏心畸变差t ?对像点的共同影响。
图2—4 径向畸变和偏心畸变[1]
2.2.3 不垂直性误差d β和比例尺不一误差ds
当对数码相机进行量测时，像点坐标的改正数（,x y ??）中还包含了坐标轴不垂直性误差d β和比例尺不一误差ds 。如图2—5，像素坐标系c xy -是非直角坐标系，其两坐标轴之间的不垂直度为d β。以像主点o 为原点有两个坐标系，分别是直角坐标系o xy -和非直角坐标系o xy -。像主点的坐标为（0,0x y ）。像点'p 在像素坐标系中的坐标是（,x y ），在非直角坐标系o xy -中的坐标为（'21,om om ）。由于受d β和ds 的影响，此坐标包含线性误差。点p 是点'p 的理想位置，它在o xy -中的坐标是（,x y ）；其中2x on =，1y on =。
假设x 向无比例误差，而y 方向比例尺规划系数为1+ds 。此时x 向的主距为x f ，y 方向的主距y f 为：
1x
y f f ds
=
+ （2—3）
图2—5 坐标不垂直性误差和比例误差
比例尺误差ds 可以认为是所用的像素坐标系x 轴和y 轴的单位长度不一致及摄影材料的不均匀变形引起的。不正交性误差d β认为像素坐标x 轴和y 轴的不垂直性引起的。因此x ?和y ?的改正应为：
00()s i n ()x y y d y y y ds
β?≈-?≈- （2—4）
3 直接线性变换解法
非量测相机检校，常用的方法有光学实验室检校法，试验场检校法和在任检校法。其中后面两种方法皆是依据物方空间分布一群合理的高质量控制点，并依据单片空间后方交会解法或者多片空间后方交会解法，解求像片内外方位元素及各类光学畸变系数。本文所采用的方法是试验场检校法。
单片后方交会的解法主要有：基于共线条件方程式的单像空间后方交会、基于共面条件方程的空间后方交会解法、直接线性变化解法和基于角锥原理的空间后方交会解法。
基于共线条件方程式的单像空间后方交会解法缺点为解算前必须知道外方位元素的初始值，因而较适用于量测相机；基于共面条件方程的后方交会解法利用了像片间内在的几何关系，因而对控制点的要求比较少，但是这种解法精度稍低[2]。基于角锥原理的空间后方交会，能够很好的避免控制点的构想范围较小时候，外方位元素的不稳定，但是这种解法比较复杂。综合上述考虑，本文采用直接线性变（DLT ）的解法，该方法建立了像平面坐标和空间坐标之间的关系，计
算中无需内外方位元素的初始值，因此特别适用于非量测相机的检校。直接线性变换的模型如下：
1234
910115678
910110
1
1
l X l Y l Z l x l X l Y l Z l X l Y l Z l y l X l Y l Z ++++
=+++++++=+++ （3—1）
引入径向畸变差、偏心畸变差、坐标轴不垂直误差和比例尺误差并线性化后，其
误差方程为：
212349101110422201020025678910111042220100201
[(()()[2()]2()())1
[(()()2()()[2()])x y v l X l Y l Z l l Xx l Yx l Zx A k x x r A
k x x r p r x x p x x y y x v l X l Y l Z l l Xy l Yy l Zy A k y y r A
k y y r p x x y y p r y y y
=-+++++++-+
-++-+--+=-+++++++-+
-+--++-+ （3—2）
根据间接平差的原理，此误差方程式及其相应的法方程式可以写成：
1()T T V ML W L M M M W
-=-= （3—3）
这个运算过程是一个迭代的过程，以x f 相邻两次迭代运算的差值是否小于0.01mm 作为判断。 3.1 存在的问题及解决方案 3.1.1 存在的问题
若初始解算的控制点分布在一个平面上，则式(3—3)中的矩阵M 不是满秩矩阵[2]。直接线性变换解法不能用于控制点分布或者近似在一个平面的场合。因此当用非量测相机对近似在一个平面上的地物（建筑物）进行摄影时，使用DLT 解法就出现了问题。在使用十一个方程直接解算系数l 时，l 的初值无法求解。
在平差理论中，如果有多余观测值，在解求观测值或者待定值的平差值时，通常会在等式两边同时乘以系数矩阵的转秩矩阵。如果在解求l 的初值时，也在方程两边同时乘以系数矩阵M ，是否能够解决l 初值的问题呢？ 3.1.2 解决方案
首先，用n （n>=6）个控制点列出2n 个的方程式，即：
211
111
2n n M L W ???=- （3—4）
所以，
1
111
2112
11211
21()()T n n n n L M M
M W -??
???=- （3—5） 式（3—1）中有11个系数l ，解算系数l 的初值至少需要十一个方程，即需要6
个控制点。如果直接使用十一个方程解算，由于矩阵M （1,1）、M （3,3）、M （5,5）、M （7,7）都等于0；因此使用一般的矩阵求逆方法不能求解M 的逆矩阵。使用
（3—5）式，将矩阵T M 与M 相乘，能够避免矩阵M 中的上述四个元素为0，因而一般矩阵求逆方法也能使用。由于直接线性变化是一个迭代的过程，若所选择的6个控制点在一个近似的平面上，计算的l 的初值不稳定，但是将所有控制点都代入计算，在迭代的过程中对l 的值进行修正，因而最后解算的结果比较稳定。
3.2 数据分析
控制点的像点坐标和三维坐标如下：
控制点的像素和三维坐标 表3—1
点号 x （像素） y （像素）
X(m)
Y(m) Z(m) 20001 514.1 3040.4 497.5045 446.2537 295.6873 20002 496.4 2413.5 497.4540 446.2532 299.8948 20003 487.2 1758.5 497.4357 446.2199 304.2994 20004 505.2 1136.5 497.5883 446.2048 308.5020 20005 1287.8 2947.5 504.8410 454.1712 297.0586 20006 1191.2 1576.2 504.4446 454.3183 304.8016 20008 1867.9 2694.2 508.8083 457.6446 298.6852 20009 1836.7 1883.5 508.7136 457.6433 302.8787 20010 1845.8 1077.9 508.8224 457.6420 307.0718 20011 1622.1 906.5 506.0524 449.9470 309.2177 20012 2718.2 3052.2 512.9311 454.3017 296.3991 20013 2703 2322.5 512.9222 454.3034 300.4483 20014 2677.5 1561.7 512.8625 454.3071 304.6817 20016 3926.4 2322.1 520.3013 446.2188 300.1135 20017 3779.1 1682.9 519.4427 446.2059 304.2976 20020 5116.3 2943 527.9277 446.2150 295.9279 20021 5124.6 2297.3 528.0793 446.2067 300.1111 20022 5046.3 1653 527.6832 446.2027 304.3088 20023 5038.0 1011.0 527.7622 446.1945 308.5126
本文实验摄站的三维坐标的初始值为（500,500,300）；比较表3—1中的Y 坐标值，可知最大值与最小值之差=457.6446-446.2027=11.4491（m ）；摄站的Y
坐标与控制点的平均Y 坐标值的约为500-450=50（m ）。因此控制点相对于摄站并不能看作在一个平面上或者近似的一个平面上。控制点的分布及定义的独立坐标系O-XYZ 如图：
图3—1 控制点的分布及定义的独立坐标系O-XYZ
下面分四种种情况进行比较：
1）由于受目标点的限制，在布设控制点时候，只有四个平面，如表3—1中的Y 坐标所示。因此近似选择六个Y 坐标相差较大的点（2001、2002、2005、2010、2011、2020）用来解算l 的初始值，最后将所有的控制点代入方程，解算
l 值和S X ，S Y ，S Z ，?，ω，κ，0x ，0y ，x f ， y f ，1k ，2k ，1p ，2p ，d β，ds ；
2）用控制点（2003、2006、2010、2011、2020、2023）列十一个方程，解算l 的初值，最后将所有的控制点代入方程，解算l 值和s X ，S Y ，S Z ，?，ω，
κ，0x ，0y ，x f ， y f ，1k ，2k ，1p ，2p ，d β，ds ；
3）任意代入六个控制点（如2001、2002、2003、2004、2005、2006）解算
l 的初始值，最后将所有的控制点代入方程，解算l 的值和X s ，S Y ，S Z ，?，ω
，
X
Z
Y
O
κ，0x ，0y ，x f ， y f ，1k ，2k ，1p ，2p ，d β，ds ；
4）所选择控制点（2001、2002、2003、2020、2021、2022、2023）近似在一个平面上，解算l 的初始值，最后将所有的控制点代入方程，解算l 值和s X ，
S Y ，S Z ，?，ω，κ，0x ，0y ，x f ， y f ，1k ，2k ，1p ，2p ，d β，ds 。
上面四种情况所解算的外方位元素、内方位元素及,ds d β值如表3—2：
各个待求参数的最终解 表3—2
控制点不在一个平面上 十一个方程解算 任意六个
控制点
控制点在一 个 近似的平面 最大值与最小值之差
单
位 S
X
516.2523
516.2510
516.2559 516.2528
0.0013
m S Y 498.9267 498.9266 498.9262 498.9180 0.0082 m S
Z
301.5005 301.5005 301.4973 301.5000 0.0027
m
? 1.087465 1.088100 1.095114 1.085520 0.001946 rad ω 1.504859 1.505283 1.502772 1.504463 0.001691 rad κ
1.107982 1.108615 1.115540 1.106022 0.009518 rad
0x
16.88 19.58 0.67 14.71 18.91 像素 0y 12.95 14.79 9.79
10.73 5.00
像素 x f
8051.31 8051.27 8049.73 8048.50 2.81 像素 y f
8050.47
8050.15
8050.44 8047.79
2.69
像素
ds 8.187E-05 6.843E-05 1.015E-04 9.938E-05 3.305E-05 m d β
1.043E-04 1.381E-04 -8.725E-05 8.837E-05
2.253E-04 m 1k 1.253E-09 1.243E-09 1.185E-09 1.135E-09 1.177E-10 m 2k
3.135E-17 3.300E-17 3.692E-17
4.601E-17 1.466E-17 m 1p -2.519E-09 -1.415E-08 6.490E-08 1.113E-08 7.905E-08 m 2p
3.109E-08 2.189E-08 6.503E-08 2.927E-08
4.315E-08 m
分析表3—2可知：
1）比较表3—2的前面六行的各列数据，发现外方位元素解算值非常的接近，其中摄站坐标中相差最大的是S Y 为0.0082m ；外方位角元素相差最大的是κ为0.009518（rad ）=0°32′43″，不到一个厘米，这是空间后方交会自身的缺点导致的；内方位元素中0x 、
0y 、x f 、y f 最大值与最小值的差值最大的是
0x =18.919（像素）=18.919?6.41?610-=116.2?610-（m ），后面的六个参数的值相差更小。
2）表3—2种各个参数变换比较大的主要是第四列（任意六个控制点），其原因主要是2001到2006这六个点集中分布在像片的左半边（如图3—1），没有均匀分布在像片的构想范围内。
3）如果采用式（3—5），即使用多于方程解算系数l 的初值，虽然解算出来系数l 的初值相差比较大，不够稳定，但是如果有足够多的控制点，即使起算控制点都集中分布在一个近似的平面上面，对最后内方位元素、外方位元素、ds 、
d β、1k 、2k 、1p 、2p 解算的结果影响不是很大，但是应当注意初始起算数据
的控制点最好均匀在整幅像片中。
4 数据分析
以往对于非量测相机精度的探讨和分析都是基于同一个焦距下利用非量测相机获取的像片进行研究。如果将非量测相机直接用在各项工程中，由于野外环境比较恶劣，如摄影时候光线较差，摄站点与目标间的距离受限，在对目标摄影过程中，目标在影像上很小等，如果还是使用同一个主距，在应用中势必受到限制。因此有必要研究非量测相机在不同的焦距下，其内方位元素和物镜构像畸变参数是否会发生改变。 4.1试验方案
为了探讨数码相机检校的精度，现在利用Canon EOS-1Ds Mark III 对西南交大测量试验场进行摄影，采用佳能24-75毫米的变焦镜头，摄影主距分别为24mm ，50mm ，和60mm 。每个主距摄影三张以上的像片。为了使试验数据稳定，采用6个均匀分布在像片上的控制点（2001、2006、2010、2016、2020、2023）作为起算数据，像点坐标的量测采用的是更加精确的LISA 进行量测，其标称精度为0.2个像素。
在表4—1—1中像片1808与1803、1804大致在同一个位置拍摄，但在拍摄1808时，首先将主距调到50mm 拍摄另外一组像片后在调回24mm ，这样能够检验当相机主距重新调整后在调回去，内方位元素是否会发生改变。像片1830
是同1826、1827、1828大致在同一个位置拍摄，但1830使用的是60mm主距。
解算参数值（主距24mm）表4—1—1
1808 1803 1804 单位
XS 508.1348 508.3799 508.3757 m
YS 499.3719 498.8445 498.9346 m
ZS 301.5944 301.5934 301.599 m
phi 1.2995 -0.8602 -0.741 rad
omega 1.5628 1.5281 1.5279 rad
kappa -1.2884 0.8733 0.7538 rad
x0 -22.33 -44.11 -34.43 像素
y0 29.99 27.55 19.89 像素
fx 3816.54 3797.67 3799.67 像素
fy 3814.98 3796.44 3799.02 像素
ds -1.14E-04 3.48E-04 4.13E-04
Dbeta 4.09E-04 3.24E-04 1.71E-04
k1 1.12E-08 8.44E-09 8.57E-09 m
k2 -5.64E-16 3.34E-16 1.55E-16 m
P1 1.99E-07 7.60E-07 7.04E-07 m
P2 -3.54E-07 -5.24E-07 -4.34E-07 m
解算参数值（主距50mm）表4—1—2
1826 1827 1828 单位
XS 503.4188 503.4285 503.4041 m
YS 501.5148 501.4438 501.5247 m
ZS 301.6393 301.6225 301.6176 m
phi -1.3776 -1.4594 -1.5185 rad
omega 1.4322 1.4372 1.4449 rad
kappa 1.3771 1.462 1.522 rad
x0 -10.23 -12.21 -8.93 像素
y0 13.37 33.72 46.19 像素
fx 7614.63 7606.8 7624.53 像素
fy 7613.97 7608.42 7624.67 像素
ds 4.34E-04 3.36E-04 2.42E-04
Dbeta 8.62E-05 -2.14E-04 -1.81E-05
k1 -1.55E-09 -1.98E-09 -1.57E-09 m
k2 -1.14E-16 -6.00E-17 -1.04E-16 m
P1 2.67E-07 2.53E-07 1.44E-07 m
P2 -3.83E-07 -4.85E-08 -2.06E-07 m
解算参数值（主距60mm）表4—1—3
1830 1832 1833 单位
XS 503.4332 520.0471 520.0455 m
YS 501.424 501.8119 501.7926 m
ZS 301.624 301.5953 301.593 m
phi -1.4194 1.5009 1.5499 rad
omega 1.418 1.4437 1.4523 rad
kappa 1.4262 -1.5046 -1.5461 rad
x0 -18.95 -26.67 -15.86 像素
y0 55.78 69.58 81.8 像素
fx 9464.02 9490.79 9480.35 像素
fy 9461.72 9493.98 9480.7 像素
ds 1.84E-04 2.57E-04 3.61E-04
Dbeta 2.44E-04 -3.35E-04 -3.66E-05
k1 -2.61E-09 -2.66E-09 -2.55E-09 m
k2 -3.41E-18 8.45E-18 -7.33E-18 m
P1 2.32E-07 2.62E-07 1.90E-07 m
P2 -2.67E-07 2.75E-07 -7.79E-07 m
4.2 解算精度分析
4.2.1 摄站坐标精度分析
由于摄影时候是手持摄影，对每一张像片摄影的过程中，不可能保持同一个位置。而1826、1827、1828、1830是在间隔更短时间内拍摄的，拍摄瞬间位置基本上没变，因此这几张影像的摄站坐标（XS,YS,ZS）能够更好的反应解算的精度。其中XS、YS、ZS的最大值与最小值的偏差分别为3cm、10cm、2cm。其中YS的偏差最大是由于后方交会解法本身的缺陷所致，其解算的精度受到交会角大小的限制。
4.2.2 内方位元素精度分析
内方位元素的均值和方差（单位：像素）表4—2
24mm 50mm 60mm
均值中误差均值中误差均值中误差x-33.62 ±8.91 -12.36 ±1.35 -20.49 ±4.55
y25.81 ±4.30 15.06 ±13.53 69.05 ±10.63
f3804.63 ±8.46 7615.32 ±7.25 9478.39 ±11.02
x
f3803.48 ±8.20 7615.69 ±6.74 9478.80 ±13.24
y
在表4—2中内方位元素（0x ，0y ，x f ，y f ）最大中误差=±13.53
像素
=±0.0867mm 。比较不同主距的0x ，0y 可知，在不同的主距下，相机的内方位元素是变化的，因此使用对同一个相机，使用不同的主距，其具有不同的内方位元素。比较表4—1—1中1808与1803、1804内方位元素（0x ，0y ），可知重新调焦后对相机的内方位元素（0x ，0y ）的影响比较大。 4.2.3 相机参数变化
各项畸变参数均值和方差 表4—3
24mm 50mm
60mm
均值 中误差
均值 中误差
均值 中误差
ds 2.16E-04 ±2.35E-04 2.25E-04 ±7.84E-05 2.67E-04 ±7.27E-05 Dbeta 3.01E-04 ±9.86E-05 2.00E-04 ±1.24E-04 -4.25E-05 ±2.36E-04 k1 9.40E-09 ±1.27E-09 5.34E-09 ±1.98E-10 -2.61E-09 ±4.58E-11 k2 -2.51E-17 ±3.88E-16
1.81E-16 ±
2.35E-17
-7.62E-19 ±6.71E-18
P1 5.54E-07 ±2.52E-07 4.03E-07 ±5.50E-08 2.28E-07 ±2.96E-08 P2
-4.37E-07
±6.94E-08
-1.84E-07
±1.37E-07
-2.57E-07
±4.30E-07
（a ） ds 、Dbeta 随主距变化图
（b ） k1随主距变化图
（c ） P1、P2随焦距变化图 图4—1 各项畸变参数随焦距变化图
图4—1的（a ）、（b ）、（c ）三幅图可以看出，随着焦距的变化ds ，d β，1k ，
1p ，2p 的值也是变化的。ds ，2p 随着焦距的增加而变大；d β，1k ，1p 随着主距的变化而减小；由于2k 参数的值太小，在主距发生改变时候，可以看作是不变的。
从上面可以看出，无论是内方位元素（0x ，0y ）还是畸变参数，都将会随着主距的改变而改变。 5 结论与展望 5.1 结论
（1）本文采用直接线性变换（DLT ）算法，解算数码相机的内方位元素及相关参数的方法是可行的，通过实验证明计算结果是正确的；
（2）物方控制点在一个近似平面的情况下，采用间接平差的方法，亦能够得到精度较高的内方位元素和相关参数的值；
（3）该方法可用于定焦、变焦等多种非量测相机的检校；
（4）变焦镜头在变焦后对内方位元素影响较大，因此，建议尽可能使用定焦镜
头，如不得不使用变焦镜头时，应对镜头常用焦距段的内方位元素进行标定。5.2 展望
随着近景摄影测量应用的不增多和低空航空摄影测量技术的成熟，需要寻求简便可靠的数码相机的检校方法。数码相机检校可以从如下几个方面开展研究工作：
（1）为了给野外近景摄影测量工作带来简便，相机的检校应该力求无需控制场；
（2）相机检校的精度取决于数码相机畸变差改正的好坏，因此，为了提高检校精度，有必要探讨更好的畸变模型；
（3）为了提高控制点像素坐标的量测精度，应能够实现控制点像素坐标的自动获取；
（4）探讨将多种相机检校方法结合，如光束法，基于共线条件方程方法，基于角锥体原理的方法等，取长补短。
主要参考文献
[1] 冯文灏. 近景摄影测量——物体外形与运动状态的摄影法测定. 武汉：武汉大学出版社
[2] LPS帮助文档LPS_PM.pdf
1 武汉大学测绘学院测量平差学科组. 误差理论与测量平差基础. 武汉大学出版社，2003
2 中华人民共和国行业标准《近景摄影测量规范》，2008
3 中华人民共和国行业标准《低空数字航空摄影规范》，2010
4 冯文灏. 近景摄影测量——物体外形与运动状态的摄影法测定. 武汉：武汉大学出版社
5 王佩军、徐亚明. 摄影测量学. 武汉大学出版社，2005
6 潘正风，程效军，成枢，王腾军，杨正尧.数字测量原理与方法.武汉大学出版社，2004
7 吴剑，王广志，丁海曙，骆文博. 三维测量系统中线性CCD 相机的直接线性变换.清华大学学报，2004
8 葛宝臻，李晓杰，邱实.基于共面点直接线性变换的摄像机畸变纠正.中国激光，2010
9 刘兴库，李宝全，张锦生. 近景摄影测量空间后方交会新解法. 铁路航测，1992
10 郑波，陆珏，王卫安.非量测数码相机内方位元素及畸变差检校研究. 山西建筑，2010.5
11 张祖勋. 数字摄影测量的发展与展望.测绘学报，2003
12王建荣,杨俊峰,胡莘,李晶.空间后方交会在航天相机检定中的应用，测绘学报，2002
13王艳丽.近景摄影测量中数码相机检校的探讨.河南科技，2010
燕山大学 课程设计说明书 题目：双线性变换法设计数字低通滤波器 学院（系）：电气工程学院 年级专业：检测 学号： 学生姓名：沫沫 指导教师：王娜 教师职称：讲师
电气工程学院《课程设计》任务书 课程名称：数字信号处理课程设计 说明：1、此表一式四份，系、指导教师、学生各一份，报送院教务科一份。 2、学生那份任务书要求装订到课程设计报告前面。 电气工程学院教务科
目录 第一章绪论 (4) 第二章用双线性变换法设计低通滤波器的主要原理 (5) 2.1 双线性变换法 (5) 2.1.1 双线性变换法的基本原理 (5) 2.1.2 转换关系分析 (6) 2.2 巴特沃斯低通滤波器原理 (9) 第三章用双线性变换法设计低通滤波器步骤 (11) 第四章 MATLAB程序 (11) 第五章程序中命令介绍 (13) 第六章运行结果及波形 (14) 第七章结果分析 (16) 第八章心得体会 (17) 参考文献 (18)
第一章 绪论 数字滤波器是对数字信号实现滤波的线性时不变系统。数字滤波实质上是一种运算过程，实现对信号的运算处理。输入数字信号（数字序列）通过特定的运算转变为输出的数字序列，因此，数字滤波器本质上是一个完成特定运算的数字计算过程，也可以理解为是一台计算机。描述离散系统输出与输入关系的卷积和差分方程只是给数字信号滤波器提供运算规则，使其按照这个规则完成对输入数 据的处理。时域离散系统的频域特性:)()()(ωωωj j j e H e X e Y =,其中）（ω j e Y 、 )(ωj e X 分别是数字滤波器的输出序列和输入序列的频域特性（或称为频谱特 性）,)(ωj e H 是数字滤波器的单位取样响应的频谱，又称为数字滤波器的频域响 应。输入序列的频谱)(ωj e X 经过滤波后)()(ω ωj j e H e X ,因此，只要按照输入信号频谱的特点和处理信号的目的， 适当选择)(ωj e H ，使得滤波后的) ()(ω ωj j e H e X 满足设计的要求，这就是数字滤波器的滤波原理。 数字滤波器根据其冲激响应函数的时域特性，可分为两种，即无限长冲激响应(IIR)数字滤波器和有限长冲激响应(FIR)数字滤波器。IIR 数字滤波器的特征是，具有无限持续时间冲激响应，需要用递归模型来实现，其差分方程为： ∑∑==-+-=N i i N i i i n y b i n x a n y 1 )()()( 系统函数为： 计IIR 滤波器的任务就是寻求一个物理上可实现的系统函数H(z)，使其频率 响应H(z)满足所希望得到的频域指标，即符合给定的通带截止频率、阻带截止频率、通带衰减系数和阻带衰减系数。 ∑ ∑=-=-+= N k k k M r r r Z a Z b z H 1 01)(
一、 填空题（每题2分，共10题） 1、 1、 对模拟信号（一维信号，是时间的函数）进行采样后，就是 信号，再 进行幅度量化后就是 信号。 2、 2、 )()]([ω j e X n x FT =，用)(n x 求出)](Re[ω j e X 对应的序列 为 。 ３、序列)(n x 的N 点DFT 是)(n x 的Z 变换在 的N 点等间隔采样。 ４、)()(5241n R x n R x ==，只有当循环卷积长度L 时，二者的循环卷积等于线性卷积。 ５、用来计算N ＝16点DFT ，直接计算需要_________ 次复乘法，采用基2FFT 算法，需要________ 次复乘法，运算效率为__ _ 。 ６、FFT 利用 来减少运算量。 ７、数字信号处理的三种基本运算是： 。 ８、FIR 滤波器的单位取样响应)(n h 是圆周偶对称的，N=6， 3)3()2(2 )4()1(5.1)5()0(======h h h h h h ，其幅度特性有什么特性？ ，相位有何特性？ 。 ９、数字滤波网络系统函数为 ∑=--= N K k k z a z H 111)(，该网络中共有 条反馈支路。 １０、用脉冲响应不变法将)(s H a 转换为)(Z H ，若)(s H a 只有单极点k s ，则系统)(Z H 稳定的条件是 （取s T 1.0=）。 二、 选择题（每题3分，共6题） 1、 1、 )6 3()(π-=n j e n x ，该序列是 。 A.非周期序列 B.周期 6π = N C.周期π6=N D. 周期π2=N 2、 2、 序列 )1()(---=n u a n x n ，则)(Z X 的收敛域为 。 A. a Z < B. a Z ≤ C. a Z > D. a Z ≥ 3、 3、 对)70()(≤≤n n x 和)190()(≤≤n n y 分别作20点DFT ，得)(k X 和)(k Y ， 19,1,0),()()(Λ=?=k k Y k X k F ，19,1,0)],([)(Λ==n k F IDFT n f ， n 在 范围内时，)(n f 是)(n x 和)(n y 的线性卷积。 A.70≤≤n B.197≤≤n C.1912≤≤n D.190≤≤n 4、 4、 )()(101n R n x =，) ()(72n R n x =，用DFT 计算二者的线性卷积，为使计算量尽可 能的少，应使DFT 的长度N 满足 。 A.16>N B.16=N C.16 线性变换思想在中学数学中的应用 摘要：本文首先给出了线性变换的定义以及中学数学中涉及到的几种特殊的线性变换，包括其表达式及特征等。然后介绍了这几种线性变换在中学几何中的意义, 它是普通线性变换的一个自然推广,同时研究了线性变换在几何中的应用。最后,给出了具体实例说明了利用线性变换解决中学中平面几何题的方法以及线性变换思想在中学数学中的影响。 关键词：线性变换中学数学几何应用 随着社会的进步和时代的发展，针对我国中学数学课程现状，制定和实施新 的课程标准势在必行。2003年颁布了《普通高中数学课程标准(实验)》(以下 简称《标准》)。由参考文献[1]、[2]、[3]、[4]可知： 《标准》规定的课程与以往的课程相比，内容上发生很大的变化，尤其在选修系列中，增加了矩阵与变换、数列与差分、初等数论初步、优选法与试验设计初步、统筹法与图论、风险与决策、开关电路与布尔代数等内容，矩阵与变换是选修系列的内容。 矩阵是代数学的基本内容之一，变换是几何中的基本内容之一。对于中学数学教材改革来说，认真研究怎样把应用广泛的矩阵内容融入代数教材，以及如何进一步用变换的观念来处理几何教材，最终用矩阵来表示线性变换可以更有效地学习和运用这部分知识。中学数学引入矩阵初步知识，主要是为表达数据提供新的工具。矩阵作为研究图形（向量）变换的基本工具，有着广泛的应用，许多数学模型都可以用矩阵来表示。由矩阵建立的线性变换就是 平面上的坐标变换，其中，矩阵起着“对应法则”的作用，用二阶矩阵 a b c d ?? ?? ?? 确定的变换, 就是构造映射,使平面上的点(向量) x y ?? ?? ?? 变成(对应)点(向量)1 1 x y ?? ?? ?? = a b c d ?? ?? ?? x y ?? ?? ?? ,这个映射 吉林建筑大学 电气与电子信息工程学院 数字信号处理课程设计报告 设计题目：IIR数字滤波器的设计 专业班级：信工111 学生姓名：马超 学号：10211115 指导教师：高晓红王超 设计时间：2014.01.06－2014.01.10 目录 一、设计目的 (1) 二、设计内容 (1) 三、设计原理 (1) 3.1 数字低通滤波器的设计原理 (1) 3.2 变换方法的原理 (2) 四、设计步骤 (8) 五、数字低通滤波器 MATLAB 编程及幅频特性曲线 (10) 5.1 MATLAB语言编程 (10) 5.2 幅频特性曲线 (12) 六、总结 (13) 七、参考文献 (13) 一、设计目的 课程设计是理论学习的延伸，是掌握所学知识的一种重要手段，对于贯彻理论联系实际、提高学习质量、塑造自身能力等于有特殊作用。本次课程设计一方面通过MATLAB仿真设计内容，使我们加深对理论知识的理解，同时增强其逻辑思维能力，另一方面对课堂所学理论知识作一个总结和补充。 二、设计内容 已知通带截止频率f p =0.2kHz,通带最大衰减α p =1dB,阻带截止频率 f s =0.3kHz，阻带最小衰减α s =25dB，T=1ms，按照以上技术要求，用脉冲响应不 变法和双线性变换法设计巴特沃斯数字低通滤波器，并观察所设计数字滤波器的幅频特性曲线。 三、设计原理 3.1 数字低通滤波器的设计原理 滤波器是自动控制、信号处理和通信领域的重要组成部分,广泛地应用于各种系统中。数字滤波器是指输入、输出均为数字信号，通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。与模拟滤波器相比具有很多突出的优点，例如它可以满足滤波器对幅度和相位特性的严格要求，可以避免模拟滤波器所无法克服的电压漂移和噪声问题。 设计数字滤波器，首先要按一定规则将给出的数字滤波器的技术指标转换为模拟低通滤波器的技术指标。根据转换后的技术指标设计模拟低通滤波器G(s)，再按一定的规则将G(s)转换成H(z)。若是高通、带通或带阻数字滤波器则将它们的技术指标先转化为低通模拟滤波器的技术指标，然后设计出低通G(s)，再将G(s)转换为所需的H(z)将系统函数Ha(s)从s平面转换到z平面的方法有多种，但工程上常用的是脉冲响应不变法和双线性变换法。这两种方法有着各自的优缺点，在设计巴特沃斯数字低通滤波器时需要根据相应要求或想要达成的效果从中选择。 浅谈高中数学线性变换的解题技巧 在新课改之后，要求高中生不仅要学会灵活运用学科基础知识解决问题，还要利用课余时间学习自身兴趣的知识点，使得每个人都能得到全面发展和锻炼。高中线性变换虽然作为选修章节，但是其所蕴含的内容是衔接高中与大学的关键点，掌握线性变换的基础知识也就是提前了解和学习了大学所要接触的高等数学知识模块，即矩阵问题。因此，笔者立足于高中选修的重要知识点——线性变换，先阐述其概念及性质，然后来探究如何巧妙解决高中数学中线性变换的难题，从而为初等数学过渡到高等数学做提前的准备。 标签：数学线性变换解题技巧 一、高中数学线性变换的概述 1.线性变换的概念 线性变换一般是指，在构建的xOy坐标系内，存在至少一个点或多个点的集合A与另一个相对应的至少一个或多个点的集合B两者之间按照一定规则可以相互变换，且不同的点与所转变后的点不相同，即在平面直角坐标系中，把形如进行几何变换，这就叫做线性变换。 2.线性变换的基本性质 线性变换具有三个基本性质，第一个性质是任何向量乘于零都为零，数学表达式为：T（0）=0；第二个性质是任何向量乘于任何一个负向量等于两个向量相乘的负数，数学表达式为：T（-a）=-T（a）；第三个性质是线性变换满足乘法交换律、结合律，即，其中A是一般矩阵，是平面直角坐标系内任意的两个向量，是任意实数。 二、高中数学线性变换的解题技巧 1.数形结合 例1：在平面直角坐标系xOy中，已知平面区域A={（x，y）|x + y≤1，且x≥0，y≥0}，求平面区域B={（x + y，x - y）|（x，y）∈A}的面積。 解析：本题考察的是线性变换结合不等式的应用难点，解决该问题首先要分析题干信息，根据题目给出的信息列出平面区域A的不等式条件。由于本题平面区域B存在与平面区域A相重合的未知数，因此要假设两个新的未知数替代B的条件，再将新的未知数条件代入A中就能很快确定B的向量表示，最后快速建立平面直角坐标系画出平面区域B的图形就能的出其面积的大小。 设：未知数u=x+y，v=x-y 实验五 IIR 数字滤波器设计 一、实验目的 （1）熟悉用双线性变换法设计IIR 数字滤波器的原理与方法； （2）学会调用MATLAB 信号处理工具箱中滤波器设计函数设计各种IIR 数字滤波器，学会根据滤波需求确定滤波器指标参数。 （3）掌握IIR 数字滤波器的MATLAB 实现方法。 （4）通过观察滤波器输入输出信号的时域波形及其频谱，建立数字滤波的概念。 二、实验原理 设计IIR 数字滤波器一般采用间接法（脉冲响应不变法和双线性变换法），应用最广泛的是双线性变换法。基本设计过程是：①先将给定的数字滤波器的指标转换成过渡模拟滤波器的指标； ②设计过渡模拟滤波器；③将过渡模拟滤波器系统函数转换成数字滤波器的系统函数。MATLAB 信号处理工具箱中的各种IIR 数字滤波器设计函数都是采用双线性变换法。第六章介绍的滤波器设计函数butter 、cheby1 、cheby2 和ellip 可以分别被调用来直接设计巴特沃斯、切比雪夫1、切比雪夫2和椭圆模拟和数字滤波器。本实验要求读者调用如上函数直接设计IIR 数字滤波器。 本实验的数字滤波器的MATLAB 实现是指调用MATLAB 信号处理工具箱函数filter 对给定的输入信号x(n)进行滤波，得到滤波后的输出信号y(n ）。 三、实验内容及步骤 （1）调用信号产生函数mstg 产生由三路抑制载波调幅信号相加构成的复合信号st ，该函数还会自动绘图显示st 的时域波形和幅频特性曲线，如图5.1所示。容易看出，图5.1中三路调幅信号的载波频率分别为250Hz 、500Hz 、1000Hz 。 抑制载波单频调幅信号的数学表示式为 0001()cos(2)cos(2)[cos(2())cos(2())]2 c c c s t f t f t f f t f f t ππππ==-++ 其中，cos(2)c f t π称为载波，f c 为载波频率，0cos(2)f t π称为单频调制信号， DLT 直接线性变换解法程序介绍 一、程序综合介绍：DLT结算程序 程序功能介绍：应用6个已知点计算左右片l 系数；然后应用已经求得的l系数求解物方空间坐标系坐标 程序名：SuYGDLT 程序界面： 程序界面有四个按钮，分别为读取文件，左片l系数计算，右片系数计算，物放坐标解算程序界面有四个编辑框，分别用来输出文件信息，左片l系数、右片l系数、以及无妨坐标结果 截图如下 程序使用介绍： 必须先点击导入文件按钮，导入文件方可进行正确的计算，如果未导入文件就点击左片平差或右片平差或无妨坐标解算就会弹出如下对话框： 读取数据后点击其它按钮进行其它计算。 程序文件格式： 数据文件分为两部分，KnownPoint，UNKnownPoint，分别代表已知点信息和待求点信息当文件读取程序读到“KnownPoint”时开始读取已知点信息，已知点信息格式如下 GCP1,1214.0000,1032.0000,1046.5180,1071.6652,9.201742,-9.672384,-2.726064 分别代表点名、左片相片X坐标、左片相片y坐标、右片相片x坐标、右片相片y坐标物方坐标X、Y、Z; 当文件读取到“END KnownPoint”时结束已知坐标的读取 待求点信息类似：文件格式截图如下： 程序运行结果与评估： 本程序区1-10号点作为已知点计算l近似值11-20号点作为未知点解求其物方三维坐标； 程序运行结果与所给参考值相似，应该可以证明其运算是正确的，运行结果截图如下： 二、程序编程思想及相关代码 程序编程思想及相关函数： 本程序设计DLTCalculation类作为l系数结算主程序，其成员变量及成员函数与作用介绍如下： CSuLMatrix LL;//左片L系数矩阵 CSuLMatrix RL;//右片L系数矩阵 int m_iKnownPointCount;//已知点个数 CControlPoint *m_pKnownPoint;//已知点 int m_iUnKnownPointCount;//未知点个数 CControlPoint *m_pUnKnownPoint;//未知点 public: CString LoadData(const CString& strFileName);//读取文件函数 int ifLoda;//判断是否导入数据 CString Datainfor;//文件信息存储 CString *SplitString(CString str,char split, int& iSubStrs); //分割函数 void LFormApproL(CSuLMatrix &LL);//计算左片L系数近似值 void RFormApproL(CSuLMatrix &RL);//计算右片L系数近似值 void FormLErrorEquations(CSuLMatrix LL,CMatrix &LM,CMatrix &LW);//组成左片系数矩阵和常数项矩阵 void LAdjust();//左片平差主函数 void FormRErrorEquations(CSuLMatrix RL,CMatrix &RM,CMatrix &RW);//组成右片系数矩阵和常数项矩阵 void RAdjust();//右片平差主函数 void Output(const CString& strFileName);//输出结果主程序 第２２卷第７期２００３年７月 大学 ＣＯＩ，ＬＥＧＥ 物理 ＰＨＹＳＩＣＳ Ｖ０ｌ２２Ｎｏ７ Ｊｕｌｙ２００３线性常微分方程的保线性变换及其应用 柬仁贵１，束萱２，李珍１ （１．东北师范大学物理系，吉林长春１３００２４；２清华大学经管学院．北京１０００８４） 摘要：研究了线性常微分方程的保线性变换，得到任意两个二阶线性常微分方程等价的条件，并用于求解一类二阶线性变系数齐次常微分方程对数学物理方法教学中怎样通过适当的变换把给定的二阶线性变系数齐次常微分方程化为可解的方程给出了合理解释 关键词：二阶线性常微分方程；保线性变换 中圈分类号：ｏ４１ｌ文献标识码：Ａ文章编号：１０００．０７１２（２００３）０７—００１１—０５ ｌ引言 二阶线性常微分方程在物理学及科学技术中有广泛的应用“。２ｏ根据线性常微分方程的一般理论，求任何线性非齐次常微分方程的解都归结为求相应齐次常微分方程的基本解组”…然而，即使对二阶线性变系数常微分方程，至今也没有求出其基本解组的一般方法在通常的高等数学及数理方法教科书中都给出了二阶线性常系数常微分方程及一些特殊函数方程求解方法，如果能够通过适当的变换把一个给定的二阶线性变系数常微分方程化为这些可解的方程，则很容易得出该二阶线性变系数常微分方程的通解问题在于如何知道这个方程能否化为可解的方程并通过怎样的变换才能化为可解的方程 本文基于线性常微分方程的一般特性”Ｊ，研究了常见的保线性变换．得到联系两个二阶线性常微分方程的等价关系，进而导出一个给定的二阶线性变系数常微分方程可化为已知可解的方程的判别式并确定相应的变换这就扩大了二阶线性变系数常微分方程求解的途径，同时可对数学物理方法教材中一些常见的变换“探本求源”，给出合理的解释 ２常见的保线性变换 在常微分方程的一般理论中，线性性质尤为重要能保持微分方程的线性性质的变换统称为保线性变换下面讨论最常见的两种保线性变换 ２．１变换，（ｚ】．“（ｔ）ｅ“”的保线性性质 为确定起见．在以下讨论中规定一般的二阶线性变系数常微分方程的标准方程为 挈岬（ｚ）塞＋Ｑ（ｍ＝０（１）对方程（１）作变换得 ｙ（ｚ）＝“（ｆ）ｅ“”（２）即同时作未知函数的齐次线性变换和自变量的变换，经过通常的微分运算方程（１）可化为未知函数“关于自变量为￡的二阶线性常微分方程： ；（￡）２＋二［塞＋ｚ二（老）２＋Ｐ害］＋“［；（塞）２＋；（睾＋，盎）＋ｉ２（￡）２＋Ｑ］＝。 （３）可见变换（２）能保持微分方程的线性性质，是保线性变换 ２．２二阶线性齐次常微分方程对自变量求ｎ阶导数的线性不变性 为便于讨论，设一般的二阶线性变系数常微分方程为 ｎｏ（ｚ）ｙ”（ｚ）＋ｎｌ（。）Ｊ７（ｚ）＋口２（ｚ）√（Ｔ）＝０（４） 收稿日期：２００２０７１ｓ 基盒项目：东北师范大学“优师工程”资助项目 作者简介：束仁责（１９４３一），男．江苏镇距人，东北师范大学物理系剐教授，硕士，主要从事教学物理研究与教学，数学物理方法课主讲教师 第六章 线性变换 映射：,X Y ≠?≠?，如果有一个法则σ，它使得X 中每个元素α，在Y 中有唯一确定的元素β与之对应，则称σ为X 到Y 的一个映射，记作:X Y σ→，()σαβ=，β称为α在σ下的象，α称为β在σ下的原象。 注：()(),X στασατα=??∈=对。 变换：一个集合到自身的映射。 线性变换的定义与性质 定义 设V 是数域F 上的线性空间，σ是V 的一个变换，如果满足条件： （1）()()()βσασβασV,α,β+=+∈?； （2）()()k F,αV,k αk σασ?∈?∈=， 则称σ是V 上的线性变换或线性算子。 (1), (2)等价于条件：,,,k l F V αβ?∈∈ ()()()σk αl βk σαl σβ+=+。 例：设σ：n n R R →，定义为()c αασ=，c 为常数。-----数乘 变换或位似变换。 c =0-----零变换，记为o 。 c =1-----恒等变换，记为ε。 例：设σ是把平面上的向量绕坐标原点逆时针旋转θ角的变换 设()()(),,,T T x y x y ασα''==，则 cos sin sin cos x x y y x y θθ θθ'=-??'=+? 记cos sin sin cos A θθθ θ-?? =??? ? ，则()A σαα=是一个线性变换。 例：判断下列变换是否是线性变换 (1) ()()12323,,1,,T T a a a a a σ=; (2) ()()12323,,0,,T T a a a a a σ=; (3) ()()12312231,,2,,T T a a a a a a a a σ=-+; (4) ()()212312 3,,,,3T T a a a a a a σ=. 线性变换的基本性质 （1）()θθσ=； （2）()()ασασ-=-； （3）线性变换保持向量的线性组合关系不变，即若s s αk αk αk β+++=Λ2211，则1122s s βk αk αk ασσσσ=+++L ； 若θ=+++s s αk αk αk Λ2211，则θσσσ=+++s s αk αk αk Λ2211。 （4）线性变换将线性相关的向量组映成线性相关的向量组。 线性变换的运算 ()V L ----线性空间V 上所有线性变换的集合。 实验三 用双线性变换法设计 IIR 滤波器 一．实验目的 1．了解工程上两种最常用的变换方法：脉冲响应不变法和双线性变换法。 2．掌握双线性变换法设计IIR 滤波器的原理及具体设计方法，熟悉用双线性设计法设计低通、带通和高通IIR 数字滤波器的计算机程序。 3．观察用双线性变换法设计的滤波器的频域特性，并与脉冲响应不变法相比较，了解双线性变换法的特点。 4．熟悉用双线性变换法设计数字Butterworth 和Chebyshev 滤波器的全过程。 5．了解多项式乘积和多项式乘方运算的计算机编程方法。 二．实验原理与方法 从模拟滤波器设计IIR 数字滤波器具有四种方法：微分－差分变换法、脉冲响应不变法、双线性变换法、z 平面变换法。工程上常用的是其中的两种：脉冲响应不变法、双线性变换法。脉冲响应不变法需要经历如下基本步骤：由已知系统传输函数H(S)计算系统冲激响应h(t)；对h(t)等间隔采样得到h(n)=h(nT)；由h(n)获得数字滤波器的系统响应H(Z)。这种方法非常直观，其算法宗旨是保证所设计的IIR 滤波器的脉冲响应和模拟滤波器的脉冲响应在采样点上完全一致。而双线性变换法的设计准则是使数字滤波器的频率响应与参考模拟滤波器的频率响应相似。脉冲响应不变法一个重要的特点是频率坐标的变换是线性的（w = ΩT ），其确定是有频谱的周期延拓效应，存在频谱混叠的现象。为了克服脉冲响应不变法可能产生的频谱混叠，提出了双线性变换法，它依靠双线性变换式： s=1111--+-z z ，z=s s -+11，其中 s = s + j Ω，z = re jw 建立其S 平面和Z 平面的单值映射关系，数字域频率和模拟域频率的关系是： Ω = tg (w / 2)， w = 2arctg (Ω) （3－1） 由上面的关系式可知，当Ω →∞时，w 终止在折叠频率w = p 处，整个j Ω轴单值的对应于单位圆的一周。因此双线性变换法不同于脉冲响应不变法，不存在频谱混叠的问题。从式（3－1）还可以看出，两者的频率不是线性关系。这种非线性关系使得通带截至频率、过渡带的边缘频率的相对位置都发生了非线性畸变。这种频率的畸变可以通过预畸变来校正。用双线性变换法设计数字滤波器时，一般总是先将数字滤波器的个临界频率经过式（3－1）的频率预畸变，求得相应参考模拟滤波器的个临界频率，然后设计参考模拟滤波器的传递函数，最后通过双线性变换式求得数字滤波器的传递函数。这样通过双线性变换，正好将这些频率点映射到我们所需要的位置上。参考模拟滤波器的设计，可以按照一般模拟滤波器设计的方法，利用已经成熟的一整套计算公式和大量的归一化设计表格和曲线。这些公式、表格主要是用于归一化低通原型的。通过原型变换，可以完成实际的低通、带通和高通滤波器的设计。在用双线性变换法设计滤波器的过程中，我们也可以通过原型变换，直接求得归一化参考模拟滤波器原型参数，从而使得设计更加简化。理论课教材给出了IIR 低通、带通和高通滤波器设计双线性原型变换公式的总结，请参阅。在本实验中，我们只设计两种滤波器（Butterworth 和Chebyshev ）的设计，相应的这两种参考模拟原型滤波器的设计公式见理论课教材。 综上所述，以低通数字滤波器设计为例，可以将双线性变换法设计数字滤波器的步 骤归纳如下： 1．确定数字滤波器的性能指标。这些指标包括：通带、阻带临界频率f p , f s ；通带内的最大衰减a p ；阻带内的最小衰减a s ；采样周期T 。 目录 摘要 (1) Abstract (2) 引言 (3) 1 线性变换 (4) 1.1 线性变换的定义 (4) 1.1.1 线性变换的概念 (4) 1.1.2 线性变换的矩阵及矩阵表示 (4) 1.2 矩阵的相似对角化问题 (5) 1.2.1 相似对角化问题 (5) 1.2.2 矩阵的特征值与特征向量 (5) 2 线性变换的对角化 (7) 2.1 线性变换的对角化 (7) 2.1.1 线性对角化的提出 (7) 2.1.2 线性对角化的定义 (7) 2.2 线性变换的特征值与特征向量 (7) 2.2.1 线性变换的特征值与特征向量的概念 (7) 2.2.2 线性变换的特征多项式 (7) 2.3 线性变换对角化与矩阵对角化之间的联系 (8) 2.3.1 特征值与特征向量的联系 (8) 2.3.2 线性变换对角化与矩阵相似对角化之间的关系 (9) 2.3.3 线性变换可对角化的充要条件及推论 (9) 2.3.4 求线性变换对角化的方法和步骤 (10) 3 线性对角化问题的相关题目 (14) 总结 (16) 参考文献 (17) 致谢 (18) 摘要 线性变换是贯穿高等代数的重要内容之一，其研究价值不言而喻。本文尝试通过探讨矩阵对角化的知识点类比线性变换对角化的知识点，再通过矩阵的特征值与特征向量，以线性对角化问题为主要线索，着手研究线性变换特征值与特征向量的求解步骤以及线性对角化的基本条件，并且总结说明线性变换的对角化与矩阵对角化的联系，更进一步的，加深了解矩阵对角化与线性对角化的内容及要点。 关键词：线性变换的对角化问题；矩阵；特征值；特征向量 Linear transformation is an important part of higher algebra through its research value is self-evident. This paper attempts to explore the matrix diagonalization by knowledge points of analog linear transformation diagonalization knowledge, and through the eigenvalues and eigenvectors of the matrix, linear diagonalization problem as the main clue, started studying linear transformations eigenvalues and eigenvectors steps to solve the basic conditions and linear keratosis, and summary description of the linear transformation matrix diagonalization diagonalization with links to further deepen understanding of linear matrix diagonalization diagonalization content and points. Keywords: Changing existing diagonalization；Matrix；Eigenvalues；Eigenvectors 双线性变换法IIR 数字滤波器设计 一、实验目的 1、掌握用双线性变换法设计低通IIR 数字滤波器的基本原理和算法。 2、掌握用双线性变换法设计高通和带通IIR 数字滤波器的基本原理和算法。 3、进一步了解数字滤波器和模拟滤波器的频率响应特性。 二、实验原理与计算方法 １、双线性变换法设计IIR 低通数字滤波器的基本原理和算法 双线性变换法设计数字滤波器，采用了二次映射的方法，就是先将整个s 平面压缩到 s 1平面的一个T j T j π π~-的横形条带范围内，然后再将这个条带映射到z 平面上，就能建 立s 平面到z 平面的一一对应关系。对于低通数字滤波器，映射关系为 z z T z z T s ++-= +-=--11211211 (1) 其中T 为抽样周期。 用双线性变换法设计低通IIR 数字滤波器的基本步骤，首先根据设计要求确定相应的模拟滤波器的传递函数)(s H a ，再应用(1)式得数字滤波器的传递函数)(z H z z T s a s H z H ++-==112)()( (2) 通常可以给定的参数为：低通数字滤波器通带边界频率p p f πΩ21=、阻带边界频率 s s f πΩ21=和对应的通带衰减函数p α、阻带衰减函数s α。s 1平面中的模拟角频率1Ω与数字角频率ω的关系为线性关系T 1Ωω=，在计算模拟滤波器的阶数N 、极点s i 和传递函数)(s H a 之前，应作预畸变处理 2 2tan 22 tan 2 1T f T T T πΩΩ== (3) 模拟滤波器的阶数N 、极点s i 和传递函数)(s H a 的计算方法与冲激响应不变法相同，可以采用Butterworth 逼近或Chebyshev 逼近。 ２、双线性变换法设计IIR 高通、带通、带阻数字滤波器的基本原理和算法 由于双线性变换法获得的数字滤波器频率响应特性中不会出现混叠现象，因此可以适用于高通、带通和带阻滤波器的设计。IIR 数字滤波器的设计通常要借助于模拟低通滤波器的设计，由原型低通滤波器到其他形式（高通、带通、带阻）IIR 数字滤波器的频带变换有模拟频带变换法和数字频带变换法。 （1）模拟频带变换法 首先将给定的对数字滤波器(DF)的技术要求转换为一个低通模拟滤波器(AF)的技术要求，根据这种要求用某种逼近设计出原型的低通模拟滤波器(LP AF)，计算出模拟滤波器的阶数N 、极点s i 和传递函数)(s H a ，再按照双线性变换的变换关系，将模拟滤波器的传递函数)(s H a 转换为数字滤波器的传递函数)(z H 。 表8-1中列出了将给定的对数字滤波器(DF)的技术要求直接转换为对一个低通模拟滤波器(AF)的技术要求的频率预畸变校正关系和转换公式。 直接线性变换Matlab实现的程序源代码 function re=DLT(A,B) %imco为像方坐标，输入单位是像素 imco=A; %此处为控制点像方坐标，格式为2×n，单位：像素 %obco为物方坐标，输入单位是毫米 obco=B; %此处为控制点物方坐标，格式为n×3单位：毫米 imco_be=[];B=[];M=[]; for i=1:size(imco,2) imco_be=[imco_be;imco(:,i)]; end for i=1:size(imco,2) A1=[obco(i,:),1,0,0,0,0]; A2=[0,0,0,0,obco(i,:),1]; M=[M;A1;A2]; B1=obco(i,:).*imco_be(2*i-1); B2=obco(i,:).*imco_be(2*i); B=[B;B1;B2]; end M=[M,B]; N=M(1:11,:); L=N\(-imco_be(1:11,:)); X0=-((L(1)*L(9)+L(2)*L(10)+L(3)*L(11))/(L(9)*L(9)+L(10)*L(10)+L(11)*L(11))); Y0=-((L(5)*L(9)+L(6)*L(10)+L(7)*L(11))/(L(9)*L(9)+L(10)*L(10)+L(11)*L(11))); L=[L;0];M3=[];W=[]; for i=1:size(imco,2) xyz=obco(i,:); A=xyz(1)*L(9)+xyz(2)*L(10)+xyz(3)*L(11)+1; r2=(imco_be(2*i-1)-X0)*(imco_be(2*i-1)-X0)+(imco_be(2*i)-Y0)*(imco_be(2*i)-Y 0); M1=[A*(imco_be(2*i-1)-X0)*r2;A*(imco_be(2*i)-Y0)*r2]; M2=-[M(2*i-1:2*i,:),M1]/A; M3=[M3;M2]; W=[W;-[imco_be(2*i-1);imco_be(2*i)]/A]; end WP=M3'*W; NBBN=inv(M3'*M3); LP=-NBBN*WP; v=M3*LP+W; imco_be=imco_be+v; X0=-(LP(1)*LP(9)+LP(2)*LP(10)+LP(3)*LP(11))/(LP(9)*LP(9)+LP(10)*LP(10)+LP (11)*LP(11)); Y0=-(LP(5)*LP(9)+LP(6)*LP(10)+LP(7)*LP(11))/(LP(9)*LP(9)+LP(10)*LP(10)+LP (11)*LP(11)); 1 实验＿三_题目_用双线性变换法设计IIR数字滤波器第１6周星期_３_第６,7节 一.实验目得 （1)熟悉用双线性变换法设计IIR数字滤波器得原理与方法。 (２）掌握数字滤波器得计算机仿真方法。 (３)通过观察对实际心电图信号得滤波作用,获得数字滤波得感性知识。 二、实验内容、方法、设计程序及实验结果 (1)复习有关巴特沃斯模拟滤波器设计与用双线性变换法设计ＩIR数字滤 波器得内容,用双线性变换法设计数字滤波器系统函数。其中满足本实验 要求得数字滤波器系统函数为: () () ()()()2 1 2 1 2 1 6 1 2155 .0 9044 .0 1 3583 .0 0106 .1 1 7051 .0 2686 .1 1 1 0007378 .0 - - - - - - - + - + - + - + = z z z z z z z z H ???????(3、1) 式中: ?(3、２) 根据设计指标,调用MATＬＡＢ信号处理工具箱butｔｏrd与ｂutｔer,也可以得到。由公式(3、1)与（３、２)可见,滤波器由三个二阶滤波器、与级联而成,如图3-1所示。 % x=、 、、 k=1；%控制滤波循环变量 clｏｓe alｌ; %关闭全部绘图窗口 fｉgure(1);%创建绘图窗口 subplot（2,2,1)；%定位子图１ n=0:55; %横坐标 sｔem(n,x,'、'）; ％画出枝干图 ａxiｓ([0,56，-100,5０]）; %调整坐标 ｘlａbeｌ（'n'); %标注横坐标 ylaｂｅl(＇ｘ(ｎ)');%标注纵坐标 titｌｅ('心电图信号采集序列x（n)'）;%命名该子图 B=[0、090３６,2*0、09036,0、０9０36]；%H1 滤波器得分子系数矩阵 Ａ=[1、２686,-0、7０51]; %H１滤波器得分母系数矩阵 A1＝[１、0１0６，-0、３58３]；%H2滤波器得分母系数矩阵 A2=［０、9044,-0、21５5］; ％H３滤波器得分母系数矩阵 ｗhiｌe(k<=3) y=fiｌteｒ（B,A,x); %进行滤波 图3-1 滤波器得组成 三维线性变换 陈祥科 1、线性空间 (2) 1.1、线性空间的代数定义 (2) 1.2 线性空间的基和维度 (2) 2、线性变换 (2) 2.1、变换的定义 (2) 2.2、线性变换的定义 (2) 2.3线性变换的性质 (3) 2.4、线性变换下的坐标变换 (3) 2.5、线性变换的矩阵表示： (3) 3、三维图形的几何变换 (4) 3.1平移变换 (5) 3.2缩放变换 (5) 3.3绕坐标轴的旋转变换 (5) 3.4绕任意轴的旋转变换 (6) 4、三维线性变换的应用实例 (7) 4.1 三维图形变换理论 (7) 4.1.1 三维图形的几何变换 (7) 4.1.2 组合三维几何变换 (8) 4.1.3 围绕任意轴的旋转矩阵的推导 (9) 4.1.4 三维图形的轴侧投影变换 (9) 4.2 叉车稳定性试验的仿真 (10) 4.2.1 纵向稳定性试验的仿真 (10) 4.2.2 横向稳定性试验的仿真 (11) 4.3 结论 (12) 1、线性空间 1.1、 线性空间的代数定义 一个定义了加法与数乘运算，且对这些运算封闭，空间中任意向量都属于数域P ，并满足八条算律的集合为数域P 上的线性空间。 1.2 线性空间的基和维度 对于一个数域上的线性空间R ，由n 个属于R 的元素组成的一个线性无关组，如果R 中的任意一个元素都是这n 个元素的线性组合，那么这个线性空间的维度为n ，且这个线性无关组为R 的一组基。显然，三维空间的基有3个元素组成。三维线性空间的的两组基分别为(0,0,1)和(1,0,0)、(0,1,0)。 2、线性变换 2.1、变换的定义 变换是广义概念的函数，它是这样定义的，如果存在2个非空集合A 、B ，α是A 中的任意元素，如果在集合B 中必定有一个元素β与集合A 中的α元素对应，则称这个对应关系是集合A 到集合B 的一个变换，变换也称为映射，记为T ，即有等式 β=T(α) 称β为α在T 变换下的象，称α为β在T 变换下的源，集合A 称为变换T 的源集，A 在变换T 下的所有象称为象集，显然象集是B 的子集。 2.2、线性变换的定义 R 是数域F 上的线性空间，σ是R 的一个变换，并且满足 ()()()()() a k ka b a b a σσσσσ=+=+ 其中a,b ∈R ，k ∈F 则称σ是R 的一个线性变换（这是由R 到R 自身的一个映射）。线性变换定义的意义是，将R 的任意2个元素的和进行变换等同于将这2个元素分别进行变换后再求和，将R 的任意元素的数乘进行变换等同于将这个元素先进行变换再数乘。下面是线性变换的另一种表述方式： )()()(βσασβασl k l k +=+ F l k R ∈∈?,,,βα 实验_三_题目_用双线性变换法设计IIR 数字滤波器 第16周星期_3_第6,7节 一. 实验目的 (1)熟悉用双线性变换法设计IIR 数字滤波器的原理与方法。 (2)掌握数字滤波器的计算机仿真方法。 (3)通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。 二、实验内容、方法、设计程序及实验结果 (1)复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR 数字滤波 器的内容，用双线性变换法设计数字滤波器系统函数()z H 。其中满足本实验要求的数字滤波器系统函数为： ()() ()( )() 2 1212 1 6 1 2155.09044.013583.00106.117051.02686.1110007378.0-------+-+-+-+= z z z z z z z z H ()z H k k ∏==31 （3.1） 式中： ()() 3211212 12 1，，，k z C z B z z A z H k k k =--++=---- （3.2） 2155 .09044.03583.00106.17051.02686.109036 .0332211-==-==-===C B C B C B A ，，， 根据设计指标，调用MATLAB 信号处理工具箱buttord 和butter ，也可以得到()z H 。由公式（3.1）和（3.2）可见，滤波器()z H 由三个二阶滤波器()z H 1、 ()z H 2和()z H 3级联而成，如图3-1所示。 ()n y 图3-1 滤波器z H 的组成 此参数下的程序如下： %实验三，用双线性变换法设计 IIR数字滤波器 x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,... -38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,... 0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0]; k=1; %控制滤波循环变量 close all; %关闭全部绘图窗口 figure(1); %创建绘图窗口 subplot(2,2,1); %定位子图 1 n=0:55; %横坐标 stem(n,x,'.'); %画出枝干图 axis([0,56,-100,50]); %调整坐标 xlabel('n'); %标注横坐标 ylabel('x(n)'); %标注纵坐标 title('心电图信号采集序列x(n)');%命名该子图 B=[0.09036,2*0.09036,0.09036];%H1 滤波器的分子系数矩阵 A=[1.2686,-0.7051]; %H1滤波器的分母系数矩阵 A1=[1.0106,-0.3583]; %H2滤波器的分母系数矩阵 A2=[0.9044,-0.2155]; %H3滤波器的分母系数矩阵 while(k<=3) y=filter(B,A,x); %进行滤波 x=y; %重新赋值X 进行下一次滤波 k=k+1; %控制循环变量 if k==2 A=A1; else A=A2; end end subplot(2,2,3); %定位子图3 stem(n,y,'.'); axis([0,56,-100,50]); xlabel('n'); ylabel('y(n)'); title('三级滤波后的心电图信号(原坐标)'); subplot(2,2,2) stem(n,y,'.'); axis([0,56,-15,5]); xlabel('n'); ylabel('y(n)'); title('调整坐标后的心电图信号'); %求数字滤波器的幅频特性 A=[0.09036,2*0.09036,0.09036];%滤波器的分子系数矩阵 B1=[1,-1.2686,0.7051]; %H1滤波器的分母系数矩阵 B2=[1,-1.0106,0.3583]; %H2滤波器的分母系数矩阵线性变换思想在中学数学中的应用
用脉冲响应不变法和双线性变换法设计IIR数字滤波器汇总
浅谈高中数学线性变换的解题技巧
实验五 双线性变换法设计IIR数字滤波器
DLT 直接线性变换解法程序
性常微分方程的保线性变换及其应用
第六章_线性变换_68180769
用双线性变换法设计 IIR滤波器
浅谈线性变换对角化问题
双线性变换法IIR数字滤波器设计
Matlab+实现直接线性变换
实验三 用双线性变换法设计IIR数字滤波器
三维线性变换及其应用
实验三-用双线性变换法设计IIR数字滤波器
相关主题
文本预览