基于Chirp调频激光辐照的半透明材料光热特性参数检测
方法
技术领域
[0001] 本
发明涉及半透明材料光热特性参数检测方法,属于光热物性测量技术领域。
背景技术
[0002] 半透明材料有广泛的应用背景,如日常生活中的
水、玻璃、空气、
生物组织、
牙齿、
树脂镜片等,都属于半透明材料的范畴。此外,半透明材料还广泛应用在各种航空航天、国防军工等领域,如
航天器热防护材料、高超声速
飞行器陶瓷热防护瓦、柴油
发动机的陶瓷组件等。半透明材料光热
辐射传输作为
能量和信息传递的一种方式,其研究内容不断深入,应用范围不断拓宽,而且与其他学科的交叉也日益增多。
[0003] 半透明介质的光热物理性质是描述光热辐射传输过程的最基本的特性单位,对半透明介质光热参数场重建研究对于高尖科技领域组件
无损检测、生物组织医学成像、临床诊断
治疗等都具有重要的应用价值。其中吸收系数、散射系数和导热系数是表征半透明材料辐射传输和导热特性的重要参数,因此对半透明材料内部辐射和导热参数分布进行重构对于半透明材料在各种工业和医疗领域的研究具有重要的意义。
[0004] 现有半透明材料光热特性参数检测技术一般采用单纯的热波雷达成像或者基于SQP
算法,但是采用这两种技术得到的光热特性参数检测准确率不高。
发明内容
[0005] 本发明为解决现有半透明材料光热特性参数检测技术准确率不高的问题,提供了
基于Chirp调频激光辐照的半透明材料光热特性参数检测方法。
[0006] 本发明所述基于Chirp调频激光辐照的半透明材料光热特性参数检测方法,通过以下技术方案实现:
[0007] 步骤一、使用基于Chirp调制激光的热波雷达成像技术对半透明材料中内含物
位置进行
锁定;
[0008] 步骤二、根据背景材料光热特性和步骤一的锁
定位置设定内含物的物性参数,作为SQP算法的初始值;
[0009] 步骤三、求解正问题计算模型得到边界真实的
温度与辐射强度;然后采用SQP算法反演步骤一锁定的内含物位置,初步确定该内含物光热特性参数,即吸收系数、散射系数以及导热系数;
[0010] 步骤四、读取步骤三中的结果,将得到的内含物光热特性的初始分布作为下一步计算的初始值;
[0011] 步骤五、通过SQP算法对整个计算场的光热特性参数进行重建;
[0012] 步骤六、重复步骤五中的计算过程,直到满足下列条件之一计算结束,得到材料最终的光热特性参数;
[0013] (1)目标函数值达到
指定的计算
精度;所述目标函数由SQP算法反演得到的数据以及正问题计算模型得到数据共同确定;
[0015] 作为对上述技术方案的进一步阐述:
[0016] 进一步的,步骤一所述对半透明材料中内含物位置进行锁定具体包括:
[0017] 采用Chirp调制的辐射源照射材料,在单个Chirp调制周期中,照射到材料表面的激光强度由下式表示:
[0018]
[0019] qS=qam (2)
[0020]
[0021] 其中,qlaser表示入射激光的功率
密度;
[0022] qam表示入射激光峰值;
[0023] qS表示入射激光静态分量;
[0024] qD表示入射激光动态分量;
[0025] f0表示Chirp调制
信号的初始
频率;
[0026] fe表示Chirp调制信号的终止频率;
[0027] Ts表示Chirp调制信号的扫描周期;t表示调制时间;
[0028] 入射激光的动态分量引起材料表面产生热波雷达信号,热波雷达信号T(n′)的Chirp锁相
相位与Chirp锁相幅值由下式计算得到:
[0029]
[0030]
[0031]
[0032]
[0033] 其中,SChirp-cos表示Chirp同相相关函数;
[0034] SChirp-sin表示Chirp
正交相关函数;
[0035] A表示热波雷达信号的Chrip锁相幅值;
[0036] 表示热波雷达信号的Chrip锁相相位;
[0038] Ns表示热波雷达信号长度或
图像采集数,Ns=Ts×fs;n′=1,…,Ns。
[0039] 进一步的,步骤三中所述正问题计算模型的求解过程如下:
[0040] 用辐射导热耦合换热描述半透明材料
传热过程,边界为漫反射灰体边界,同时为
对流换热边界,
环境温度为Ta,对流换热系数为hw,所述计算场的左部表面受红外激光照射,能量转换辐射导热耦合方程用下式描述:
[0041]
[0042] 其中,ρ、cp、λ以及T分别为材料的密度、
比热容、导热系数以及温度,qr为由辐射传热引起的辐射源项,表示哈密顿算子;
[0043] 能量方程的初始条件和边界条件分别为:
[0044] T|t=0=T0 (9)
[0045] τqlaser+qr,w+qc,w=hw(Tw-Ta) (10)
[0046] 其中,τ为边界透射率,qlaser表示入射激光的功率密度,下标w表示半透明材料的边界,qr,w表示边界处辐射换热热流;qc,w表示边界处对流换热热流;Tw表示材料边界处温度;T0表示材料的初始温度;
[0047] 辐射源项能够用下述辐射传输方程求解:
[0048]
[0049] 其中,I(r,Ω)表示r位置和Ω方向的辐射强度,βe、κa、κs分别表示材料的衰减系数、吸收系数、散射系数,βe=κa+κs,Ib=σT4/π表示在温度T下的
黑体辐射强度,σ为黑体辐射常数,Φ(Ω',Ω)为散射相函数,Ω'表示入射方向;
[0050] 在直
角坐标系(x,y)下,采用离散坐标法对辐射传递方程(11)进行离散,得到:
[0051]
[0052] 其中,ξm表示为x轴方向的方向余弦,ηm表示y轴方向的方向余弦,wl表示立体角l上的方向权重,上角标l,m表示空间方向离散的第l个和第m个立体角;l、m=1,2,3,…,NΩ;NΩ为4π空间方向离散的立体角总数,Il(x,y)表示第l个立体角(x,y)处的辐射强度;Φ(Ωm,Ωl)为散射相函数;Ib(r)表示r位置处的黑体辐射强度;
[0053] 用下标e、w、s、n表示控制体P的各边界,则上式(12)变为下式:
[0054]
[0055] 其中, 表示控制体P内立体角m上的辐射强度;Δx、Δy分别表示控制体在x轴y轴的长度; 表示边界w上立体角m内的辐射强度; 表示边界s上立体角m内的辐射强度;Ωm表示第m个立体角内的入射方向;Ib,P表示控制体P内的黑体辐射强度;wm表示立体角m上的方向权重;
[0056] 半透明材料表面的辐射传输方程边界条件能够用下式表示:
[0057]
[0058] 其中,n1和n0分别表示环境和材料的折射率,γ表示壁面反射率,nw表示壁面外法向单位向量, 表示边界处立体角m内的辐射强度;
[0059] 对能量方程式(8)进行离散:
[0060]
[0061] 采用全隐格式,上式(15)左侧非稳态项的积分能够表示为:
[0062]
[0063] 其中,TP表示t+Δt时刻的控制体P的温度值, 表示t时刻该控制体P的温度值;
[0064] 方程(15)右侧扩散项变为:
[0065]
[0066] 其中,TE、TW、TS、TN分别表示控制体e、w、s、n边界的温度值;λe、λw、λs、λn分别表示控制体e、w、s、n边界的导热系数;δxe、δxw、δys、δyn分别表示控制体e、w、s、n边界的长度值;
[0067] 用S表示能量方程(15)中的源项,并将源项线性化得到下式:
[0068]
[0069] 其中,S0=κaG, G表示投射辐射;Δz表示控制体在z轴的长度;
[0071] aPTP=aETE+aWTW+aNTN+aSTS+b (19)
[0072] 其中:
[0073]
[0074]
[0075]
[0076]
[0077]
[0078] (ρcp)P表示控制体P的密度和定压比
热容的乘积;
[0079] 求解式(19),得到t+Δt时刻的控制体P的温度值TP与SP。
[0080] 进一步的,步骤三中所述SQP算法过程具体包括以下过程:
[0081] 考虑如下形式的非线性规划问题:
[0082]
[0083] 其中,F(x)是将要被优化的目标函数,ci(x)表示约束条件,m和me分别表示总约束和等式约束的数量;
[0085]
[0086] 其中,dk表示当前迭代中的搜索方向,xk表示当前的重建参数,Hk是如下拉格朗日方程的Hessian矩阵的近似:
[0087]
[0088] 式中ui为朗格朗日乘子;u表示表示拉格朗日乘子向量;
[0089] 引入如下罚函数:
[0090]
[0091] 式中ψ表示罚因子;
[0092] 重建参数能够更新为下式:
[0093] xk+1=xk+αkdk (24)
[0094] 式中αk表示步长,k表示迭代次数。
[0095] 进一步的,所述步长αk满足下式:
[0096]
[0097] 其中,ι表示一个正常数,跟据经验取0.1≤ι≤0.2;
[0098] 式(25)中:
[0099]
[0100] 进一步的,步骤三SQP算法在更新重建参数xk时,为了避免
马洛托斯效应,考虑下面二阶近似:
[0101]
[0102] 重建参数xk和搜索步长αk根据下式更新:
[0103]
[0104]
[0105] 其中, 表示搜索方向;
[0106] 值得注意的是仅仅在同时满足下式的时候才会考虑公式(27)描述的子问题:
[0107]
[0108]
[0109] 其中,ε表示给定的极小值;μ表示大于0小于1的正常数。
[0110] 进一步的,所述目标函数用下式表示:
[0111]
[0112] 其中,Iest(i′,j)表示边界反演的辐射强度,Iexa(i′,j)=SP表示边界真实的辐射强度;j=1,…,Nt;Nt表示正问题计算模型中的采样时间;i′=1,…,Nd;Nd表示边界探测点的数量;对于材料导热系数的重建,目标函数能够用下式表示:
[0113]
[0114] 其中,Test(i′,j)表示边界反演的温度,Texa(i′,j)=TP表示边界真实的温度;
[0115] 引入如下归一化均方误差NRMSE衡量重建结果精度:
[0116]
[0117] 其中,xj′表示材料真实的光热特性参数; 表示材料反演的光热特性参数,j′=1,…,N;N表示重建参数的数量。
[0118] 本发明最为突出的特点和显著的有益效果是:
[0119] 本发明所涉及的基于Chirp(
鸟声信号)调频激光辐照的半透明材料光热特性参数检测方法,结合了基于Chirp调制激光的热波雷达成像技术(TWRI)和序列二次规划算法(SQP),因此既有热波雷达成像技术快速定位内含物位置的优点,又有SQP算法准确重建材料光热特性的优点;本发明方法可以精确地对材料中的内含物的吸收系数、散射系数和导热系数进行同时重建,相比单纯的热波雷达成像技术和SQP算法更有效更准确,半透明材料光热特性参数检测准确率高达95%。
附图说明
[0120] 图1为本发明中对半透明材料检测光热特性参数的物理模型示意图;
具体实施方式
[0122] 具体实施方式一:结合图2对本实施方式进行说明,本实施方式给出的基于Chirp调频激光辐照的半透明材料光热特性参数检测方法,具体包括以下步骤:
[0123] 步骤一、使用基于Chirp(鸟声信号)调制激光的热波雷达成像技术对半透明材料中内含物位置进行锁定;
[0124] 步骤二、根据背景材料光热特性和步骤一的锁定位置设定内含物的相关物性参数,作为SQP算法的初始值;
[0125] 步骤三、求解正问题计算模型得到边界真实的温度与辐射强度;然后采用SQP算法反演步骤一锁定的内含物位置,初步确定该内含物光热特性参数,即吸收系数、散射系数以及导热系数;
[0126] 步骤四、读取步骤三中的结果,将得到的内含物光热特性的初始分布作为下一步计算的初始值;
[0127] 步骤五、通过SQP算法对整个计算场的光热特性参数进行重建;
[0128] 步骤六、重复步骤五中的计算过程,直到满足下列条件之一计算结束,得到材料最终的光热特性参数;
[0129] (1)目标函数值达到指定的计算精度;所述目标函数由SQP算法反演得到的数据以及正问题计算模型得到数据(边界真实的温度与辐射强度)共同确定;
[0130] (2)迭代步数达到最大值。
[0131] 本发明方法进行光热特性参数检测时,能够采用如图1所示物理模型;其中的实线箭头表示入射激光,虚线箭头表示测量信号,测量信号包括
光信号与热信号。
[0132] 具体实施方式二:本实施方式与具体实施方式一不同的是,步骤一所述对半透明材料中内含物位置进行锁定具体包括:
[0133] 采用Chirp调制的辐射源照射材料,可以在材料表面获得同样具有Chirp调制信号特征的热响应信号,并可以通过改变Chirp的初始频率和终止频率来获得更多的测量信号,根据测量信号可以确定材料的热物理性质和光学特性。在单个Chirp调制周期中,照射到材料表面的激光强度可由下式表示:
[0134]
[0135] qS=qam (2)
[0136]
[0137] 其中,qlaser表示入射激光的功率密度,W/m2;
[0138] qam表示入射激光峰值,W/m2;
[0139] qS表示入射激光静态分量,W/m2;
[0140] qD表示入射激光动态分量,W/m2;
[0141] f0表示Chirp调制信号的初始频率,Hz;
[0142] fe表示Chirp调制信号的终止频率,Hz;
[0143] Ts表示Chirp调制信号的扫描周期,秒;t表示调制时间;
[0144] 入射激光的静态分量引起材料表面温度的持续上升,入射激光的动态分量引起材料表面温度的振荡变化,其产生的表面温度信号又可称为热波雷达信号,采用Chirp锁相算法来提取热波雷达信号的幅值和相位信息,这是一种采用双路Chirp信号(同相/正交Chirp信号)作为参考信号的热波雷达信号特征信息提取算法,热波雷达信号T(k)的Chirp锁相相位与Chirp锁相幅值由下式计算得到:
[0145]
[0146]
[0147]
[0148]
[0149] 其中,SChirp-cos表示Chirp同相相关函数;
[0150] SChirp-sin表示Chirp正交相关函数;
[0151] A表示热波雷达信号的Chrip锁相幅值;
[0152] 表示热波雷达信号的Chrip锁相相位;
[0154] Ns表示热波雷达信号长度或图像采集数,Ns=Ts×fs;n′=1,…,Ns。
[0155] 其他步骤及参数与具体实施方式一相同。
[0156] 具体实施方式三:本实施方式与具体实施方式二不同的是,步骤三中所述正问题计算模型的求解过程如下:
[0157] 用辐射导热耦合换热描述半透明材料传热过程,边界为漫反射灰体边界,同时为对流换热边界,环境温度为Ta,对流换热系数为hw,如图1所示,包括内含物与背景材料的所述计算场的左部表面受红外激光照射,能量转换辐射导热耦合方程可以用下式描述:
[0158]
[0159] 其中,ρ、cp、λ以及T分别为材料的密度、比热容、导热系数以及温度,qr为由辐射传热引起的辐射源项,表示哈密顿算子;
[0160] 能量方程的初始条件和边界条件分别为:
[0161] T|t=0=T0 (9)
[0162] τqlaser+qr,w+qc,w=hw(Tw-Ta) (10)
[0163] 其中,τ为边界透射率,qlaser表示入射激光的功率密度,下标w表示半透明材料的边界,qr,w表示边界处辐射换热热流;qc,w表示边界处对流换热热流;Tw表示材料边界处温度;T0表示材料的初始温度;
[0164] 辐射源项能够用下述辐射传输方程求解:
[0165]
[0166] 其中,I(r,Ω)表示r位置和Ω方向的辐射强度,βe、κa、κs分别表示材料的衰减系数、吸收系数、散射系数,βe=κa+κs,Ib=σT4/π表示在温度T下的黑体辐射强度,σ为黑体辐射常数,Φ(Ω',Ω)为散射相函数,Ω'表示入射方向;
[0167] 在直角坐标系(x,y)下,采用离散坐标法对辐射传递方程(11)进行离散,得到:
[0168]
[0169] 其中,ξm表示为x轴方向(辐射传输方向)的方向余弦,ηm表示y轴方向的方向余弦,wl表示立体角l上的方向权重,上角标l,m表示空间方向离散的第l个和第m个立体角;l、m=1,2,3,…,NΩ;NΩ为4π空间方向离散的立体角总数,Il(x,y)表示第l个立体角(x,y)处的辐射强度;Φ(Ωm,Ωl)为散射相函数;Ib(r)表示r位置处的黑体辐射强度;
[0170] 用下标e、w、s、n表示控制体P的各边界,则上式(12)变为下式:
[0171]
[0172] 其中, 表示控制体P内立体角m上的辐射强度;Δx、Δy分别表示控制体在x轴y轴的长度; 表示边界w上立体角m内的辐射强度; 表示边界s上立体角m内的辐射强度;Ωm表示第m个立体角内的入射方向;Ib,P表示控制体P内的黑体辐射强度;wm表示立体角m上的方向权重;
[0173] 半透明材料表面的辐射传输方程边界条件能够用下式表示:
[0174]
[0175] 其中,n1和n0分别表示环境和材料的折射率,γ表示壁面反射率,nw表示壁面外法向单位向量, 表示边界处立体角m内的辐射强度;
[0176] 对能量方程式(8)进行离散:
[0177]
[0178] 采用全隐格式,上式(15)左侧非稳态项的积分能够表示为:
[0179]
[0180] 其中,TP表示t+Δt时刻的控制体P的温度值, 表示t时刻该控制体P的温度值;
[0181] 方程(15)右侧扩散项变为:
[0182]
[0183] 其中,TE、TW、TS、TN分别表示控制体e、w、s、n边界的温度值;λe、λw、λs、λn分别表示控制体e、w、s、n边界的导热系数;δxe、δxw、δys、δyn分别表示控制体e、w、s、n边界的长度值;
[0184] 用S表示能量方程(15)中的源项,并将源项线性化得到下式:
[0185]
[0186] 其中,S0=κaG, G表示投射辐射;Δz表示控制体在z轴的长度;
[0187] 整理以上结果得:
[0188] aPTP=aETE+aWTW+aNTN+aSTS+b (19)
[0189] 其中:
[0190]
[0191]
[0192]
[0193]
[0194]
[0195] (ρcp)P表示控制体P的密度和定压比热容的乘积;
[0196] 求解式(19),得到t+Δt时刻的控制体P的温度值TP与SP。
[0197] 其他步骤及参数与具体实施方式二相同。
[0198] 具体实施方式四:本实施方式与具体实施方式一、二或三不同的是,步骤三中所述SQP算法过程具体包括以下过程:
[0199] 考虑如下形式的非线性规划问题:
[0200]
[0201] 其中,F(x)是将要被优化的目标函数,ci(x)表示约束条件,m和me分别表示总约束和等式约束的数量;
[0202] 在SQP算法优化过程中,优化任务转化成一系列二次规划(QP)子问题,SQP算法通过求解这些QP子问题超线性地收敛到最优。方程(20)能够转化成如下形式:
[0203]
[0204] 其中,dk表示当前迭代中的搜索方向,xk表示当前的重建参数,Hk是如下拉格朗日方程的Hessian矩阵的近似:
[0205]
[0206] 式中ui为朗格朗日乘子;u表示表示拉格朗日乘子向量;
[0207] 为了提高SQP算法的全局收敛能
力,引入如下罚函数:
[0208]
[0209] 式中ψ表示罚因子;
[0210] 重建参数能够更新为下式:
[0211] xk+1=xk+αkdk (24)
[0212] 式中αk表示步长,k表示迭代次数。
[0213] 其他步骤及参数与具体实施方式一、二或三相同。
[0214] 具体实施方式五:本实施方式与具体实施方式四不同的是,所述步长αk满足下式:
[0215]
[0216] 其中,ι表示一个正常数,跟据经验取0.1≤ι≤0.2;
[0217] 式(25)中:
[0218]
[0219] 其他步骤及参数与具体实施方式四相同。
[0220] 具体实施方式六:本实施方式与具体实施方式五不同的是,步骤三SQP算法在更新重建参数xk时,为了避免马洛托斯(Maratos)效应,考虑下面二阶近似:
[0221]
[0222] 重建参数xk和搜索步长αk根据下式更新:
[0223]
[0224]
[0225] 其中, 表示式(27)的搜索方向;
[0226] 值得注意的是仅仅在同时满足下式的时候才会考虑公式(27)描述的子问题:
[0227]
[0228]
[0229] 其中,ε表示给定的极小值;μ表示大于0小于1的正常数,0<μ<1。
[0230] 其他步骤及参数与具体实施方式五相同。
[0231] 具体实施方式七:本实施方式与具体实施方式六不同的是,步骤六中在优化过程提出多阶段重建技术,根据不同的目标函数对半透明材料的光学参数和热物性参数进行分阶段重建,对于材料吸收系数和散射系数的重建,所述目标函数可用下式表示:
[0232]
[0233] 其中,Iest(i′,j)表示边界反演的辐射强度,Iexa(i′,j)=SP表示边界真实的辐射强度;j=1,…,Nt;Nt表示正问题计算模型中的采样时间;i′=1,…,Nd;Nd表示边界探测点的数量;对于材料导热系数的重建,目标函数能够用下式表示:
[0234]
[0235] 其中,Test(i′,j)表示边界反演的温度,Texa(i′,j)=TP表示边界真实的温度;
[0236] 引入如下归一化均方误差NRMSE衡量重建结果精度:
[0237]
[0238] 其中,xj′表示材料真实的光热特性参数; 表示材料反演的光热特性参数,j′=1,…,N;N表示重建参数的数量。
[0239] 其他步骤及参数与具体实施方式六相同。
[0240] 本发明还可有其它多种
实施例,在不背离本发明精神及其实质的情况下,本领域技术人员当可根据本发明作出各种相应的改变和
变形,但这些相应的改变和变形都应属于本发明所附的
权利要求的保护范围。