CN111221048B - 基于跨孔电阻率ct多尺度反演的孤石边界识别与成像方法 - Google Patents

基于跨孔电阻率ct多尺度反演的孤石边界识别与成像方法 Download PDF

Info

Publication number
CN111221048B
CN111221048B CN202010071585.1A CN202010071585A CN111221048B CN 111221048 B CN111221048 B CN 111221048B CN 202010071585 A CN202010071585 A CN 202010071585A CN 111221048 B CN111221048 B CN 111221048B
Authority
CN
China
Prior art keywords
inversion
characteristic parameters
boundary
boulder
cross
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010071585.1A
Other languages
English (en)
Other versions
CN111221048A (zh
Inventor
刘征宇
庞永昊
聂利超
王宁
闫冰
邓朝阳
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shandong Bai20 Huitong Engineering Technology Co ltd
Original Assignee
Shandong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shandong University filed Critical Shandong University
Priority to CN202010071585.1A priority Critical patent/CN111221048B/zh
Publication of CN111221048A publication Critical patent/CN111221048A/zh
Application granted granted Critical
Publication of CN111221048B publication Critical patent/CN111221048B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本公开提供了一种基于跨孔电阻率CT多尺度反演的孤石边界识别与成像方法，根据初步探测结果，确定地质构造与环境状态信息，并根据探测结果划分精细探测区域，即可能的孤石群赋存区；在确定的精细探测区域，根据地表场地条件布置钻孔，钻孔内布置测线和电极，形成跨孔电阻率CT的探测空间；依据孤石群探测精度需求，确定电极间距以及正反演模型的网格大小；使用全空间跨孔电阻率CT方法进行数据采集；提取迭代模型的边界特征参数，构造针对边界特征参数的多尺度反演目标函数及反演方程，求解反演方程，得到边界特征参数的反演结果，实现孤石边界识别。

Description

基于跨孔电阻率CT多尺度反演的孤石边界识别与成像方法
技术领域
本公开属于孤石边界识别与成像技术领域，涉及一种基于跨孔电阻率CT多尺度反演的孤石边界识别与成像方法。
背景技术
本部分的陈述仅仅是提供了与本公开相关的背景技术信息，不必然构成在先技术。
盾构掘进机是地铁隧道建设中常用机械。但是，盾构机在遇到孤石群等复杂地质条件时，如果不能提前探明，会发生卡机甚至塌方等灾害。因此，保证施工安全的一个必要措施是：对孤石边界准确的识别与成像。
跨孔电阻率CT是一种孔中精细探测方法，与地表电阻率探测方法相比，在精细化成像方面更具优势：①电极深入围岩，背景电磁干扰小；②深部区域的信号强度大；③数据量大，有效信息多。正是基于这种优势，跨孔电阻率CT方法被认为在隐蔽不良地质精细探测领域具有更好的应用和发展前景。
然而，跨孔电阻率CT作为一种直流电阻率法，与其他地球物理方法(如地震波法、雷达法等)相比，由于其固有的物理限制而具有较低的地质异常体界面识别能力，对中小不良地质体的空间定位和形态识别能力亟待提高和改善。
据发明人了解，目前，跨孔电阻率CT仍存在以下关键问题，尚未解决：
①由于数据不完整、不一致和噪声数据，跨孔电阻率CT数据的反演和解释通常是不准确的、非唯一的。对于跨孔电阻率CT反演，非唯一性通常导致反演陷入局部极小并产生虚假异常，这可能导致错误甚至误导地质解释。
在缺少先验信息的情况下，正则化通常采用在最小二乘反演中施加普通光滑约束的形式，光滑约束反演产生平滑变化的电阻率分布模型，以不牺牲太多成像精度的方式来稳定反演过程和降低反演的非唯一性，但地质异常界面的定位难以精确确定。
发明内容
本公开为了解决上述问题，提出了一种基于跨孔电阻率CT多尺度反演的孤石边界识别与成像方法，本公开通过提取电阻率模型的边界特征参数，将反演待求量由电阻率模型转换为边界特征参数，提升边界信息在反演中的重要程度，更有针对性的对孤石边界准确识别与成像。
根据一些实施例，本公开采用如下技术方案：
一种基于跨孔电阻率CT多尺度反演的孤石边界识别方法，包括以下步骤：
根据初步探测结果，确定地质构造与环境状态信息，并根据探测结果划分精细探测区域，即可能的孤石群赋存区；
在确定的精细探测区域，根据地表场地条件布置钻孔，钻孔内布置测线和电极，形成跨孔电阻率CT的探测空间；依据孤石群探测精度需求，确定电极间距以及正反演模型的网格大小；
使用全空间跨孔电阻率CT方法进行数据采集；
提取迭代模型的边界特征参数，构造针对边界特征参数的多尺度反演目标函数及反演方程，求解反演方程，得到边界特征参数的反演结果，实现孤石边界识别。
作为可选择的实施方式，进行探测，确定地质构造与环境状态信息时，利用高密度电法或地震勘探进行探测。可以粗粒度的获得地质构造信息和地质信息。
作为可选择的实施方式，根据地表场地条件布置钻孔时，尽量在目标体附近打孔，孔间距通常不超过孔深的0.75倍。
作为可选择的实施方式，使用全空间跨孔电阻率CT方法进行数据采集过程中，采用四极法。
作为可选择的实施方式，使用小波变换提取迭代模型的边界特征参数。
作为可选择的实施方式，针对边界特征参数的多尺度反演目标函数为：
Figure GDA0002708867400000031
等式右端第一项是数据误差项，采用最小二乘方法约束实际观测数据与理论观测数据的差距；G代表正演过程，W-1表示逆小波变换；
Figure GDA0002708867400000033
表示敏感度矩阵，其求解公式为：
Figure GDA0002708867400000032
数据d为跨孔电阻率CT方法中MN的电位差或电阻；Wm为特征参数的加权矩阵，用于调整不同尺度上边界特征参数的权重。
作为可选择的实施方式，针对边界特征参数的多尺度反演目标函数求导可得反演方程为：
Figure GDA0002708867400000041
μI为阻尼项，μ为常数，I为单位矩阵，λ为元素值均为1的列向量。
一种基于跨孔电阻率CT多尺度反演的孤石边界成像方法，包括以下步骤：
根据初步探测结果，确定地质构造与环境状态信息，并根据探测结果划分精细探测区域，即可能的孤石群赋存区；
在确定的精细探测区域，根据地表场地条件布置钻孔，钻孔内布置测线和电极，形成跨孔电阻率CT的探测空间；依据孤石群探测精度需求，确定电极间距以及正反演模型的网格大小；
使用全空间跨孔电阻率CT方法进行数据采集；
提取迭代模型的边界特征参数，构造针对边界特征参数的多尺度反演目标函数及反演方程，求解反演方程，得到边界特征参数的反演结果，实现孤石边界识别；
将所述反演结果由边界特征参数转换为电阻率模型参数，绘制探测区域的电阻率分布图像，结合探测区域的地质情况，解译孤石群赋存情况。依据电阻率分布图，通常对于高电阻率的小区域可认为是孤石。
与现有技术相比，本公开的有益效果为：
1.本公开提出了一种基于跨孔电阻率CT多尺度反演的孤石边界准确识别与成像方法，通过提取电阻率模型的边界特征参数，将反演待求量由电阻率模型转换为边界特征参数，提升边界信息在反演中的重要程度，更有针对性的对孤石边界准确识别与成像。
2.本公开在多尺度反演的目标函数中，对不同尺度的边界特征系数进行了加权，可以根据需求，更有针对的获得小尺度边界(孤石边界)和大尺度边界(光滑背景)。
3.本公开也可通过动态调整模型项权重矩阵，实现从大尺度到小尺度的递进反演，即提高了全局搜索能力，又没有丢失孤石边界的局部刻画能力。
附图说明
构成本公开的一部分的说明书附图用来提供对本公开的进一步理解，本公开的示意性实施例及其说明用于解释本公开，并不构成对本公开的不当限定。
图1是基于跨孔电阻率CT多尺度反演的孤石边界准确识别与成像方法流程图；
图2是进行数值模拟时使用的地电模型图；
图3是根据基于跨孔电阻率CT多尺度反演的孤石边界准确识别与成像方法成像结果。
具体实施方式：
下面结合附图与实施例对本公开作进一步说明。
应该指出，以下详细说明都是例示性的，旨在对本公开提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本公开所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本公开的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
一种基于跨孔电阻率CT多尺度反演的孤石边界准确识别与成像方法，如图1所示，包括以下步骤：
(1)利用其它手段(如高密度电法、地震勘探等)进行普查，确定地质构造与环境状态信息，并根据探测结果划分精细探测区域，即可能的孤石群赋存区。
(2)在选定的精细探测区域，根据地表场地条件布置钻孔，钻孔内布置测线和电极，形成跨孔电阻率CT的探测空间。依据孤石群探测精度需求，确定电极间距以及正反演模型的网格大小。
(3)使用全空间跨孔电阻率CT方法进行数据d采集。
(4)提取迭代模型的边界特征参数
Figure GDA0002708867400000061
构造针对边界特征参数的多尺度反演目标函数及反演方程，求解反演方程，得到边界特征参数的反演结果。
(5)将(4)中反演结果由边界特征参数转换为电阻率模型参数，绘制探测区域的电阻率分布图像，结合地质情况，解译孤石群赋存情况。
所述步骤(2)中的尽量在目标体附近打孔，孔间距通常不超过孔深的0.75倍。
所述步骤(3)中，数据采集通常为四极法，如bipole-bipole(AM-BN)、pole–tripole(A-BMN和A-NMB)、dipole-dipole(AB-MN)。
所述步骤(4)中，提取迭代模型的边界特征参数通常使用小波变换，公式表示为：
Figure GDA0002708867400000071
所述步骤(4)中，所述针对边界特征参数的多尺度反演目标函数具体为：
Figure GDA0002708867400000072
等式右端第一项是数据误差项，采用最小二乘方法约束实际观测数据与理论观测数据的差距。G代表正演过程，W-1表示逆小波变换。
Figure GDA0002708867400000073
表示敏感度矩阵，其求解公式为：
Figure GDA0002708867400000074
数据d通常为步骤(3)中MN的电位差或电阻(MN的电位差除以此时的供电电流)。Wm为特征参数的加权矩阵，用于调整不同尺度上边界特征参数的权重。实例中Wm取对角阵，元素为5×10-6。
所述步骤(4)中，对针对边界特征参数的多尺度反演目标函数求导可得反演方程如下(模型项取L1范数)：
Figure GDA0002708867400000075
μI为阻尼项，μ为常数，I为单位矩阵，实例中μ＝6×10-5，λ为元素值均为1的列向量。
A.地电模型图如附图2所示，反演区域设置为8.0m*16.0m,背景电阻率为100Ω·m。模型中添加了两个不同大小、不同深度的矩形状的高阻体(代表孤石)，其电阻率为1000Ω·m。
B.正演使用矩形单元双线性插值的有限元方法，采用自然边界条件。网格大小为0.5m*0.5m,电极间距设为0.5米,两根测线均布置32个电极。所用数据为bipole-bipole(AM-BN)、pole–tripole(A-BMN和A-NMB)、dipole-dipole(AB-MN)三种形式的电极排列组合。
C.得到反演结果，如附图3所示。该数值模拟结果表明，基于跨孔电阻率CT多尺度反演的识别与成像方法能够较准确的刻画孤石目标体的形态和定位。
本领域内的技术人员应明白，本公开的实施例可提供为方法、系统、或计算机程序产品。因此，本公开可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本公开可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本公开是参照根据本公开实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述仅为本公开的优选实施例而已，并不用于限制本公开，对于本领域的技术人员来说，本公开可以有各种更改和变化。凡在本公开的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。
上述虽然结合附图对本公开的具体实施方式进行了描述，但并非对本公开保护范围的限制，所属领域技术人员应该明白，在本公开的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本公开的保护范围以内。

Claims (5)

1.一种基于跨孔电阻率CT多尺度反演的孤石边界识别方法，其特征是：包括以下步骤：
根据初步探测结果，确定地质构造与环境状态信息，并根据探测结果划分精细探测区域，即可能的孤石群赋存区；
在确定的精细探测区域，根据地表场地条件布置钻孔，钻孔内布置测线和电极，形成跨孔电阻率CT的探测空间；依据孤石群探测精度需求，确定电极间距以及正反演模型的网格大小；
使用全空间跨孔电阻率CT方法进行数据采集；
使用小波变换提取迭代模型的边界特征参数，公式表示为：
Figure FDA0002708867390000011
针对边界特征参数的多尺度反演目标函数为：
Figure FDA0002708867390000012
等式右端第一项是数据误差项，采用最小二乘方法约束实际观测数据与理论观测数据的差距；G代表正演过程，W-1表示逆小波变换；
Figure FDA0002708867390000013
表示敏感度矩阵，其求解公式为：
Figure FDA0002708867390000014
数据d为跨孔电阻率CT方法中MN的电位差或电阻；Wm为特征参数的加权矩阵，用于调整不同尺度上边界特征参数的权重；
针对边界特征参数的多尺度反演目标函数求导可得反演方程为：
Figure FDA0002708867390000015
μI为阻尼项，μ为常数，I为单位矩阵，λ为元素值均为1的列向量；
提取迭代模型的边界特征参数，构造针对边界特征参数的多尺度反演目标函数及反演方程，求解反演方程，得到边界特征参数的反演结果，实现孤石边界识别。
2.如权利要求1所述的一种基于跨孔电阻率CT多尺度反演的孤石边界识别方法，其特征是：进行探测，确定地质构造与环境状态信息时，利用高密度电法或地震勘探进行探测。
3.如权利要求1所述的一种基于跨孔电阻率CT多尺度反演的孤石边界识别方法，其特征是：根据地表场地条件布置钻孔时，尽量在目标体附近打孔，孔间距通常不超过孔深的0.75倍。
4.如权利要求1所述的一种基于跨孔电阻率CT多尺度反演的孤石边界识别方法，其特征是：使用全空间跨孔电阻率CT方法进行数据采集过程中，采用四极法。
5.一种基于跨孔电阻率CT多尺度反演的孤石边界成像方法，其特征是：包括以下步骤：
根据初步探测结果，确定地质构造与环境状态信息，并根据探测结果划分精细探测区域，即可能的孤石群赋存区；
在确定的精细探测区域，根据地表场地条件布置钻孔，钻孔内布置测线和电极，形成跨孔电阻率CT的探测空间；依据孤石群探测精度需求，确定电极间距以及正反演模型的网格大小；
使用全空间跨孔电阻率CT方法进行数据采集；
使用小波变换提取迭代模型的边界特征参数，公式表示为：
Figure FDA0002708867390000021
针对边界特征参数的多尺度反演目标函数为：
Figure FDA0002708867390000031
等式右端第一项是数据误差项，采用最小二乘方法约束实际观测数据与理论观测数据的差距；G代表正演过程，W-1表示逆小波变换；
Figure FDA0002708867390000032
表示敏感度矩阵，其求解公式为：
Figure FDA0002708867390000033
数据d为跨孔电阻率CT方法中MN的电位差或电阻；Wm为特征参数的加权矩阵，用于调整不同尺度上边界特征参数的权重；
针对边界特征参数的多尺度反演目标函数求导可得反演方程为：
Figure FDA0002708867390000034
μI为阻尼项，μ为常数，I为单位矩阵，λ为元素值均为1的列向量；
提取迭代模型的边界特征参数，构造针对边界特征参数的多尺度反演目标函数及反演方程，求解反演方程，得到边界特征参数的反演结果，实现孤石边界识别；
将所述反演结果由边界特征参数转换为电阻率模型参数，绘制探测区域的电阻率分布图像，结合探测区域的地质情况，解译孤石群赋存情况。
CN202010071585.1A 2020-01-21 2020-01-21 基于跨孔电阻率ct多尺度反演的孤石边界识别与成像方法 Active CN111221048B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010071585.1A CN111221048B (zh) 2020-01-21 2020-01-21 基于跨孔电阻率ct多尺度反演的孤石边界识别与成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010071585.1A CN111221048B (zh) 2020-01-21 2020-01-21 基于跨孔电阻率ct多尺度反演的孤石边界识别与成像方法

Publications (2)

Publication Number Publication Date
CN111221048A CN111221048A (zh) 2020-06-02
CN111221048B true CN111221048B (zh) 2021-04-06

Family

ID=70831235

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010071585.1A Active CN111221048B (zh) 2020-01-21 2020-01-21 基于跨孔电阻率ct多尺度反演的孤石边界识别与成像方法

Country Status (1)

Country Link
CN (1) CN111221048B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113034638B (zh) * 2021-03-22 2022-06-17 山东大学 一种基于定向平滑约束反演的场地污染物刻画方法及系统
CN113960695A (zh) * 2021-09-18 2022-01-21 山东大学 一种复杂城市环境富水岩溶精细探查方法
CN115166842B (zh) * 2022-06-24 2024-08-09 山东大学 基于可变网格的隧道直流电阻率自适应反演方法及系统
CN116927178B (zh) * 2023-07-21 2024-05-28 广东省水利水电第三工程局有限公司 一种花岗岩孤石地质地下连续墙施工方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6489772B1 (en) * 1995-01-23 2002-12-03 The Regents Of The University Of California Borehole induction coil transmitter
CN104407394B (zh) * 2014-11-26 2015-07-22 山东大学 基于电阻率跨孔ct的地铁盾构区间孤石精细化探测方法
CN104504721A (zh) * 2015-01-08 2015-04-08 中国科学院合肥物质科学研究院 一种基于Gabor小波变换纹理描述的非结构化道路检测方法
CN105467461B (zh) * 2015-12-01 2017-11-17 长江地球物理探测（武汉）有限公司 一种利用二维视电阻率数据识别地质异常体的方法
CN105807336A (zh) * 2016-03-11 2016-07-27 山东大学 适用于多地球物理探测方法的综合孤石探测模型试验装置及试验方法

Also Published As

Publication number Publication date
CN111221048A (zh) 2020-06-02

Similar Documents

Publication Publication Date Title
CN111221048B (zh) 基于跨孔电阻率ct多尺度反演的孤石边界识别与成像方法
Dannemiller et al. A new method for determination of magnetization direction
Beiki et al. Eigenvector analysis of gravity gradient tensor to locate geologic bodies
EP2235568B1 (en) Geophysical data processing systems
AU2011256959B2 (en) Processing geophysical data
CA2733729C (en) Tracking geologic object and detecting geologic anomalies in exploration seismic data volume
Singh et al. Fuzzy constrained Lp-norm inversion of direct current resistivity data
CN115327663A (zh) 深部矿产资源勘查用空-地-井立体地球物理探测方法
Wisén et al. Combination of 1D laterally constrained inversion and 2D smooth inversion of resistivity data with a priori data from boreholes
Monteiro et al. Comparative study of local versus global methods for 1D joint inversion of direct current resistivity and time‐domain electromagnetic data
AU2013356645B2 (en) Method and system of using a data weighted electromagnetic source for time-lapse monitoring of a reservoir production or hydraulic fracturing
CN105242317B (zh) 一种纵波速度的确定方法及装置
Kim et al. 3D inversion of irregular gridded 2D electrical resistivity tomography lines: Application to sinkhole mapping at the Island of Corfu (West Greece)
Al-Hadithi et al. Using source parameter imaging technique to the aeromagnetic data to estimate the basement depth of Tharthar Lake and surrounding area in Central Iraq
Kamiński et al. Geophysical model and displacement of active landslide—an example from Jastrzębia Góra Cliff (Northern Poland)
Elwaseif et al. Reconstruction of discrete resistivity targets using coupled artificial neural networks and watershed algorithms
Hickey Application of Land-Based Controlled Source EM Method to Hydraulic Fracture Monitoring
Martinelli et al. Different techniques for the assessment of geoelectrical data errors to improve the electrical images obtained at an industrial plant
Lin et al. RANSAC-MP: A Robust Plane Fitting Algorithm for Noisy Microseismic Data During Hydraulic Fracturing
Vértesy et al. Ground-based geophysical surveys at the Bátaapáti (Üveghuta) Site and in its vicinity
Paine Developments in geophysical inversion in the last decade
Ulla et al. Automation of the SLUTH method: a novel approach to airborne magnetic data interpretation
CN117369007A (zh) 一种电磁反演数据弱信息提取方法和装置
CN116400412A (zh) 地震层位提取方法、装置、电子设备及计算机存储介质
CN118551203A (zh) 一种基于重力数据的多源信息的边界增强方法及系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20220311

Address after: 250101 Room 601, 6 / F, building a7-5, Hanyu Golden Valley, No. 7000, Jingshi Road, Jinan area, China (Shandong) pilot Free Trade Zone, Jinan City, Shandong Province

Patentee after: Shandong bai20 Huitong Engineering Technology Co.,Ltd.

Address before: 250061, No. ten, No. 17923, Lixia District, Ji'nan City, Shandong Province

Patentee before: SHANDONG University