Tikhonov正则化方法在航空γ测量数据处理中的应用.pdf

T i k h o n o v正则化方法在航空测量数据处理中的应用李华锋,盛 伟,韩 斌,王雪梅,王志慧,刘文彪,李国辉(西北核技术研究所,西安7 1 0 0 2 4)摘 要:在航空测量中,当地形平坦且飞行高度变化不大时,采用传统的高度修正方法可得到良好的结果。实际测量中,地形崎岖不平、飞行高度突变等常见现象会导致高度修正后的结果仍存在偏差甚至错误。基于条带模型构设崎岖地形条件下的反演方程组,应用T i k h o n o v正则化方法求解可得到地面放射性核素的面活度浓度。T i k h o n o v正则化方法的应用结果表明:随机噪声较小时,采用广义交叉验证(g e n e r a l i z e d c r o s s v a l i d a t i o n,G C V)方法选取正则化参数得到的反演结果较好;随机噪声较大时,采用L曲线方法选取正则化参数得到的反演结果较好。与传统高度修正方法的计算结果相比该方法好,且适用于飞行高度变化较大的情形。关键词:T i k h o n o v正则化;反演;航空测量;高度修正中图分类号:T L 8 1;O 5 7 1 文献标志码:AD O I:1 0.1 2 0 6 1/j.i s s n.2 0 9 5 6 2 2 3.2 0 2 4.0 1 0 2 0 6收稿日期:2 0 2 3 0 2 1 0;修回日期:2 0 2 3 1 1 2 0基金项目:国家重点研发计划资助项目(2 0 2 0 Y F A 0 7 0 9 8 0 0)作者简介:李华锋(1 9 8 9-),男,湖北松滋人,工程师,硕士,主要从事航空放射性测量技术研究。E-m a i l:h f l i_n i n t 1 6 3.c o mA p p l i c a t i o n o f T i k h o n o v R e g u l a r i z a t i o n M e t h o d i n A i r b o r n eG a mm a R a y S p e c t r o m e t e r S u r v e y D a t a P r o c e s s i n gL I H u a f e n g S HE NG W e i HAN B i n WANG X u e m e i WANG Z h i h u i L I U W e n b i a o L I G u o h u i N o r t h w e s t I n s t i t u t e o f N u c l e a r T e c h n o l o g y X i a n 7 1 0 0 2 4 C h i n a A b s t r a c t I n a i r b o r n e g a mm a-r a y m e a s u r e m e n t s t r a d i t i o n a l h e i g h t c o r r e c t i o n m e t h o d c a n y i e l d g o o d r e s u l t s w h e n t h e t e r r a i n i s f l a t a n d t h e f l i g h t a l t i t u d e d o e s n o t v a r y s i g n i f i c a n t l y H o w e v e r c o mm o n p h e n o m e n a s u c h a s r u g g e d t e r r a i n a n d a b r u p t c h a n g e s i n f l i g h t a l t i t u d e c a n l e a d t o d e v i a t i o n s o r e v e n e r r o r s i n t h e a c t u a l m e a s u r e m e n t r e s u l t s a f t e r a p p l y i n g h e i g h t c o r r e c t i o n B y c o n s t r u c t i n g a s t r i p m o d e l u n d e r r u g g e d t e r r a i n c o n d i t i o n s a n i n v e r s i o n e q u a t i o n s e t c a n b e d e r i v e d t o e s t i m a t e t h e s u r f a c e a c t i v i t y c o n c e n t r a t i o n o f g r o u n d r a d i o a c t i v e n u c l i d e s b y u s i n g t h e T i k h o n o v r e g u l a r i z a t i o n m e t h o d T h e a p p l i c a t i o n o f t h e T i k h o n o v r e g u l a r i z a t i o n m e t h o d s h o w s t h a t w h e n t h e r a n d o m n o i s e i s s m a l l t h e i n v e r s i o n r e s u l t s o b t a i n e d b y s e l e c t i n g t h e r e g u l a r i z a t i o n p a r a m e t e r s u s i n g t h e g e n e r a l i z e d c r o s s v a l i d a t i o n G C V m e t h o d a r e b e t t e r Wh e n t h e r a n d o m n o i s e i s l a r g e t h e i n v e r s i o n r e s u l t s o b t a i n e d b y s e l e c t i n g t h e r e g u l a r i z a t i o n p a r a m e t e r s u s i n g t h e L-c u r v e m e t h o d a r e b e t t e r T h e 1-602010第1 5卷 第1期2 0 2 4年2月现 代 应 用 物 理MO D E R N A P P L I E D P HY S I C SV o l.1 5,N o.1F e b.2 0 2 4m e t h o d y i e l d s b e t t e r c o m p u t a t i o n r e s u l t s c o m p a r e d t o t r a d i t i o n a l h e i g h t c o r r e c t i o n m e t h o d s a n d i s s u i t a b l e f o r s i t u a t i o n s w h e r e t h e f l i g h t a l t i t u d e v a r i e s c o n s i d e r a b l y K e y w o r d s T i k h o n o v r e g u l a r i z a t i o n i n v e r s i o n a i r b o r n e g a mm a r a y s p e c t r o m e t e r s u r v e y h e i g h t c o r r e c t i o n 航空能谱测量的数据处理,本质是数学物理反演问题,即研究能谱测量数据与地面放射性核素面活度的映射关系。通过对计算模型的简化和近似,基于最简单的地球物理模型和射线辐射传输几何参数,国际原子能机构(I n t e r n a t i o n a l A t o m i c E n e r g y A g e n c y,I A E A)航空能谱测量技术报告中给出了标准数据处理方法1。该方法实现了对测量能谱的逐点处理,避免了对反演问题的描述和求解;但其高度修正只适用于理想情形,对于稍复杂的情形,如崎岖地形或低空飞行测量等,需对结果进行适当的修正2。D r u k e r基于条带模型,给出了一种航空能谱测量的正演问题描述,提出了离散线性不适定反演问题3。求解离散线性不适定问题,理论上最完备的是正则化方法,其中,与其他正则化方法 相 比,如T S V D(t r u n c a t e d s i n g u l a r v a l u e d e c o m p o s i t i o n)正则化,T i k h o n o v正则化方法对正则化参数的敏感性更低且正则化解更为平滑,已在地球物理反演问题中得到广泛应用4。本文首先基于条带模型和蒙特卡罗方法构设航空能谱测量的离散线性方程组,再应用T i k h o n o v正则化方法求解该离散线性不适定反演问题,通过对正则化解与标准数据处理方法的解和面活度浓度参考值的对比,对T i k h o n o v正则化方法在航空测量数据处理中的应用效果进行了分析和讨论。1航空能谱测量的离散线性不适定反演航空能谱测量的条带模型3,如图1所示。图1航空能谱测量条带模型F i g.1 S t r i p m o d e l f o r a i r b o r n e s p e c t r o m e t e r该模型假设放射源由横贯飞行航线的一系列条带组成,条带的沿飞行航线(x轴)的水平片段长度为a,b,其垂直飞行航线(y轴)的长度为2w,每个条带中的放射性核素为均匀分布。在水平面上对矩形面积la bl2w积分,得到某一条带中的放射性核素发射的射线在高度h处探测器中全能峰净计数率Nh可表示为Nh=baw-wAsP4 r2hre-arSE F A(E)R(E,)dxdy(1)其中,As为条带中放射性核素的面活度,B qm-2;E为射线的能量,k e V;P为射线的绝对强度;r为面积微元距探测器的距离,m;a 为射线在空气中的高度衰减系数,m-1;SE F A(E)为探测器的有效前面积,m2;R(E,)为探测器的角响应函数;为微元与探测器正面法线的夹角,r a d;为方位角,r a d。和分别可表示为=a r c s i nx2+y2r ,=a r c t a nyx (2)在该模型中调整各条带的高度即可近似模拟崎岖地形,条带的法线方向仍沿z轴。探测器对各条带的探测效率为Nh/As,m2B q-1s-1,对于某一飞行航线而言,可建立离散线性方程组,表示为A x=b+(3)其中:A为系数矩阵,矩阵中的元素Ai j为探测器对各条带的探测效率,m2B q-1s-1;x为活度向量,向量中的 元素xi为条 带中放射性 核素的面活 度,B qm-2;b为计数率向量,向量中的元素bi为全能峰净计数率,s-1;为计数率向量附加的随机噪声,s-1。式(3)为离散形式的1维第一类F r e d h o l m方程,该方程组的解是不连续依赖于右端项的,是不具有稳定性的不适定问题5。2T i k h o n o v正则化方法第1节中建立的离散线性不适定方程组为不定方程组,法方程组系数矩阵ATA不可逆,无法通过直接求解法方程组得到最小二乘解。正则化方法的基2-602010第1 5卷现 代 应 用 物 理本思想是,通过寻找矩阵A的广义逆矩阵A+近似矩阵A#,将不适定问题转化为一个近似的适定问题,从而得到逼近原方程组真实解的近似解,即x=A#(b+)(4)满足x=a r g m i nA x-(b+)22 (5)T i k h o n o v正则化方法是在求解过程中利用附加信息对解进行限制,使解不过分地偏离原问题的真实解的方法。理论上,有许多方法可以将解的附加信息加入到方程中,但目前最主要的一类方法是将解的2-范数最小作为附加条件加入到原问题中,同时还可以将解的初始估计值x0加入到原问题中,相当于在原问题中加入约束条件,表示为(x)=Ni=0iLi(x-x0)22(6)其中:Li为单位矩阵I或行满秩的pn阶带状矩阵;ai为约束条件的权重。则原问题的解满足x=a r g m i nA x-(b+)22+2Ni=0iLi(x-x0)22(7)其中,为正则化参数。式(3)存在最佳近似解式(4)的充分必要条件是其离散P i c a r d准则成立。在满足离散P i c a r d准则的条件下,对于任意的0,式(7)存在唯一解67。2.1T i k h o n o v正则化解对于离散形式的1维第一类F r e d h o l m方程求解,式(7)中一般取N=2,L0为单位阵,L1和L2分别为解一阶偏导算子和二阶偏导算子的离散近似,具体如下8。为保证解的光滑性,取=(1,7,1)3。L0=111;L1=1-11-11-1;L2=1-211-211-2 1 根据2-范数的定义,有Li(x-x0)22=(x-x0)TLTiLi(x-x0)(8)令矩阵L满足LTL=Ni=0aiLTiLi(9)则式(6)可变换为求解一般形式的T i k h o n o v正则化问题x=a r g m i nA x-(b+)22+2L(x-x0)22(1 0)其解为x=(2LTL+ATA)-12LTL x0+AT(b+)(1 1)初始估计值x0=0的解为x=(2LTL+ATA)-1AT(b+)(1 2)对应的正则化矩阵为A#=(2LTL+ATA)-1AT(1 3)显然,对于 非零的正则 化参数而言,矩阵2LTL+ATA可逆。对于单条飞行航线而言,建立的系数矩阵A为小规模矩阵,可通过奇异值分解(s i n g u l a r v a l u e d e c o m p o s i t i o n,S V D)法准确求出矩阵2LTL+ATA的逆矩阵,因此式(1 2)可方便地计算出。2.2正则化参数的选取正则化方法的实质是将原问题的求解转化为在数据的拟合程度和解的范数最小之间达到某种折中的问题。需通过算法确定合适的正则化参数,从而使正则化解收敛于真实解,确保解的精确性和稳定度。对于航空测量数据反演的正则化问题,正则化参数的选取应采用与先验未知误差信息无关的方法。目前研究最多的两种方法是广义交叉验证(g e n e r a l i z e d c r o s s v a l i d a t i o n,G C V)方法和L曲线方法91 0。确定正则化参数的方法各有利弊,需要在实际的不适定问题求解过程中,视情况对方法进行优选。2.2.1 G C V方法在无法确定噪声的情况下,计算正则化参数一个较好的方法是G C V方法,其基本思想起源于统计学。当计数率向量b中一个分量被去掉时(相当于方程组中去掉一个方程),由此产生的新方程组的解,能够较好地预测b中被去掉的分量。G C V方法不依赖于b的排列顺序。基于上述思想,关于T i k h o n o v正则化参数的G C V函数表示为fG C V()=(b+)-A x22t r(I-A A#)2=(I-A A#)(b+)22t r(I-A A#)2(1 4)利用G C V方法求正则化参数是指确定正则化参数,使G C V函数fG C V达到极小。3-602010 李华锋 等:T i k h o n o v正则化方法在航空测量数据处理中的应用第1期2.2.2“L”曲线方法“L”曲线方法是一种不依赖先验噪声的参数选择 方 法。该 方 法 通 过 画 出 正 则 化 解 的 范 数L x2沿着残差的范数A x-(b+)2的曲线确定 正 则 化 参 数。“L”曲 线 上 点 的 坐 标 为(l o gA x-(b+)2,l o gL x2)。“L”曲 线通常由水平部分和竖直部分构成,水平部分对应的值过大,正则化解为过度拟合的解,受残差的范数影响较大;竖直部分对应的值过小,正则化解为欠拟合的解,受解的范数影响较大。水平部分和竖直部分的交点是曲线的拐点,残差的范数和解的范数取得了折中,对应的为最优正则化参数,所求的解x为最优正则化解。拐点处的曲率最大,曲率可表示为k()=p q-p q(p)2+(q)2 32(1 5)其中:p=l o g(A x-(b+)2);q=l o gL x2。3应用结果与讨论对于在崎岖地形下6 0 m飞行高度航测1 3 7C s的场景,单 条 测 线 长 度 取 为6 k m。条 带 宽 度 取 为3 0 m,条带长度取为2 0 0 m,所测区域沿测线均匀分为2 0 0个条带。单个测量点的净计数率由1 0个条带的线性叠加已足够3。因此,可构设由1 9 1个方程、2 0 0个未知数组成的方程组。系数矩阵A由蒙特卡罗模拟计算得到,探测器在建模时选用最简单的参数。根据系数矩阵A及活度向量的参考值xt r u e即可计算出计数率向量b。叠加不同幅度的随机噪声后即可得到式(3),构设的方程组中活度向量x为待求量。计算中使用的地形参数及飞行高度如图2所示。图2地形参数和飞行高度F i g.2 T e r r a i n a n d f l i g h t a l t i t u d e分别针对等高飞行和仿地飞行附加较小随机噪声(最大振幅为计数率向量b的5%)和附加较大随机噪声(最大振 幅为计数率 向量b的2 5%)。对T i k h o n o v正则化解与标准数据处理方法的解的进行讨论。3.1正则化参数选取根据第2.2节,分别计算了等高飞行和仿地飞行的G C V函数和“L”曲线。G C V值如图3所示,“L”曲线如图4所示。(a)C o n s t a n t a l t i t u d e f l i g h t(b)G r o u n d i m i t a t i o n f l i g h t图3等高飞行和仿地飞行的G C V随的变化关系F i g.3fG C V v s.(a)C o n s t a n t a l t i t u d e f l i g h t4-602010第1 5卷现 代 应 用 物 理(b)G r o u n d i m i t a t i o n f l i g h t图4等高飞行和仿地飞行的“L”曲线F i g.4“L”c u r v e o f c o n s t a n t a l t i t u d e f l i g h ta n d g r o u n d i m i t a t i o n f l i g h t由图3可见,随机噪声越大,G C V函数在极小值处越平坦,采用G C V方法可能造成选取的正则化参数并非最优值。由图4可见,随机噪声越大,“L”曲线的“拐点”越显著,采用“L”曲线方法可较好地选取最优正则化参数。3.2T i k h o n o v正则化解根据第3.1节选取的正则化参数,等高飞行的正则化解如图5所示,仿地飞行的正则化解如图6所示。(a)5%r a n d o m n o i s e,t h e r e g u l a r i z a t i o n p a r a m e t e r i sd e t e r m i n e d b y G C V m e t h o d,=2.5 4 6(b)5%r a n d o m n o i s e,t h e r e g u l a r i z a t i o n p a r a m e t e r i sd e t e r m i n e d b y L c u r v e m e t h o d,=4.0 4 0(c)2 5%r a n d o m n o i s e,t h e r e g u l a r i z a t i o n p a r a m e t e r i sd e t e r m i n e d b y G C V m e t h o d,=6.9 3 1(d)2 5%r a n d o m n o i s e,t h e r e g u l a r i z a t i o n p a r a m e t e r i sd e t e r m i n e d b y L c u r v e m e t h o d,=1 4.9 5图5等高飞行的正则化解F i g.5 R e g u l a r i z a t i o n s o l u t i o n o f c o n s t a n t a l t i t u d e f l i g h t(a)5%r a n d o m n o i s e,t h e r e g u l a r i z a t i o n p a r a m e t e r i sd e t e r m i n e d b y G C V m e t h o d,=1.8 4 6(b)5%r a n d o m n o i s e,t h e r e g u l a r i z a t i o n p a r a m e t e r i sd e t e r m i n e d b y L c u r v e m e t h o d,=2.9 8 25-602010 李华锋 等:T i k h o n o v正则化方法在航空测量数据处理中的应用第1期(c)2 5%r a n d o m n o i s e,t h e r e g u l a r i z a t i o n p a r a m e t e r i sd e t e r m i n e d b y G C V m e t h o d,=6.3 8 7(d)2 5%r a n d o m n o i s e,t h e r e g u l a r i z a t i o n p a r a m e t e r i sd e t e r m i n e d b y L c u r v e m e t h o d,=9.6 3 9图6仿地飞行的正则化解F i g.6 R e g u l a r i z a t i o n s o l u t i o n o f g r o u n d i m i t a t i o n f l i g h t采用标准化欧氏距离(s t a n d a r d i z e d E u c l i d e a n d i s t a n c e)来度量正则化解和参考值的距离,距离越小表示正则化解越接近于参考值。正则化解和参考值的标准化欧氏距离可表示为d=nk=1xr e g,k-xt u r e,ksk 2(1 6)其中:xr e g,k为正则化解的第k个分量;xt u r e,k为参考值的第k个分量;sk为第k个分量的标准差。图5和图6中各正则化解与参考值的标准化欧氏距离如表1所列。表1正则化解与参考值的标准化欧氏距离T a b.1 S t a n d a r d i z e d E u c l i d e a n d i s t a n c e o fr e g u l a r i z a t i o n s o l u t i o n a n d r e f e r e n c e v a l u eF l i g h tm o d eR a n d o mn o i s e/1 0-2R e g u l a r i z a t i o np a r a m e t e rs e l e c t m e t h o ddC o n s t a n ta l t i t u d ef l i g h t5G C V m e t h o d4.6 6L-c u r v e m e t h o d4.7 92 5G C V m e t h o d8.7 0L-c u r v e m e t h o d7.4 3G r o u n di m i t a t i o nf l i g h t5G C V m e t h o d4.9 5L-c u r v e m e t h o d5.1 72 5G C V m e t h o d7.8 6L-c u r v e m e t h o d7.3 3 由表1可知,随机噪声较小时(5%),采用G C V方法选取的正则化参数对应的正则化解较好;随机噪声较大时(2 5%),采用“L”曲线方法选取的正则化参数对应的正则化解较好。该结论与第3.1节的结论一致。3.3正则化解与传统高度修正方法计算结果的对比通过蒙特卡罗模拟计算得到传统方法的高度修正函数,利用该函数和全能峰净计数率可计算得到面活度浓度,传统高度修正方法计算结果如图7所示。(a)5%r a n d o m n o i s e,c o n s t a n t a l t i t u d e f l i g h t(b)5%r a n d o m n o i s e,g r o u n d i m i t a t i o n f l i g h t(c)2 5%r a n d o m n o i s e,c o n s t a n t a l t i t u d e f l i g h t(d)2 5%r a n d o m n o i s e,g r o u n d i m i t a t i o n f l i g h t图7传统高度修正方法计算结果F i g.7 S o l u t i o n o f t r a d i t i o n a l h e i g h t c o r r e c t i o n m e t h o d6-602010第1 5卷现 代 应 用 物 理参考式(1 6),传统高度修正方法计算结果与参考值的标准化欧氏距离如表2所列。表2传统高度修正方法计算结果与参考值的标准化欧氏距离T a b.2 S t a n d a r d i z e d E u c l i d e a n d i s t a n c e o f t r a d i t i o n a lh e i g h t c o r r e c t i o n s o l u t i o n a n d r e f e r e n c e v a l u eF l i g h tm o d eR a n d o mn o i s e/1 0-2dC o n s t a n t a l t i t u d ef l i g h t51 0.7 42 51 2.4 5G r o u n d i m i t a t i o nf l i g h t51 2.6 32 51 4.3 3由表1和表2可知,正则化解较传统高度修正方法结果好。由图2,图5-图7可见,在高度起伏较大的区域(n=6 0附近),正则化解与参考值符合较好,而传统高度修正方法计算结果较参考值的偏差为2 02 5%,说明T i k h o n o v正则化方法适用于飞行高度变化较大的情形。4结论本文基于条带模型和蒙特卡罗方法构设航空能谱测量的离散线性方程组,应用T i k h o n o v正则化方法求解离散线性不适定反演问题,采用标准化欧氏距离来度量正则化解和参考值的距离。计算结果表明:随机噪声较小时,采用G C V方法选取的正则化参数对应的正则化解较好;随机噪声较大时,采用“L”曲线方法选取的正则化参数对应的正则化解较好。该方法计算结果比传统高度修正方法的好,且适用于飞行高度变化较大的情形。下一步可基于无人机航测系统在天然放射性核素动态测试带上开展变高飞行实验,并对该方法的应用效果进行进一步评估。参考文献 1 A i r b o r n e g a mm a r a y s p e c t r o m e t e r s u r v e y I A E A-T E C D O C-3 2 3 R V i e n n a I A E A 1 9 9 1 2 G u i d e l i n e s f o r r a d i o e l e m e n t m a p p i n g u s i n g g a mm a r a y s p e c t r o m e t r y d a t a I A E A-T E C D O C-1 3 6 3 R V i e n n a I A E A 2 0 0 3 3 D RUK E R E A i r b o r n e g a mm a-r a y s p e c t r o m e t r y d a t a u s i n g 1 5 D i n v e r s i o n J J o u r n a l o f E n v i r o n m e n t a l R a d i o a c t i v i t y 2 0 1 7 1 7 7 1 3 2 3 4 刘继军 不适定问题的正则化方法及应用 M 北京 科学出版社 2 0 0 8 L I U J i-j u n R e g u l a r i z a t i o n M e t h o d f o r I l l-p o s e d P r o b l e m a n d T h e i r A p p l i c a t i o n s M B e i j i n g S c i e n c e P r e s s 2 0 0 8 5 王宝娥 反问题中离散不适定问题的数值求解方法 D 西安 西安理工大学 2 0 0 6 WAN G B a o-e S t u d y o n n u m e r i c a l m e t h o d o f d i s c r e t e i l l-p o s e d p r o b l e m i n i n v e r s e p r o b l e m D X i a n X i a n U n i v e r s i t y o f T e c h n o l o g y 2 0 0 6 6 M I L L E R E L KA R L W C F u n d a m e n t a l s o f i n v e r s e p r o b l e m s E B O L h t t p s e c e n o r t h e a s t e r n e d u f a c-e c e e l m i l l e r e c e g 3 9 8 f 0 3 n o t e s p d f 7 F O X C N I C HO L L S G K T AN S M A n i n t r o d u c t i o n t o i n v e r s e p r o b l e m E B O L h t t p s w w w s t a t u c h i c a g o e d u g u i l l a u m e b a l P A P E R S I n t r o d u c t i o n I n v e r s e P r o b l e m s p d f 8 P R A D E E P-KUMA R K A S HANMUGHA-S UN D A R AM G A TH I RUV E NGA D A THAN R A d v a n c e s i n d e t e c t i o n a l g o r i t h m s f o r r a d i a t i o n m o n i t o r i n g J J o u r n a l o f E n v i r o n m e n t a l R a d i o a c t i v i t y 2 0 2 0 2 1 7 1 0 6 2 1 6 9 王彦飞 反演问题的计算方法及其应用 M 北京 高等教育出版社 2 0 0 7 WAN G Y a n-f e i C o m p u t a t i o n a l M e t h o d s f o r I n v e r s e P r o b l e m s a n d T h e i r A p p l i c a t i o n s M B e i j i n g H i g h e r E d u c a t i o n P r e s s 2 0 0 7 1 0 肖庭延 于慎根 王彦飞 反问题的数值解法 M 北京 科学出版社 2 0 0 3 X I AO T i n g-y a n YU S h e n-g e n WAN G Y a n-f e i N u m e r i c a l S o l u t i o n s o f I n v e r s e P r o b l e m s M B e i j i n g S c i e n c e P r e s s 2 0 0 7 7-602010 李华锋 等:T i k h o n o v正则化方法在航空测量数据处理中的应用第1期