反演基本问题
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
1 病态矩阵
1.1概念----与奇异阵的区别
病态矩阵[1]是指求解方程组时对数据的小扰动很敏感的矩阵。解线性方程组Ax=b时,若对于系数矩阵A及右端项b的小扰动δA、δb,方程组(A+δA)χ=b+δb 的解χ与原方程组Ax=b的解差别很大,则称矩阵A为病态矩阵。
方程组的近似解χ一般都不可能恰好使剩余r=b-Aχ为零,这时χ亦可看作小扰动问题Aχ=b-r(即δA=0,δb=-r)的解,所以当A为病态时,即使剩余很小,仍可能得到一个与真解相差很大的近似解。
奇异阵就是行列式为零的矩阵。
1.2判断
A的最小奇异值可以衡量A与奇异值矩阵集合相距有多远[2]。【区别奇异值与特征值:方阵才有特征值】
条件数cond A=A∗A−1,当该式范数为欧氏范数时,
cond A=σmax
,越大则病态程度越严重
σmin
可以使用matlab中的cond函数来判断,用法c = cond(X);
norm函数也可以,即条件数的第一种定义;
已经在matlab中验证条件数为1e8数量级的病态矩阵,用以上cond()或norm()的方法结果一致。见附录程序一。
反演程序中,cond(K)=2.6073e+09
2 矩阵除法及线性方程组的解
2.1 逆矩阵inv()
在线性代数中,没有除法,只有逆矩阵。矩阵除法是MATLAB从逆矩阵的概念引申来的。先介绍逆矩阵的定义,对于任意n´n阶方阵A,如果能找到一个同阶的方阵V,使
AV=I
其中,I为n阶的单位矩阵eye(n)。则V就是A的逆阵。数学符号表示为
V=A-1
逆阵V存在的条件是A的行列式det(A)不等于0,V的最古典的求法为高斯消去法,可参阅线性代数书。MATLAB已把它做成了内部函数inv,输入
V=inv(A)
就可得到A的逆矩阵V。如果det(A)等于或很接近于零,MATLAB会显示出出错或警告信息:“A矩阵病态(ill-conditioned),结果精度不可靠”。【病态矩阵:由于计算机软硬件的原因(精度、舍入误差什么的)对矩阵求解造成很大的误差】
注意:
1.求解方程组时很少直接采用inv(),而是用左除,无论是执行时间还是数值精度上都要优于直接求逆。
2.inv(A)中A必须为方阵,方阵才具有逆矩阵。
2.2 左除"\"和右除"/"
现在来看方程D*X=B,设X为未知矩阵,在等式两端同时左乘以inv(D),即
inv(D)*D*X = inv(D)*B
等式左端inv(D)*D=I,而I*X=X,因此,上式成为
X = inv(D)*B = D\B
把D的逆阵左乘以B,MATLAB就记作D\,称之为“左除”。从D*X=B的阶数检验可知,B与D的行数相等,因此,左除时的阶数检验条件是:两矩阵的行数必须相等。
如果原始方程的未知矩阵在左而系数矩阵在右,即
X*D = B
则按上述同样的方法可以写出
X = B*inv(D) = B/D
把D的逆阵右乘以B,记作/D,称之为“右除”。同理,右除时的阶数检验条件是:两矩阵的列数必须相等。
6x=3 3x=6
6\3=0.5 6/3=2
2.3 常定、超定与欠定方程组
矩阵除法可以用来方便地解线性方程组。例如要求下列方程组的解x=[ x1; x2; x3]。
6 x1 + 3 x2 + 4 x3 = 3
-2 x1 + 5 x2 + 7 x3 = -4
8 x1 - 4 x2 - 3 x3 = -7
此式可写成矩阵形式 Ax=B,求解的MATLAB程序为
A = [6,3,4;-2,5,7;8,-4,-;
B = [3;-4;-;x = A\B
得 x = 0.6000
7.0000
-5.4000
MATLAB中的除法还可以用来解方程数不等于未知数个数的情况。比如再加上一个方程
x1+5 x2 -7x3 = 9
这时系数矩阵A1的阶数为4´3。不难看出A1的行数nA1是方程数,其列数mA1是未知数的个数,nA1>mA1,说明方程组是超定的,方程无解。照样列MATLAB程序
A1 = [6,3,4;-2,5,7;8,-4,-3;1,5,-;B1 = [3;-4;-7;;x1 = A1\B1
答案为 x1 = -0.1564
1.0095
-0.6952
它并未显示出错信息,却给出了解,这怎么可能呢?实际上,这时MATLAB 给出的是最小二乘解。把这个x1代入方程组,肯定任何一个方程都不满足,都可得出1个误差,把这4个误差的平方相加开方,称为均方差。解x1保证比其他任何解所得的均方差都小。
MATLAB 中的除法还可以用来解方程数少于未知数个数的情况,A1矩阵的nA1 3 为什么不直接得x=A +b ? 这个问题实际上是问为什么不用matlab 中的A\b ?反演遇到的一般为超定方程,从2.3节超定方程组分析看出,用matlab 直接得出的解使残差r=||b-Ax||最小;当遇到A 为严重病态时,由1.1节,r 最小仍可能得到一个与真值相差很大的近似解。 4 信噪比的定义 (1)SNR 定义为第一个回波的幅度值除以误差矢量r (ˆr b Ax =- )的标准差[3]。 (2)SNR 定义为采集回波串获取的首个数据1()y t 除以测量误差的标准差[4]。 1()/SNR y t = (3)"CAL_0526_2D09_band3"程序提出,信噪比定义为第三个回波幅度除以噪声的标准差。 5 回波串累加法提高信噪比 由于 NMR 信号强度随着累加次数α增加α 因此α次累加后,信噪比 SNR ( signal-to-noise ) 当我们把数据累加 100 次,SNR 就增加了 10 倍。这也是目前运用最广的提高核磁共振回波串信噪比的方法。 将多次测量的回波串信号进行累加,得到适合的信噪比 。另外,对于长T2分量,由于其衰减得很慢,单次测量有时会漏掉这类信息,增加回波串个数,经过多次累加后,增强了对衰减慢的长T2分量的分辨能力。