超时间分辨冲击波速度计算方法

摘要

本发明公开了一种超时间分辨冲击波速度计算方法，包括以下步骤：S1：推导出冲击波速度与干涉条纹相位之间的解析关系；S2：推导出冲击波速度与干涉条纹相位之间的数值关系；S3：滤除数值方程解的噪声。采用以上方法，从波动方程解出发，精密证明条纹移动与速度变化之间的关系，为冲击波产生瞬间或追赶瞬间的速度计算提供了数值计算方法，将为研究极端条件下辐射冲击问题、流体力学不稳定性行为瞬时观测提供了与VISAR成像系统配套的数值计算方法。

权利要求

1.一种超时间分辨冲击波速度计算方法，其特征在于，包括以下步骤：S1：推导出冲击波速度与干涉条纹相位之间的解析关系，具体步骤包括：S11：表征条纹相机狭缝处的波动方程的解：条纹相机狭缝处波动方程的解为

为干涉条纹振幅，Aj(x,y,t)为第j支路在条纹相机狭缝上的最大光振幅，i2＝-1，t0为激光器发出光的时刻， 为光源的初始相位， 为物面上入射光与反射光之间的相位差，υ0为探针光的初始频率，Tj为支路j产生信息光的对应时刻，sj为Tj时刻光源与冲击波面之间的光程，λ0为探针光的初始波长，支路1为延时支路，支路2为非延时支路；

式(1)中

u为冲击波速度，c为光速，T条纹相机记录时刻，T0为物处于静态的任意时刻，sj为Tj时刻光源与冲击波面之间的光程，s0为T0时刻光源与物面之间的光程，Lj为Tj时刻条纹相机表面与物面之间的光程，L0为T0时刻条纹相机表面与物面之间的光程，τ为动态探测时探针光在标准具中的延迟时间，k为VISAR 成像系统的纵向放大倍率，(xj,yj)为支路j在条纹相机狭缝处动态像的任意点坐标，(x0j,y0j)为支路j在条纹相机狭缝处静态像的坐标；

S12：表征相位方程解的形式：

光场强度I(x,y,t)满足以下关系：

A1为支路1在条纹相机狭缝上的最大光振幅，A2为支路2在条纹相机狭缝上的最大光振幅，由于实验条纹图样到计算使用的条纹图样，会经历傅里叶滤波过程，此时计算使用的条纹图样分布主要由相位变化决定，故A1＝A2；

令x轴与条纹相机狭缝重合，将联立旋转变换公式 以及式(1)代入式(2)得到的相位方程解的形式：其中，(x0,y0)为像面传播方向与条纹相机狭缝完全垂直时的静态像面坐标，(0,m)为支路1与支路2的静态像面交点；

S13：将式(12)和式(13)进行化简得：其中

令标准具修正量常数

则

h为标准具厚度，n为动态探测时探针光在标准具中的折射率，n0为静态探测时探针光在标准具中的折射率，τ0为静态探测时探针光在标准具中的延迟时间，λ为动态探测时探针光的波长；

对式(4)求微分得 则

其中， 为T1时刻的冲击波速度， 为T2时刻的冲击波速度；

S14：用解析式表征冲击波速度与干涉条纹相位的关系：对式(3)求微分

上式中，由于中括号内的省略引起的相位变化远小于条纹相机记录的分辨率，中括号外的系数省略使得相位变化小于千分之一，因此在数值计算中直接使用约等号之后的形式，并在时间标度上对齐：t为激光器发出光的时间变量， 为t+τ0时刻的冲击波速度，ut为t时刻的冲击波速度；

S2：推导出冲击波速度与干涉条纹相位之间的数值关系，具体步骤包括：S21：令τ0＝4.378ι，ι为条纹相机的时间分辨率，先对式(5)求时间分辨上的算术平均，然后在整个时间域上求和并化简：S22：由于在相位记录过程中，存在着数据涨落情况，故将涨落计入到相位值中，真实相位等于记录相位减去涨落，方程(6)的差分形式重新写为：ω为随机涨落，将速度、相位值分别看成以时间变化的向量，上式写为：其中 d和p均为正整数；

故冲击波速度与干涉条纹相位的数值关系为

S3：滤除式(7)中解的噪声：

对式（7）即 进行离散傅里叶变换：

其中 m、f和ξ均为正整数， 表示第ξ列的振幅，利用形如 的窗口进行滤波。

2.根据权利要求 1所述的超时间分辨冲击波速度计算方法，其特征在于：光源相位在诊断时间内连续变化。

说明书

超时间分辨冲击波速度计算方法

技术领域

[0001] 本发明属于激光测量技术领域，具体涉及一种超时间分辨冲击波速度计算方法。

背景技术

[0002] 冲击波的传播经常伴随着介质的压强、温度、密度等物理性质的跳跃式改变。测量得到冲击波速度将为物质状态方程求解提供边界条件，是广泛用于物态研究的物理参量。尤其在惯性约束聚变领域，烧蚀内爆过程产生的冲击波是一种高速流逝过程，其产生传播、反射、衰减、消失全过程只在纳秒、皮秒量级，且其冲击传播难免受到激光等离子体相互作用、流体力学不稳定性的影响。因此，具备更高冲击波速度求解方法将为此类物理过程观察提供途径，可以为研究者进一步完善基础物理理论模型，以及从事激光脉冲整形、调束提供可靠实验参数。

[0003] 任意反射面速度干涉仪(VISAR)是用于冲击波速度信息采集的重要工具。其原理可以概括为：探针激光在冲击波界面发生反射，由于光学多普勒效应，反射光频率发生变化，因而承载冲击波速度信息；VISAR成像系统收集信息光，然后信息光在干涉仪中被分为两束，其中一束被延时，一束不被延时；最后两束光在相机阴极附近发生干涉，通过观察干涉条纹的移动便可以反推出冲击波速度。为满足不同条件下的诊断需要，VISAR装置在不断发展的同时，关于VISAR的数值计算分析也一直在完善。数值计算分析对象包含成像系统收光立体角、窗口材料、标准具的光学延时、噪声涨落等等。对于现有的VISAR技术，其时间分辨取决于系统本身的时间分辨(即标准具延迟时间，大于百皮秒)、像面记录装置的时间分辨、信号传输系统响应时间、记录装置的响应时间。随着工艺的不断改进，后面三项对VISAR装置时间分辨影响只有几皮秒，于是系统本身的时间分辨成为制约VISAR探测技术时间分辨提升的主要因素，经研究发现，系统本身的时间分辨其实是由现有算法的局限性导致，解决以上问题成为当务之急。

发明内容

[0004] 为解决现有算法的局限性限制系统本身时间分辨的技术问题，本发明提供一种超时间分辨冲击波速度计算方法，大幅提高了冲击波速度时间分辨的计算精度。

[0005] 为实现上述目的，本发明技术方案如下：

[0006] 一种超时间分辨冲击波速度计算方法，包括以下步骤：

[0007] S1：推导出冲击波速度与干涉条纹相位之间的解析关系，具体步骤包括：

[0008] S11：表征条纹相机狭缝处的波动方程的解：

[0009] 条纹相机狭缝处波动方程的解为

[0010]

[0011] 为干涉条纹振幅，Aj(x,y,t)为第j支路在条纹相机狭缝上的最大光振幅，i2＝-1，t0为激光器发出光的时刻， 为光源的初始相位， 为物面上入射光与反射光之间的相位差，υ0为探针光的初始频率，Tj为支路j产生信息光的对应时刻，sj为Tj时刻光源与冲击波面之间的光程，λ0为探针光的初始波长，支路1为延时支路，支路2为非延时支路；

[0012] 式(1)中

[0013] u为冲击波速度，c为光速，T条纹相机记录时刻，T0为物处于静态的任意时刻，sj为Tj时刻光源与冲击波面之间的光程，s0为T0时刻光源与物面之间的光程，Lj为Tj时刻条纹相机表面与物面之间的光程，L0为T0时刻条纹相机表面与物面之间的光程，τ为动态探测时探针光在标准具中的延迟时间，k为VISAR成像系统的纵向放大倍率，(xj,yj)为支路j在条纹相机狭缝处动态像的任意点坐标，(x0j,y0j)为支路j在条纹相机狭缝处静态像的坐标；

[0014] S12：表征相位方程解的形式：

[0015] 光场强度I(x,y,t)满足以下关系：

[0016]

[0017] A1为支路1在条纹相机狭缝上的最大光振幅，A2为支路2在条纹相机狭缝上的最大光振幅，由于实验条纹图样到计算使用的条纹图样，会经历傅里叶滤波过程，此时计算使用的条纹图样分布主要由相位变化决定，故A1＝A2；

[0018] 令x轴与条纹相机狭缝重合，将联立旋转变换公式 以及式(1)代入式(2)得到的相位方程解的形式：

[0019]

[0020] 其中，(x0,y0)为像面传播方向与条纹相机狭缝完全垂直时的静态像面坐标，(0,m)为支路1与支路2的静态像面交点；

[0021] S13：将式(12)和式(13)进行化简得：

[0022]

[0023] 其中

[0024] 令标准具修正量常数

[0025] 则

[0026] h为标准具厚度，n为动态探测时探针光在标准具中的折射率，n0为静态探测时探针光在标准具中的折射率，τ0为静态探测时探针光在标准具中的延迟时间，λ为动态探测时探针光的波长；

[0027] 对式(4)求微分得 则

[0028] 其中， 为T1时刻的冲击波速度， 为T2时刻的冲击波速度；

[0029] S14：用解析式表征冲击波速度与干涉条纹相位的关系：

[0030] 对式(3)

[0031] 上式中，由于中括号内的省略引起的相位变化远小于条纹相机记录的分辨率，中括号外的系数省略使得相位变化小于千分之一，因此在数值计算中直接使用约等号之后的形式，并在时间标度上对齐：

[0032]

[0033] t为激光器发出光的时间变量， 为t+τ0时刻的冲击波速度，ut为t时刻的冲击波速度；

[0034] S2：推导出冲击波速度与干涉条纹相位之间的数值关系，具体步骤包括：

[0035] S21：令τ0＝4.378ι，ι为条纹相机的时间分辨率，先对式(5)求时间分辨上的算术平均，然后在整个时间域上求和并化简：

[0036]

[0037]

[0038]

[0039]

[0040]

[0041] S22：由于在相位记录过程中，存在着数据涨落情况，故将涨落计入到相位值中，真实相位等于记录相位减去涨落，方程(6)的差分形式重新写为：

[0042]

[0043] ω为随机涨落，将速度、相位值分别看成以时间变化的向量，上式写为：

[0044] 其中 d和p均为正整数；

[0045] 故冲击波速度与干涉条纹相位的数值关系为

[0046] S3：滤除式(7)中解的噪声。

[0047] 作为优选：步骤S3中，对式(7)即 进行离散傅里叶变换：

[0048]

[0049] 其中Wm＝e(-2πi)/m，m、f和ξ均为正整数， 表示第ξ列的振幅，利用形如 的窗口进行滤波。

[0050] 采用以上方法， 中，相当于对函数 进行了加权，因此某一频谱的噪声在计算中被加强，利用形如 的窗口进行滤波，故提出了一种新的速度谱的滤波方法，通过理论计算将噪声对速度谱影响的实质展现出来，并设计滤波窗口对速度谱进行去噪。

[0051] 作为优选：光源相位在诊断时间内连续变化。因为需要保证在静态探测时具有稳定的干涉条纹，这样动态条纹相位的变化才仅与冲击波速度相关。既满足连续探测的要求，又在建立相位变化与速度的关系时有了实际意义。

[0052] 与现有技术相比，本发明的有益效果是：

[0053] 采用本发明提供的超时间分辨冲击波速度计算方法，从波动方程解出发，精密证明条纹移动与速度变化之间的关系，为冲击波产生瞬间或追赶瞬间的速度计算提供了数值计算方法，将为研究极端条件下辐射冲击问题、流体力学不稳定性行为瞬时观测提供了与VISAR成像系统配套的数值计算方法。

附图说明

[0054] 图1为VISAR成像系统的原理图；

[0055] 图2为实验数据解相位过程示意图；

[0056] 图3为窗口滤波在速度处理中的演示实例；

[0057] 图4为本发明的边值问题对数值计算的影响示意图；

[0058] 图5为本发明与传统算法的速度分布示意图。

具体实施方式

[0059] 以下结合实施例和附图对本发明作进一步说明。

[0060] 如图1所示，VISAR成像系统的基本原理如下：

[0061] 探针光通过分束镜照射在冲击波面上，由于光学多普勒效应，反射光携带冲击波速度信息，反射光经成像系统进入干涉仪，在干涉仪内被分为两束，其中一束被延时，一束不被延时，最后两束光在条纹相机狭缝附近发生干涉，干涉条纹移动反映速度信息，因此可以求解速度与条纹移动(相位变化)的关系。

[0062] 一种超时间分辨冲击波速度计算方法，包括以下步骤：

[0063] S1：推导出冲击波速度与干涉条纹相位之间的解析关系，具体步骤包括：

[0064] S11：表征条纹相机狭缝处的波动方程的解：

[0065] 条纹相机狭缝处波动方程的解为

[0066]

[0067] 为干涉条纹振幅，Aj(x,y,t)为第j支路在条纹相机狭缝上的最大光振幅，i2＝-1，t0为激光器发出光的时刻， 为光源的初始相位， 为物面上入射光与反射光之间的相位差，υ0为探针光的初始频率，Tj为支路j产生信息光的对应时刻，sj为Tj时刻光源与冲击波面之间的光程，λ0为探针光的初始波长，支路1为延时支路，支路2为非延时支路，光源相位在诊断时间内变化连续；

[0068] 式(1)中

[0069] u为冲击波速度，c为光速，T条纹相机记录时刻，T0为物处于静态的任意时刻，sj为Tj时刻光源与冲击波面之间的光程，s0为T0时刻光源与物面之间的光程，Lj为Tj时刻条纹相机表面与物面之间的光程，L0为T0时刻条纹相机表面与物面之间的光程，τ为动态探测时探针光在标准具中的延迟时间，k为VISAR成像系统的纵向放大倍率，(xj,yj)为支路j在条纹相机狭缝处动态像的任意点坐标，(x0j,y0j)为支路j在条纹相机狭缝处静态像的坐标；

[0070] S12：表征相位方程解的形式：

[0071] 光场强度I(x,y,t)满足以下关系：

[0072]

[0073] A1为支路1在条纹相机狭缝上的最大光振幅，A2为支路2在条纹相机狭缝上的最大光振幅，由于实验条纹图样到计算使用的条纹图样，会经历傅里叶滤波过程，此时计算使用的条纹图样分布主要由相位变化决定，故A1＝A2。

[0074] 请参见图2，图2为一次实验数据解相位过程，其中图2(a)为获取图样，图2(b)为取样，图2(c)为滤波，图2(d)为解相位。

[0075] 令x轴与条纹相机狭缝重合，将联立旋转变换公式 以及式(1)代入式(2)得到的相位方程解的形式：

[0076]

[0077] 其中，(x0,y0)为像面传播方向与条纹相机狭缝完全垂直时的静态像面坐标，(0,m)为支路1与支路2的静态像面交点；

[0078] 以上存在同一时间标度下的三个时刻，条纹相机记录时刻T，波面1的作用时刻T1，波面2的作用时刻T2，需要将这三个时刻统一起来，以下是三者的关系：

[0079] S13：将式(12)和式(13)进行化简得：

[0080]

[0081] 其中

[0082] 令标准具修正量常数

[0083] 则

[0084] h为标准具厚度，n为动态探测时探针光在标准具中的折射率，n0为静态探测时探针光在标准具中的折射率，τ0为静态探测时探针光在标准具中的延迟时间，λ为动态探测时探针光的波长；

[0085] 因此，当T1为均匀变化时，对式(4)求微分得：

[0086]

[0087] 则

[0088] 其中， 为T1时刻的冲击波速度， 为T2时刻的冲击波速度，能够测到的均取VISAR成像系统能够测得的最大速度200km/s，△T和△τ均为通过VISAR成像系统测得的实验数据；

[0089] S14：用解析式表征冲击波速度与干涉条纹相位的关系：

[0090] 对式(3)

[0091] 可以计算

[0092] 由于中括号内的省略引起的相位变化远小于条纹相机记录的分辨率，中括号外的系数省略使得相位变化小于千分之一，因此在数值计算中直接使用约等号之后的形式，只有当记录系统达到0.001的相位分辨时，该约等将会被重新考虑。

[0093] 进一步在时间标度上对齐：

[0094]

[0095] t为激光器发出光的时间变量， 为t+τ0时刻的冲击波速度，ut为t时刻的冲击波速度；

[0096] S2：推导出冲击波速度与干涉条纹相位之间的数值关系，具体步骤包括：

[0097] S21：解的时间精度原则上不能小于条纹相机时间分辨，因此对式(5)求时间分辨上的算术平均，然后在整个时间域上求和并化简，其中令τ0＝4.378ι，ι为条纹相机的时间分辨率。

[0098]

[0099]

[0100]

[0101]

[0102]

[0103] S22：在真实的相位记录过程中，总是存在着数据涨落情况，这些可能来自于条纹相机记录响应的涨落，也有可能来自于噪声涨落，还有可能来自与相位解缠过程中引入的算法涨落等等。这些涨落相对于相位值来说很小，使用滤波的方法也难以消除，但是，它们会对最终的数值运算解产生影响，涨落中的特定频率会被迭代求解放大，导致最终速度解失真。将涨落计入到相位值中，真实相位等于记录相位减去涨落，式(6)的差分形式重新写为：

[0104]

[0105] ω为随机涨落，将速度、相位值分别看成以时间变化的向量，上式写为：

[0106] 其中 d和p均为正整数；

[0107] 故冲击波速度与干涉条纹相位的数值关系为

[0108] S3：滤除式(7)中解的噪声，具体如下：

[0109] 对式(7)即 进行离散傅里叶变换：

[0110]

[0111] 其中Wm＝e(-2πi)/m，m、f和ξ均为正整数， 的绝大多数列的离散傅里叶变换振幅相等。

[0112] 请参见图3，图3(a)从上到下分别是图二采样的第200、400、600行的值，对其做快速傅里叶变换，傅里叶变换频谱如图3(b)。式(6)求出来的速度以及对速度的傅里叶变换频谱如图3(c)。图3(d)是两种频谱的对数坐标下的叠加。

[0113] 因此，用第ξ列对振幅项进行表示，即 ei2πζ(s)是叠加后的相位分布，ζ(s)≈0。当涨落能通过标定测量得到，式(8)可以精确求解。涨落引起的相位变化相对于真实φ的值，在动态过程可以忽略，而静态过程(静态条纹)对真实φ的获取影响大。但是静态过程的速度大小理想条件下是0，无需求解运算。

[0114] 其中， 相当于对函数 进行了加权，故某一频谱的噪声在计算中被加强。因此，使用形如 的窗口对第一项进行滤波。

[0115] 请参见图4，在滤波过程中有两个问题需要注意：

[0116] 一是 最初的几行傅里叶变换频谱与“标准窗口”差异大，比如图4(a)是第10行的值与傅里叶变换结果。导致用标准谱滤波后最初时间和最后时间附近的速度值出现较大震荡，如图4(b)，即该方法其实牺牲了边界值附近的精确速度解。

[0117] 二是窗口的系数选择。系数值与随机变量Σωp相关。观察静态条纹中Σωp取值随时间的额变化，在前100个像素点内变化了约0.1，在前200个像素点内变化了约0.5,于是Σωp总体取值可在0.6-2.0之间。图4(b)列出了第60行几种不同的滤波情况下的最终速度曲线。在不能完全确定条纹相机的涨落情况下，不妨对Σωp取为1。

[0118] 请参见图5，比较本发明提供的算法与传统算法的速度分布示意图，在不考虑新速度图的边值问题，图5(a)和图5(b)求得的速度基本一致，符合理论预期，但是区别在于新的数值计算方法将冲击加载一瞬间的速度分布求了出来(体现了时间分辨的提升)，这为以后空间内辐射冲击转换为机械冲击、冲击波追赶或加载瞬间的速度分布求解提供了途径，甚至为精密观测流体力学不稳定行为提供了观测基础；为精密诊断冲击加载过程尤其是速度变化较大情况时的速度值精密求解提供了可能；可以精确观测得到具体哪个时段产生了不稳定性的增长，以及不稳定性的演化。

[0119] 需要指出的是，本超时间分辨冲击波速度计算方法基于以下假设：

[0120] 1.观测点反射回的中心光束频率替代返回的所有光的频率。实际收集到的每一点的返回探针光立体角在8度以内，计算该点返回频率无需进一步校正。

[0121] 2.波动方程解建立在条纹相机的狭缝附近。因为干涉条纹形成区域是所关心的。

[0122] 3.狭缝附近的光场传播方向即为像面传播方向。冲击压缩中，像面始终处于景深范围内，因而像面相对位置不产生变化，但相位分布随时间和位置变化。

[0123] 最后需要说明的是，上述描述仅仅为本发明的优选实施例，本领域的普通技术人员在本发明的启示下，在不违背本发明宗旨及权利要求的前提下，可以做出多种类似的表示，这样的变换均落入本发明的保护范围之内。