一种基于探地雷达全波形反演的树木缺陷定量检测方法

本发明属于地球物理及林业检测，具体涉及一种基于探地雷达全波形反演的树木缺陷定量检测方法。

背景技术：

1、活立名古树木是历史文化重要载体和人类社会的宝贵资源，在经济、社会、生态、艺术多个方面扮演着重要角色，由于自然灾害、生物侵蚀和环境变化的影响，活立木会出现质量损失、内部结构组织密度和水分的变化，具体表现为树干内部出现开裂、结疤、腐朽和中空，内部缺陷的存在能够导致活立木发生塌陷或死忙，危害活立木的健康，进而给国家带来不可挽回的损失。为了及时发现树木内部缺陷并跟踪树木的健康状况，对活立木树干进行定期检测和评估成为重中之重。

2、早期的树木缺陷检测技术，主要依赖于树干外部视觉线索进行评估，包括断枝、茎膨出或茎变色，评估方法的精度、效率有限；探头技术局部检测精度高，但会对活立木造成二次伤害；为了不影响活立木正常生长，无损检测技术应运而生，包括应力波检测、超声波检测、射线检测；上述检测方法各具特点，且存在亟需解决的问题，如应力波检测对数据采集设备要求极高，超声波检测易受外界干扰，x射线检测成本过高，且易对人体造成辐射性危害。

3、探地雷达(ground penetrating radar，gpr)作为一种新兴的非介入式地球物理探测手段，近年来备受广大林业专家的青睐，其优势包括精度高、无污染、完全无损伤、操作便捷；通过向树干内部发射高频电磁波，gpr能够获取到与树干内部有关的高精度信息，实现对树干缺陷及结构的检测，具有实际意义；但由于树干表面不平整、内部结构复杂及较强的个体差异性，导致雷达波对于树干内部病害数据解译困难；为了提高对树木检测精度并实现病害的高精度成像，邸向辉和王立海详细介绍了基于gpr的树木无损检测基本步骤；肖夏阳等实现了对颐和园柳木试件的缺陷定位与二维成像；alani等构建了一种针对不规则树木的检测框架，可直接使用商用gpr进行检测。

4、尽管探地雷达用于树木病害检测已有许多成功案例，但迄今为止，基于gpr的树木缺陷检测技术大多局限于定位检测，关于缺陷高精度定量成像及缺陷定性分析工作少见报道，在一定程度上限制了探地雷达在树木检测中的应用。

技术实现思路

1、本发明的目的在于提供一种精度高、稳定性强、鲁棒性好的基于探地雷达全波形反演的树木缺陷定量检测方法。

2、本发明提供的这种基于探地雷达全波形反演的树木缺陷定量检测方法，包括如下步骤：

3、s1.针对待检测树木进行检测预处理，获取树木检测数据，同时通过探地雷达获取树木的雷达观测数据；

4、s2.采用步骤s1获取的检测数据，以及雷达观测数据，建立树木检测的初始模型；

5、s3.采用步骤s2建立得到的初始模型，通过正演模拟处理获取模拟数据；

6、s4.采用步骤s1获取的检测数据、雷达观测数据，通过确定反演频率，进而确定多频率策略；

7、s5.采用步骤s1获取的雷达观测数据，和步骤s3获取的模拟数据，构建wiener滤波器，获取步骤s4确定的反演频率下的观测数据和模拟数据，实现步骤s4确定的多频率策略；

8、s6.采用步骤s5获取的观测数据和模拟数据，构建树木内部电性模型，针对模型的反演目标函数进行求解处理，获取目标函数的梯度信息；

9、s7.采用步骤s6获取的梯度信息，计算树木内部电性模型的更新参数；

10、s8.采用步骤s7计算获得的更新参数，构建树木检测模型，完成对树木内部结构与缺陷的定量检测。

11、步骤s1所述的针对待检测树木进行检测预处理，获取树木检测数据，同时通过探地雷达获取树木的雷达观测数据，具体包括：

12、(1-1)测量待检测的树木周长，根据待检测树木情况选取设定频率的雷达天线；

13、(1-2)布置设定数量的测线，并确定检测点距；检测点距选取通过步骤(1-1)测量获得的树木周长、测量精度确定；

14、(1-3)将探地雷达底部贴紧树干表面，环绕树干一周，采集包含树干内部信息的雷达观测数据；

15、步骤s2所述的采用步骤s1获取的检测数据，以及雷达观测数据，建立树木检测的初始模型，具体包括：

16、(2-1)采用下述公式计算树干直径：

17、l＝πd

18、其中，l为树干周长，d为树干直径；

19、(2-2)通过雷达波的双程走时公式计算雷达波在树干内部的传播速度v，计算公式如下所示：

20、t＝2d/v

21、其中，t为雷达波双程旅行时间；

22、(2-3)采用下述公式计算介质的相对介电常数εr：

23、

24、其中，c＝3×108m/s为电磁波在真空介质中传播的速度；

25、(2-4)基于步骤(2-3)获取的εr，通过经验公式设置电导率σ，同时结合树木横截面形状，构建对应的横截面模型，将内部介质的相对介电常数和电导率设置为求取的设定值，获得对应的相对介电常数和电导率的初始模型；

26、步骤s3所述的采用步骤s2建立得到的初始模型，通过正演模拟处理获取模拟数据，具体包括：

27、采用基于交错网格的时域有限差分算法完成正演模拟处理，模拟实际数据的采集情况，获取模拟数据，具体包括：

28、根据二维maxwell方程与本构方程式，采用下述公式表示二维时域tm模式的解耦方程：

29、

30、其中，hx和hz均是平行于xz平面的磁场分量，ey是垂直于xz的电场分量，参数εr(f/m)是树干内部的相对介电常数，σ(h/m)是电导率，μ(s/m)是磁导率；

31、采用下述简略公式表示二维时域的tm模式的maxwell方程：

32、lu＝j

33、其中，l为正演算子，用于计算各点的波场值，计算公式如下所示：

34、

35、其中，

36、

37、u为波场向量，计算公式如下所示：

38、u＝(hx hz ey)t

39、u是与时间t和模型参数m有关的变量，可以简写为如下形式：

40、u(m,t)

41、j为电流源向量，计算公式如下所示：

42、j＝(0 0 jy)t

43、其中，jy是电流源；

44、采用下述公式计算模拟数据：

45、ei(m,rj,t)＝rui(m,t)

46、其中，ei(m,rj,t)是模拟数据，r是当前激励源的接收点的对应位置矩阵，ui是第i个激励源下的全波场；

47、步骤s4所述的采用步骤s1获取的检测数据、雷达观测数据，通过确定反演频率，进而确定多频率策略，具体包括：

48、采用步骤s1中获取的树干直径数据、观测数据频率，确定n个反演频率，按照由低频至高频的顺序排序处理，获得反演频率序列f＝(f1,f2,…,fn)，其中，f1＜f2＜…＜fn；根据反演频率序列的获得，进而确定多频率策略；

49、步骤s5所述的采用步骤s1获取的雷达观测数据，和步骤s3获取的模拟数据，构建wiener滤波器，获取步骤s4确定的反演频率下的观测数据和模拟数据，实现步骤s4确定的多频率策略，具体包括：

50、采用下述公式构建滤波器：

51、

52、其中，fwiener是频率域中的wiener低通滤波器；woriginal(ω)是原始源子波的频谱，对应为原始观测数据或模拟数据源子波的频谱；上标*表示共轭；wtarget(ω)是期望目标源子波的频谱，对应为步骤s4中确定的反演频率源子波的频谱；δ为常小数，用于预防分母为0时数值溢出；在本发明方法中，δ＝10-4；

53、利用构建的低通滤波器，分别对步骤s1获取的雷达观测数据、步骤s3获取的模拟数据进行滤波处理，获得多个反演频率下的观测数据和模拟数据，依次进行反演，实现多频率策略；

54、步骤s6所述的采用步骤s5获取的观测数据和模拟数据，构建树木内部电性模型，针对模型的反演目标函数进行求解处理，获取目标函数的梯度信息，具体包括：

55、(6-1)构建反演目标函数：

56、采用步骤s5获取的观测数据和模拟数据的差的二范数构建反演目标函数，具体的构建过程采用下述公式进行表示：

57、

58、通过对数据差的平方在源、接收点和时间三个方向上进行积分与求和，实现构建目标函数的构建；

59、其中，φd(m)表示数据目标函数项；模型介质向量m＝[εr(r),σ(r)]t表示物质电性参数，介电常数εr与电导率σ的分布，t表示转置算子；ei(m,rj,t)为正演模拟数据；为实际观测数据；m为源的个数；n为接收点的个数；rj是第j个接收点的空间坐标位置；t为时间项；

60、(6-2)引入相对电导率：

61、引入相对电导率解决双参数反演问题，具体公式如下所示：

62、σr＝σ/σ0

63、其中，σr为相对电导率，σ为电导率，σ0为参考电导率，且σ0＝1/η0，η0为自由空间波阻抗，在本发明方法中，η0＝120πω；

64、在反演中引入无量纲的比例因子β，重写尺度变换后的模型向量m与数据目标函数梯度向量gd(m)，具体公式如下所示：

65、

66、其中，mε表示相对介电常数模型；mσ表示电导率模型；gε为介电常数梯度；gσ为电导率梯度；ge为相对介电常数梯度；ε0为真空介电常数；

67、(6-3)引入全变差正则化项：

68、引入全变差正则化项(total variation,tv)，将引入tv后的目标函数采用下述公式重新表述：

69、φ(m)＝φd(m)+φm(m)

70、其中，φ(m)为新的目标函数，φd(m)为数据目标函数项，φm(m)为模型参数目标函数项，φm(m)的计算公式如下所示：

71、φm(m)＝λεtvη(mε)+λσtvη(mσ)

72、其中，λε是相对介电常数的正则化参数；λσ是电导率的正则化参数；mε是相对介电常数模型；mσ是电导率模型；tvη是全变差正则化算子，计算公式如下所示：

73、

74、其中，f为计算区域中的相对介电常数或电导率，表示梯度算子，η是为确保上述公式能够连续而设定的数值，ω为计算区域；

75、采用下述公式表示新的目标函数的梯度：

76、

77、其中，gd(m)是数据目标函数项的梯度，gm(m)是模型参数目标函数项的梯度；

78、(6-4)求解模型的反演目标函数，获取目标函数的梯度信息：

79、通过对模型参数m进行求导，完成对目标函数的求导处理，进而求解目标函数梯度；

80、(1)数据项目标函数的fréchet导数如下述公式所示：

81、

82、其中，δm是m的变分；δei(m,rj,t)是当r＝rj时δui的第三项；vit(m,rj,t)为残差，具体的计算公式如下所示：

83、

84、δm＝(δεr(r),δσ(r))t

85、其中，ei(m,rj,t)为正演模拟数据；为实际观测数据；参数εr(r)是相对介电常数；σ(r)是电导率；

86、(2)已知δu是δm引起的u的变化量，根据步骤与s3步骤中采用的计算公式，有如下计算关系：

87、lu＝j

88、

89、

90、忽略高阶项，得到下述公式：

91、

92、其中，δc为c的变分，δd为d的变分，具体的计算公式如下所示：

93、

94、进而得到上述公式的简化公式，并通过下述公式表示波场的变分δu与波场u的关系：

95、

96、(3)引入伴随场同时，定义算子l*为算子l的伴随算子，根据伴随作用，得到下述公式：

97、<l*w,δu>＝<w,lδu>

98、其中，<·>表示对时空进行积分，将上述公式处理为下述形式：

99、

100、其中，v是积分区域，t是时间变量，t0是积分时间上限；

101、定义的伴随场需要满足下述条件：

102、

103、wi(m,rj,t0)＝0

104、其中，iy是y方向上的单位向量，δ(r-rj)为空间delta函数；进而得到下述公式：

105、

106、利用delta函数的积分性质得到下述公式：

107、

108、进而有：

109、

110、其中，gε、gσ的计算公式如下所示：

111、

112、其中，ey为波场值，为伴随波场值；

113、基于上述公式可知，通过对源的正波场与残差反传产生的伴随波场积分，进而得到目标函数的梯度；

114、(4)上述公式中提出的伴随方程、伴随算子的推导公式如下所示：

115、

116、①将上式右边第一项进行如下公式计算处理：

117、

118、根据电磁波的衰减特性，当传播距离满足设定数值时有如下特性：

119、

120、进而使得等式右边的第一项为0，所以上述公式简化为：

121、

122、②将上式右边第二项进行如下公式计算处理：

123、

124、③根据初始条件：

125、δu|t＝0＝0

126、伴随场的终止条件：

127、

128、将上式右边第三项进行如下公式计算处理：

129、

130、基于上述步骤的处理，得到下述公式：

131、

132、因此得到伴随算子l*的计算公式，如下所示：

133、

134、步骤s7所述的采用步骤s6获取的梯度信息，计算树木内部电性模型的更新参数，具体包括：

135、(7-1)计算更新方向：

136、采用下述公式表示l-bfgs算法下的搜索方向pk：

137、

138、其中，bk为第k次迭代时hessian矩阵的逆矩阵，为第k次迭代时基于步骤s6最后求得的目标函数梯度；

139、(7-2)计算更新步长：

140、基于步骤(7-1)计算得到的更新方向，采用非精确线搜索中的强wolfe准则求取搜索步长αk；

141、强wolfe准则需要满足以下两个条件：

142、“充分减小条件”：

143、φ(αk)≤φ(0)+c1αkφ'(0)

144、其中，φ(αk)＝φ(mk+αkpk)，φ(0)表示初始目标函数值，φ(αk)表示第k次迭代后的目标函数值，参数c1取设定值，c1＝10-4；

145、充分减小条件用于保证函数值一定是递减的，保证存在全局收敛的可能性；

146、“曲率条件”：

147、|φ'(α)|≤c2|φ'(0)|

148、其中，|·|表示绝对值算子，参数c2取设定值，c2＝0.9；

149、曲率条件用于搜索最优步长区间内的最优步长，保证梯度小于，将步长限制在设定值的领域内；

150、步骤s8所述的采用步骤s7计算获得的更新参数，构建树木检测模型，完成对树木内部结构与缺陷的定量检测，具体包括：

151、采用步骤s7计算出的迭代更新步长αk与迭代更新方向pk，不断更新迭代模型参数，使得模拟数据与观测数据的残差最小，实现目标函数最小化；

152、采用下述公式描述模型的更新过程：

153、mk+1＝mk+δmk＝mk+αkpk

154、其中，mk为初始模型m0迭代k次后的模型，δmk为模型更新量，mk+1为新的初始模型；

155、基于构建的新的初始模型，判断模型是否满足如下迭代终止条件：

156、

157、其中，φk为第k次迭代的目标函数值，φ0为初始的目标函数值，αk为第k次迭代的更新步长；

158、如果上述模型满足迭代终止条件中的任意一个条件，则输出模型；否则重复上述步骤s3-s8直到反演结束；输出的模型为反演结果，反演结果为检测树干的横截面成像结果，通过反演结果实现对树干内部结构与缺陷的定位与定量分析。

159、本发明提供的这种基于探地雷达全波形反演的树木缺陷定量检测方法，通过提出多频率策略，将频率分量由低频到高频依次反演，防止高频数据占主导产生的多次散射等现象，降低反演的非线性；应用双参数反演策略，通过相对介电常数与电导率结果相互约束、相互验证，实现树木内部缺陷类型的准确识别，充分挖掘了雷达数据中蕴含的树木缺陷信息；而且本发明的精度高、稳定性强、鲁棒性好。