[go: up one dir, main page]
More Web Proxy on the site http://driver.im/

CN114895364B - 基于温压耦合电阻率约束的深部地温场预测方法及装置 - Google Patents

基于温压耦合电阻率约束的深部地温场预测方法及装置 Download PDF

Info

Publication number
CN114895364B
CN114895364B CN202210428404.5A CN202210428404A CN114895364B CN 114895364 B CN114895364 B CN 114895364B CN 202210428404 A CN202210428404 A CN 202210428404A CN 114895364 B CN114895364 B CN 114895364B
Authority
CN
China
Prior art keywords
resistivity
temperature
pressure
normalized
logging
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
CN202210428404.5A
Other languages
English (en)
Other versions
CN114895364A (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.)
China University of Geosciences
Original Assignee
China University of Geosciences
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 China University of Geosciences filed Critical China University of Geosciences
Priority to CN202210428404.5A priority Critical patent/CN114895364B/zh
Publication of CN114895364A publication Critical patent/CN114895364A/zh
Priority to PCT/CN2022/131137 priority patent/WO2023202047A1/zh
Priority to US18/032,343 priority patent/US12158557B2/en
Application granted granted Critical
Publication of CN114895364B publication Critical patent/CN114895364B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38Processing data, e.g. for analysis, for interpretation, for correction
    • EFIXED CONSTRUCTIONS
    • E21EARTH OR ROCK DRILLING; MINING
    • E21BEARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B49/00Testing the nature of borehole walls; Formation testing; Methods or apparatus for obtaining samples of soil or well fluids, specially adapted to earth drilling or wells
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01KMEASURING TEMPERATURE; MEASURING QUANTITY OF HEAT; THERMALLY-SENSITIVE ELEMENTS NOT OTHERWISE PROVIDED FOR
    • G01K13/00Thermometers specially adapted for specific purposes
    • G01K13/10Thermometers specially adapted for specific purposes for measuring temperature within piled or stacked materials
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N25/00Investigating or analyzing materials by the use of thermal means
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01RMEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R27/00Arrangements for measuring resistance, reactance, impedance, or electric characteristics derived therefrom
    • G01R27/02Measuring real or complex resistance, reactance, impedance, or other two-pole characteristics derived therefrom, e.g. time constant
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V11/00Prospecting or detecting by methods combining techniques covered by two or more of main groups G01V1/00 - G01V9/00
    • G01V11/002Details, e.g. power supply systems for logging instruments, transmitting or recording data, specially adapted for well logging, also if the prospecting method is irrelevant
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/18Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02EREDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E10/00Energy generation through renewable energy sources
    • Y02E10/10Geothermal energy

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Mining & Mineral Resources (AREA)
  • Geophysics (AREA)
  • Environmental & Geological Engineering (AREA)
  • Data Mining & Analysis (AREA)
  • Geochemistry & Mineralogy (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Operations Research (AREA)
  • Biochemistry (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Remote Sensing (AREA)
  • Probability & Statistics with Applications (AREA)
  • Immunology (AREA)
  • Algebra (AREA)
  • General Health & Medical Sciences (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Pathology (AREA)
  • Analytical Chemistry (AREA)
  • Chemical & Material Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供一种基于温压耦合电阻率约束的深部地温场预测方法及装置,包括:将研究区或相邻区域m口钻孔的测井电阻率‑上覆地层压力‑温度数据依据深度等分成N段,且归一化处理各层段电阻率;推导不同地层不同层段归一化电阻率与温度、压力精确关系表征;反演研究区电磁数据体获取电阻率分布特征,将反演电阻率等分成M段并归一化处理为归一化反演电阻率;利用重力观测数据反演密度分布并换算成上覆地层压力;基于不同层段(深度)精确关系表征及归一化反演电阻率与上覆地层压力,便可逐层逐点计算地下深部温度场的展布特征。本发明能够基于上覆地层压力约束,将地下介质宏观电阻率特征精确转换为可视化温度场分布,预测范围广且深。

Description

基于温压耦合电阻率约束的深部地温场预测方法及装置
技术领域
本发明涉及地温场预测领域,尤其涉及一种基于温压耦合电阻率约束的深部地温场预测方法及装置。
背景技术
温度是地球内部的关键特征之一,对它的了解决定了我们研究基础地球科学问题和应用地热问题的能力。因此,最大限度地准确估计地下空间温度分布特征显得极为重要。
目前,获取地球内部温度的方式主要有两大类:直接测量及间接计算。第一类方法直接测量主要是通过钻孔测井获取沿深度方向温度特征,并基于不规则分布的钻孔测温空间插值获取区域温度场,但钻孔测温成本高且通过少量钻孔测井温度插值常常导致相当大的误差,尤其在地质构造复杂区域。第二类方法间接计算主要是基于地球化学或地球物理手段预测地温场。其中,地球化学手段主要是利用地球化学地温计来预测温度场,即通过收集到的地球化学同位素或气体成分等数据反推温度特征,虽然这类间接地温计可以预测热储温度的分布范围,但它们不能估算区域性温度分布且无法将预测温度与深度匹配。而基于地球物理探测的地温计,一方面是通过构建研究区传热模型预测地温场,另一方面是通过搭建地球物理参数(如电阻率、波速等)与温度的耦合关系并基于地球物理探测反推地温场;但前者对于温度模型的建立需要准确定义地下空间热物性参数及模型边界条件,由于这些数值及先验约束条件通常只能粗略估计,因此导致预测的温度场会有较大的误差;而后者,利用地球物理参数(如电阻率)预测温度场目前主要是使用纯经验公式,其有效性被假设为不随空间位置的变化,即经验公式中的各项参数在任意地质环境及深度下均假设为定值,很明显这种方法是不合理的。因此,现有的温度估算方法无法准确预测钻井未到达深度的温度,也无法有效提供井间空间的温度展布,更无法准确预测区域性深部地温场展布特征。
发明内容
针对目前地下温度场预测尤其深部空间温度预测所建立的温度模型误差较大,且基于纯经验公式预测地温场的不合理性。本发明提出了一种基于温压耦合电阻率约束的深部地温场预测方法及装置,通过研究区或相邻区域钻孔测井电阻率-上覆地层压力-温度数据对,推导构建地下空间不同层段(不同深度)归一化电阻率与温度、压力之间的精确关系表征;开展电磁探测数据及重力场数据精细结构反演获取电阻率及密度分布特征,并分别将电阻率进行归一化处理、密度数据换算成上覆地层压力;最终结合研究区不同深度归一化电阻率与温度、压力之间精确关系表征及归一化电阻率、上覆地层压力实现地下深部温度场预测。
为了实现上述目的,根据本发明的第一方面,本发明提出了一种基于温压耦合电阻率约束的深部地温场预测方法,所述深部地温场预测方法包括以下步骤:
获取研究区或相邻区域m口钻孔的测井电阻率-密度-温度数据对,将所述测井电阻率-密度-温度数据对中的密度换算成上覆地层压力数据,构建测井电阻率-上覆地层压力-温度数据对,依据深度将所述测井电阻率-上覆地层压力-温度数据对等分成N段,并将每口井各段测井电阻率进行归一化处理,获得不同层段归一化测井电阻率-上覆地层压力-温度数据集;
根据多组所述测井电阻率-上覆压力-温度数据集构造归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP,并计算不同井位不同层位处FRTP中的压力控制系数A(i)、温度控制系数B(i)及常系数C(i),i表示第i层段;
将所述压力控制系数A(i)、温度控制系数B(i)及常系数C(i)在不同井位不同层段的取值分别与对应层段深度构建目标数据集,对所述目标数据集进行回归分析,得到所述压力控制系数A(i)、温度控制系数B(i)及常系数C(i)随深度z的变化关系,分别表示为A(z、B(z和C(z;
对研究区电磁数据体进行三维精细反演,获取地下空间各节点P(x,y,z)的反演电阻率Rinv(x,y,z)的分布特征,并将所述反演电阻率的剖面随深度等分成M段,对每段的反演电阻率进行归一化处理,得到不同层段不同节点的归一化反演电阻率RNinv(x,y,z);同时基于重力观测数据反演三维密度分布ρg(x,y,z)并换算成上覆地层压力分布Pre(x,y,z);其中,x为地下空间节点的横向距离,y为地下空间节点的纵向距离,z为地下空间节点垂向深度;
根据地下不同深度的所述归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP(i)中包含的压力控制系数A(i)、温度控制系数B(i)与常系数C(i)及所述压力控制系数A(i)、温度控制系数B(i)及常系数C(i)随深度z的变化关系A(z)、B(z)、C(z)及所述上覆地层压力分布Pre(x,y,z),得到研究区地下深部不同层段不同节点P(x,y,z)处归一化反演电阻率与温度之间的精确关系表征,根据所述精确关系表征及所述归一化反演电阻率RNinv(x,y,z)预测研究区深部地温场T(x,y,z)展布特征。
优选地,所述将所述测井电阻率-密度-温度数据对中的密度换算成上覆地层压力数据,其中,上覆地层压力的计算公式为:
Figure BDA0003610783400000031
Pover(h)为不同深度上覆地层压力;ρlog为地层测井密度;g为重力加速度;h为地层深度;Δh为测井数据间距;NP为深度h位置处至地表处测井数据节点数。
优选地,所述依据深度将所述测井电阻率-上覆地层压力-温度数据对等分成N段,并将每口井各段测井电阻率进行归一化处理,获得不同层段归一化测井电阻率-上覆地层压力-温度数据集,其中,测井电阻率归一化处理公式为:
Figure BDA0003610783400000032
各等分层段Hrange(i)表示为:
HLmin+DLseg×(i-1)≤Hrange(i)≤HLmin+DLseg×i
Figure BDA0003610783400000033
RNorm(i,j)即为每口井在第i层段内第j个测点归一化测井电阻率,i=1,2,3…N;Rlog(i,j)为每口井第i层段内第j个测点电阻率测井数据;Rlog(i)为每口井在第i层段内所有测点测井电阻率数据集;max[Rlog(i)]表示每口井在第i层段内所有测井电阻率数据集的最大值;HLmax为每口井中测井深度的最大深度;HLmin为每口井中测井深度的最小深度;DLseg为测井等分层段的厚度,
Figure BDA0003610783400000041
表示
Figure BDA0003610783400000042
四舍五入取整数。
优选地,所述归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP(i)的表达式为:
FRTP(i)=F(RNorm(i),Ti,Piover(i)
其中,RNorm(i)、Ti及Piover三者之间耦合关系表示为:
Figure BDA0003610783400000043
FRTP(i)表示第i层段归一化电阻率-上覆地层压力-温度耦合关系约束方程,为RNorm(i)、Ti及Piover(i)的函数,表示为F(RNorm(i),Ti,Piover(i));RNorm(i)表示每口井在第i层段内所有测井归一化电阻率数据集;P0表示常压;Pover(i)表示每口井在第i层段内所有上覆地层压力数据集;Ti表示每口井在第i层段内所有测井温度数据集。
优选地,所述对每段的反演电阻率进行归一化处理,得到不同层段不同节点的归一化反演电阻率RNinv(x,y,z),其中,反演电阻率归一化处理公式为:
Figure BDA0003610783400000044
x为地下空间各节点横向距离,y为地下空间各节点纵向距离,z为地下空间各节点沿垂向深度;Rinv(k,c,y,z)表示反演电阻率数据Rinv(c,y,z)在第k层段中的数据集,k=1,2,3,…M;max[Rinv(k)]为反演电阻率数据Rinv(x,y,z)在第k层段内的最大值;RNinv(k,x,y,z)为归一化反演电阻率RNinv(c,y,z)在第k层段的数据集。
优选地,所述基于重力观测数据反演三维密度分布ρg(x,y,z)并换算成上覆地层压力分布Pre(x,y,z),其中,换算公式为:
Figure BDA0003610783400000045
Ppre(x,y,z)表示基于反演三维密度分布ρg(x,y,z)计算的上覆地层压力分布;ρg(x,y,z)为重力观测数据反演三维密度分布;g为重力加速度;Δzinv为重力反演数据沿深度z方向间距;NP为深度z位置处至地表处反演密度数据节点数。
优选地,所述根据所述精确关系表征及所述归一化反演电阻率RNinv(x,y,z)预测研究区深部地温场T(x,y,z展布特征,其中,预测研究区深部地温场展布特征的预测公式为:
Figure BDA0003610783400000051
Figure BDA0003610783400000052
Figure BDA0003610783400000053
Figure BDA0003610783400000054
T(x,y,z)为各节点温度预测值;Acal表示计算深部温度时,所采用的压力控制系数;Bcal表示计算深部温度时,所采用的温度控制系数;Ccal表示计算深部温度时,所采用的常系数。
根据本发明的第二方面,一种基于温压耦合电阻率约束的深部地温场预测装置,所述深部地温场预测装置包括以下模块:
数据集获取模块,用于获取研究区或相邻区域m口钻孔的测井电阻率-密度-温度数据对,将所述测井电阻率-密度-温度数据对中的密度换算成上覆地层压力数据,构建测井电阻率-上覆地层压力-温度数据对,依据深度将所述测井电阻率-上覆地层压力-温度数据对等分成N段,并将每口井各段测井电阻率进行归一化处理,获得不同层段归一化测井电阻率-上覆地层压力-温度数据集;
约束方程构造及系数计算模块,用于根据多组所述测井电阻率-上覆压力-温度数据集构造归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP,并计算不同井位不同层段处FRTP中的压力控制系数A(i)、温度控制系数B(i)及常系数C(i),i表示第i层段;
系数回归分析模块,用于将所述压力控制系数A(i)、温度控制系数B(i)及常系数C(i)在不同井位不同层段的取值分别与对应层段深度构建目标数据集,对所述目标数据集进行回归分析,得到所述压力控制系数A(i)、温度控制系数B(i)及常系数C(i)随深度z的变化关系,分别表示为A(z)、B(z)和C(z);
反演电阻率归一化处理模块,用于对研究区电磁数据体进行三维精细反演,获取地下空间各节点P(x,y,z)的反演电阻率Rinv(x,y,z)的分布特征,并将所述反演电阻率的剖面随深度等分成M段,对每段的反演电阻率进行归一化处理,得到不同层段不同节点的归一化反演电阻率RNinv(x,y,z);
上覆地层压力换算模块,用于反演重力观测数据三维密度分布ρg(x,y,z)并换算成上覆地层压力分布Pre(x,y,z);其中,x为地下空间节点的横向距离,y为地下空间节点的纵向距离,z为地下空间节点垂向深度;
地温场分布预测模块,用于根据地下不同深度的所述归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP中包含的压力控制系数A(i)、温度控制系数B(i)与常系数C(i)及所述压力控制系数A(i)、温度控制系数B(i)与常系数C(i)随深度z的变化关系A(z)、B(z)、C(z)及所述上覆地层压力分布Pre(x,y,z),得到研究区地下深部不同层段不同节点P(x,y,z)处归一化反演电阻率与温度之间的精确关系表征,根据所述精确关系表征及所述归一化反演电阻率RNinv(c,y,z)预测研究区深部地温场T(x,y,z)展布特征。
本发明所采取的技术方案带来的有益效果是:本发明基于上覆地层压力约束建立研究区地下不同层段不同节点的归一化反演电阻率与温度之间的精确关系表征,根据所述精确关系表征将地下介质宏观电阻率特征转换为直观的温度场分布,预测值的准确率达到了85.51%。实用性强,预测范围广且深,对于地热资源评价、地热开发温度监测等具有重要意义。
附图说明
图1为本发明的一种基于温压耦合电阻率约束的深部地温场预测方法流程图;
图2为本发明的压力控制系数A(i)、温度控制系数B(i)、常系数C(i)随深度变化特征图;
图3为本发明的研究区地下空间反演电阻率剖图;
图4为本发明的研究区地下空间上覆地层压力剖图;
图5为本发明的研究区地下空间温度场预测剖面图;
图6为本发明的研究区预测温度与实测温度对比验证图;
图7为本发明的一种基于温压耦合电阻率约束的深部地温场预测装置的结构图。
具体实施方式
为了使本发明的目的、技术方案和效果更加清楚的理解,现将结合附图对本发明实施方式作进一步的描述。
实施例一:参考图1,图1为本发明的一种基于温压耦合电阻率约束的深部地温场预测方法流程图,实施例提供的一种基于温压耦合电阻率约束的深部地温场预测方法,包括以下步骤:
S1、获取研究区或相邻区域m口钻孔的测井电阻率-密度-温度数据对,将所述测井电阻率-密度-温度数据对中的密度换算成上覆地层压力数据,构建测井电阻率-上覆地层压力-温度数据对,依据深度将所述测井电阻率-上覆地层压力-温度数据对等分成N段,并将每口井各段测井电阻率进行归一化处理,获得不同层段归一化测井电阻率-上覆地层压力-温度数据集;
在本实施例中,步骤S1具体为:利用研究区或相邻区域的m=4口钻孔的测井电阻率-密度-温度数据对,利用(1)式将m=4口钻孔的测井密度换算成上覆地层压力数据,然后将m=4口钻孔的测井电阻率-上覆地层压力-温度数据依据深度将各层段等分成DLseg=200m,并利用(2)式将各层段电阻率数据进行归一化处理:
Figure BDA0003610783400000071
Figure BDA0003610783400000072
各等分层段Hrange(i)表示为:
HLmin+DLseg×(i-1)≤Hrange(i)≤HLmin+DLseg×i (3)
Figure BDA0003610783400000073
Pover(h)为不同深度上覆地层压力;ρlog为地层测井密度;g为重力加速度;h为地层深度;Δh为测井数据间距;NP为深度h位置处至地表处测井数据节点数;RNorm(i,j)即为每口井在第i层段内第j个测点归一化测井电阻率,i=1,2,3…N;Rlog(i,j)为每口井第i层段内第j个测点电阻率测井数据;Rlog(i)为每口井在第i层段内所有测点测井电阻率数据集;max[Rlog(i)]表示每口井在第i层段内所有测井电阻率数据集的最大值;HLmax为每口井中测井深度的最大深度;HLmin为每口井中测井深度的最小深度;DLseg为测井等分层段的厚度,
Figure BDA0003610783400000081
表示
Figure BDA0003610783400000082
四舍五入取整数。
S2、基于步骤S1中归一化测井电阻率-上覆地层压力-温度数据集,通过多元回归计算不同井位不同层段的归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP中的压力控制系数A(i)、温度控制系数B(i)及常系数C(i),并获取压力控制系数、温度控制系数、常系数随地层深度的关系分别为A(h),B(h),C(h)。
在本实施例中,步骤S2具体为:基于步骤S1中归一化电阻率-上覆地层压力-温度数据集,构建归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP,通过多元回归计算不同井位不同层段FRTP(i)中的压力控制系数A(i)、温度控制系数B(i)、常系数C(i);将所述压力控制系数A(i)、温度控制系数B(i)及常系数C(i)分别与对应层段深度构建目标数据集,对所述目标数据集进行回归分析,得到所述压力控制系数、温度控制系数及常系数随深度的变化关系,分别表示为A(h),B(h),C(h)(参考图2)。其中,归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP(i)的表达式为:
FRTP(i)=F(RNorm(i),Ti,piover(i))
其中,RNorm(i)、Ti及Piover三者之间耦合关系表示为:
Figure BDA0003610783400000083
FRTP(i)表示第i层段归一化电阻率-上覆地层压力-温度耦合关系约束方程,为RNorm(i)、Ti及Piover(i)的函数,表示为F(RNorm(i),Ti,Piover(i));RNorm(i)表示每口井在第i层段内所有测井归一化电阻率数据集;P0表示常压;Pover(i)表示每口井在第i层段内所有上覆地层压力数据集;Ti表示每口井在第i层段内所有测井温度数据集。
S3、对研究区电磁数据体进行三维精细反演,获取地下空间各节点P(x,y,z)的反演电阻率ρinv(x,y,z)的分布特征,并将所述反演电阻率归一化处理为RNinv(x,y,z);其中,x为地下空间节点的横向距离,y为地下空间节点的纵向距离,z为地下空间节点垂向深度;
在本实施例中,步骤S3具体为:对研究区电磁数据体进行三维精细反演,获取地下空间各节点P(x,y,z)的反演电阻率Rinv(x,y,z)的分布特征(参考图3,是对三维反演电阻率Rinv(x,y,z)进行切片,获取某条电磁测线对应的电阻率剖面)。将所述反演电阻率剖面随深度等分成M=54段,利用(7)式对每段的反演电阻率进行归一化处理,得到不同层段不同节点的归一化反演电阻率RNinv(x,y,z)。
Figure BDA0003610783400000091
x为地下空间各节点横向距离,y为地下空间各节点纵向距离,z为地下空间各节点沿垂向深度;Rinv(k,x,y,z)表示反演电阻率数据Rinv(x,y,z)在第k(k=1,2,3,…M)层段中的数据集;max[Rinv(k)]为反演电阻率数据Rinv(x,y,z)在第k层段内的最大值;RNinv(k,x,y,z)为归一化反演电阻率RNinv(x,y,z)第k层段的数据集。
S4、基于重力观测数据反演三维密度分布ρg(x,y,z)并换算成上覆地层压力分布Pre(x,y,z);
在本实施例中,步骤S4具体为:对研究区重力数据体开展精细反演,获取地下空间各节点P(x,y,z)的密度ρg(x,y,z)分布特征,然后利用(8)式将反演密度换算成上覆地层压力Ppre(x,y,z)分布(参考图4,是对三维上覆地层压力Ppre(x,y,z)进行切片,获取步骤S3中电磁测线对应剖面的上覆地层压力分布)。
Figure BDA0003610783400000092
Ppre(x,y,z)表示基于反演三维密度分布ρg(x,y,z)计算的上覆地层压力分布;ρg(x,y,z)为重力观测数据反演三维密度分布;g为重力加速度;Δzinv为重力反演数据沿深度z方向间距;NP为深度z位置处至地表处反演密度数据节点数。
S5、基于所述地下不同深度(层段)归一化电阻率-上覆地层压力-温度耦合关系约束方程FρTP及所述Pre(x,y,z)与ρNinv(x,y,z),可预测研究区深部地温场展布特征。
在本实施例中,步骤S5具体为:基于步骤S2中归一化电阻率-上覆地层压力-温度耦合关系约束方程FρTP中包含的压力控制系数A(i)、温度控制系数B(i)、常系数C(i)及三者随深度的变化A(z)、B(z)、C(z)及步骤S4中的上覆地层压力Ppre(x,y,z),可获取研究区深部不同节点处P(x,y,z)归一化电阻率与温度之间精确关系表征,利用所述精确关系表征及步骤S3中的归一化反演电阻率数据RNinv(x,y,z),可利用式(9)逐层逐点计算地下空间各节点P(x,y,z)的温度场Y(x,y,z)展布特征(参考图5,是将地下深部三维温度场T(x,y,z)进行切片,获取步骤S3中电磁测线对应剖面的地下温度场分布特征)。
Figure BDA0003610783400000101
Figure BDA0003610783400000102
Figure BDA0003610783400000103
Figure BDA0003610783400000104
T(x,y,z)为各节点温度预测值;Acal表示计算深部温度时,所采用的压力控制系数;Bcal表示计算深部温度时,所采用的温度控制系数;Ccal表示计算深部温度时,所采用的常系数。
参考图6,图6为本发明的研究区预测温度与实测温度对比验证图;在本实施例中,提取研究区预测温度剖面中钻孔测温点位处温度随深度的变化曲线,并与实际测井温度进行对比。实验结果表明,D35钻孔预测温度与实测温度的拟合优度为R2=0.8551,由此表明,该实施例温度场预测精度为85.51%,预测温度与实测温度吻合程度较高。
实施例二:参考图7,本实施例提供了一种基于温压耦合电阻率约束的深部地温场预测装置,包括以下模块:
数据集获取模块1,用于获取研究区或相邻区域m口钻孔的测井电阻率-密度-温度数据对,将所述测井电阻率-密度-温度数据对中的密度换算成上覆地层压力数据,构建测井电阻率-上覆地层压力-温度数据对,依据深度将所述测井电阻率-上覆地层压力-温度数据对等分成N段,并将每口井各段测井电阻率进行归一化处理,获得不同层段归一化测井电阻率-上覆地层压力-温度数据集;
约束方程构造及系数计算模块2,用于根据多组所述测井电阻率-上覆压力-温度数据集构造归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP,并计算不同井位不同层段处FRTP中的压力控制系数A(i)、温度控制系数B(i)及常系数C(i),i表示第i层段;
系数回归分析模块3,用于将所述压力控制系数A(i)、温度控制系数B(i)及常系数C(i)在不同井位不同层段的取值分别与对应层段深度构建目标数据集,对所述目标数据集进行回归分析,得到所述压力控制系数A(i)、温度控制系数B(i)及常系数C(i)随深度z的变化关系,分别表示为A(z)、B(z)和C(z);
反演电阻率归一化处理模块4,用于对研究区电磁数据体进行三维精细反演,获取地下空间各节点P(x,y,z)的反演电阻率Rinv(x,y,z)的分布特征,并将所述反演电阻率的剖面随深度等分成M段,对每段的反演电阻率进行归一化处理,得到不同层段不同节点的归一化反演电阻率RNinv(x,y,z);
上覆地层压力换算模块5,用于反演重力观测数据三维密度分布ρg(x,y,z)并换算成上覆地层压力分布Pre(x,y,z);其中,x为地下空间节点的横向距离,y为地下空间节点的纵向距离,z为地下空间节点垂向深度;
地温场分布预测模块6,用于根据地下不同深度的所述归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP中包含的压力控制系数A(i)、温度控制系数B(i)与常系数C(i)及所述压力控制系数A(i)、温度控制系数B(i)与常系数C(i)随深度z的变化关系A(z)、B(z)、C(z)及所述上覆地层压力分布Pre(x,y,z),得到研究区地下深部不同层段不同节点P(x,y,z)处归一化反演电阻率与温度之间的精确关系表征,根据所述精确关系表征及所述归一化反演电阻率RNinv(x,y,z)预测研究区深部地温场T(x,y,z)展布特征。
本发明实施例提供一种基于温压耦合电阻率约束的深部地温场预测方法及装置,将研究区或相邻区域m口钻孔的测井电阻率-上覆地层压力-温度数据依据深度等分成N段,且归一化处理各层段电阻率;推导不同地层不同层段归一化电阻率与温度、压力精确关系表征FRTP(i);反演研究区电磁数据体获取电阻率分布特征,将反演电阻率等分成M段并归一化处理为RNinv(x,y,z);利用重力观测数据反演密度分布并换算成上覆地层压力Pre(x,y,z);基于所述不同层段(深度)精确关系表征FRTP及RNinv(x,y,z)与Ppre(x,y,z),便可逐层逐点计算地下深部温度场T(x,y,z)的展布特征。本发明能够基于上覆地层压力约束,将地下介质宏观电阻率特征精确转换为可视化温度场分布,预测范围广且深。
需要说明的是,在本文中,术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含,从而使得包括一系列要素的过程、方法、物品或者系统不仅包括那些要素,而且还包括没有明确列出的其他要素,或者是还包括为这种过程、方法、物品或者系统所固有的要素。在没有更多限制的情况下,由语句“包括一个……”限定的要素,并不排除在包括该要素的过程、方法、物品或者系统中还存在另外的相同要素。
上述本发明实施例序号仅仅为了描述,不代表实施例的优劣。在列举了若干装置的单元权利要求中,这些装置中的若干个可以是通过同一个硬件项来具体体现。词语第一、第二、以及第三等的使用不表示任何顺序,可将这些词语解释为标识。
以上仅为本发明的优选实施例,并非因此限制本发明的专利范围,凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换,或直接或间接运用在其他相关的技术领域,均同理包括在本发明的专利保护范围内。

Claims (8)

1.一种基于温压耦合电阻率约束的深部地温场预测方法,其特征在于,所述深部地温场预测方法包括以下步骤:
获取研究区或相邻区域m口钻孔的测井电阻率-密度-温度数据对,将所述测井电阻率-密度-温度数据对中的密度换算成上覆地层压力数据,构建测井电阻率-上覆地层压力-温度数据对,依据深度将所述测井电阻率-上覆地层压力-温度数据对等分成N段,并将每口井各段测井电阻率进行归一化处理,获得不同层段归一化测井电阻率-上覆地层压力-温度数据集;
根据多组所述测井电阻率-上覆地层压力-温度数据集构造归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP,并计算不同井位不同层位处FRTP中的压力控制系数A(i)、温度控制系数B(i)及常系数C(i),i表示第i层段;
将所述压力控制系数A(i)、温度控制系数B(i)及常系数C(i)在不同井位不同层段的取值分别与对应层段深度构建目标数据集,对所述目标数据集进行回归分析,得到所述压力控制系数A(i)、温度控制系数B(i)及常系数C(i)随深度z的变化关系,分别表示为A(z)、B(z)和C(z);
对研究区电磁数据体进行三维精细反演,获取地下空间各节点P(x,y,z)的反演电阻率Rinv(x,y,z)的分布特征,并将所述反演电阻率的剖面随深度等分成M段,对每段的反演电阻率进行归一化处理,得到不同层段不同节点的归一化反演电阻率RNinv(x,y,z);同时基于重力观测数据反演三维密度分布ρg(x,y,z)并换算成上覆地层压力分布Pre(x,y,z);其中,x为地下空间节点的横向距离,y为地下空间节点的纵向距离,z为地下空间节点垂向深度;
根据地下不同深度的所述归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP(i)中包含的压力控制系数A(i)、温度控制系数B(i)与常系数C(i)及所述压力控制系数A(i)、温度控制系数B(i)及常系数C(i)随深度z的变化关系A(z)、B(z)、C(z)及所述上覆地层压力分布Pre(x,y,z),得到研究区地下深部不同层段不同节点P(x,y,z)处归一化反演电阻率与温度之间的精确关系表征,根据所述精确关系表征及所述归一化反演电阻率RNinv(x,y,z)预测研究区深部地温场T(x,y,z)展布特征。
2.如权利要求1所述的深部地温场预测方法,其特征在于,所述将所述测井电阻率-密度-温度数据对中的密度换算成上覆地层压力数据,其中,上覆地层压力的计算公式为:
Figure FDA0004110744410000021
Pover(h)为不同深度上覆地层压力;ρlog为地层测井密度;g为重力加速度;h为地层深度;Δh为测井数据间距;NP为深度h位置处至地表处测井数据节点数。
3.如权利要求1所述的深部地温场预测方法,其特征在于,所述依据深度将所述测井电阻率-上覆地层压力-温度数据对等分成N段,并将每口井各段测井电阻率进行归一化处理,获得不同层段归一化测井电阻率-上覆地层压力-温度数据集,其中,测井电阻率归一化处理公式为:
Figure FDA0004110744410000022
各等分层段Hrange(i)表示为:
HLmin+DLseg×(i-1)≤Hrange(i)≤HLmin+DLseg×i
Figure FDA0004110744410000023
RNorm(i,j)即为每口井在第i层段内第j个测点归一化测井电阻率,i=1,2,3…N;Rlog(i,j)为每口井第i层段内第j个测点电阻率测井数据;Rlog(i)为每口井在第i层段内所有测点测井电阻率数据集;max[Rlog(i)]表示每口井在第i层段内所有测井电阻率数据集的最大值;HLmax为每口井中测井深度的最大深度;HLmin为每口井中测井深度的最小深度;DLseg为测井等分层段的厚度,
Figure FDA0004110744410000024
表示
Figure FDA0004110744410000025
四舍五入取整数。
4.如权利要求1所述的深部地温场预测方法,其特征在于,所述归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP(i)的表达式为:
FRTP(i)=F(RNorm(i),Ti,Piover(i))
其中,RNorm(i)、Ti及Piover三者之间耦合关系表示为:
Figure FDA0004110744410000026
FRTP(i)表示第i层段归一化电阻率-上覆地层压力-温度耦合关系约束方程,为RNorm(i)、Ti及Piover(i)的函数,表示为F(RNorm(i),Ti,Piover(i));RNorm(i)表示每口井在第i层段内所有测井归一化电阻率数据集;P0表示常压;Pover(i)表示每口井在第i层段内所有上覆地层压力数据集;Ti表示每口井在第i层段内所有测井温度数据集。
5.如权利要求1所述的深部地温场预测方法,其特征在于,所述对每段的反演电阻率进行归一化处理,得到不同层段不同节点的归一化反演电阻率RNinv(x,y,z),其中,反演电阻率归一化处理公式为:
Figure FDA0004110744410000031
x为地下空间各节点横向距离,y为地下空间各节点纵向距离,z为地下空间各节点沿垂向深度;Rinv(k,x,y,z)表示反演电阻率数据Rinv(x,y,z)在第k层段中的数据集,k=1,2,3,...M;max[Rinv(k)]为反演电阻率数据Rinv(x,y,z)在第k层段内的最大值;RNinv(k,x,y,z)为归一化反演电阻率RNinv(x,y,z)在第k层段的数据集。
6.如权利要求1所述的深部地温场预测方法,其特征在于,所述基于重力观测数据反演三维密度分布ρg(x,y,z)并换算成上覆地层压力分布Pre(x,y,z),其中,换算公式为:
Figure FDA0004110744410000032
Ppre(x,y,z)表示基于反演三维密度分布ρg(x,y,z)计算的上覆地层压力分布;ρg(x,y,z)为重力观测数据反演三维密度分布;g为重力加速度;Δzinv为重力反演数据沿深度z方向间距;NP为深度z位置处至地表处反演密度数据节点数。
7.如权利要求1所述的深部地温场预测方法,其特征在于,所述根据所述精确关系表征及所述归一化反演电阻率RNinv(x,y,z)预测研究区深部地温场T(x,y,z)展布特征,其中,预测研究区深部地温场展布特征的预测公式为:
Figure FDA0004110744410000033
Figure FDA0004110744410000041
Figure FDA0004110744410000042
Figure FDA0004110744410000043
T(x,y,z)为各节点温度预测值;Acal表示计算深部温度时,所采用的压力控制系数;Bcal表示计算深部温度时,所采用的温度控制系数;Ccal表示计算深部温度时,所采用的常系数。
8.一种基于温压耦合电阻率约束的深部地温场预测装置,其特征在于,所述深部地温场预测装置包括以下模块:
数据集获取模块,用于获取研究区或相邻区域m口钻孔的测井电阻率-密度-温度数据对,将所述测井电阻率-密度-温度数据对中的密度换算成上覆地层压力数据,构建测井电阻率-上覆地层压力-温度数据对,依据深度将所述测井电阻率-上覆地层压力-温度数据对等分成N段,并将每口井各段测井电阻率进行归一化处理,获得不同层段归一化测井电阻率-上覆地层压力-温度数据集;
约束方程构造及系数计算模块,用于根据多组所述测井电阻率-上覆地层压力-温度数据集构造归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP,并计算不同井位不同层段处FRTP中的压力控制系数A(i)、温度控制系数B(i)及常系数C(i),i表示第i层段;
系数回归分析模块,用于将所述压力控制系数A(i)、温度控制系数B(i)及常系数C(i)在不同井位不同层段的取值分别与对应层段深度构建目标数据集,对所述目标数据集进行回归分析,得到所述压力控制系数A(i)、温度控制系数B(i)及常系数C(i)随深度z的变化关系,分别表示为A(z)、B(z)和C(z);
反演电阻率归一化处理模块,用于对研究区电磁数据体进行三维精细反演,获取地下空间各节点P(x,y,z)的反演电阻率Rinv(x,y,z)的分布特征,并将所述反演电阻率的剖面随深度等分成M段,对每段的反演电阻率进行归一化处理,得到不同层段不同节点的归一化反演电阻率RNinv(x,y,z);
上覆地层压力换算模块,用于反演重力观测数据三维密度分布ρg(x,y,z)并换算成上覆地层压力分布Pre(x,y,z);其中,x为地下空间节点的横向距离,y为地下空间节点的纵向距离,z为地下空间节点垂向深度;
地温场分布预测模块,用于根据地下不同深度的所述归一化电阻率-上覆地层压力-温度耦合关系约束方程FRTP(i)中包含的压力控制系数A(i)、温度控制系数B(i)与常系数C(i)及所述压力控制系数A(i)、温度控制系数B(i)与常系数C(i)随深度z的变化关系A(z)、B(z)、C(z)及所述上覆地层压力分布Pre(x,y,z),得到研究区地下深部不同层段不同节点P(x,y,z)处归一化反演电阻率与温度之间的精确关系表征,根据所述精确关系表征及所述归一化反演电阻率RNinv(x,y,z)预测研究区深部地温场T(x,y,z)展布特征。
CN202210428404.5A 2022-04-22 2022-04-22 基于温压耦合电阻率约束的深部地温场预测方法及装置 Active CN114895364B (zh)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CN202210428404.5A CN114895364B (zh) 2022-04-22 2022-04-22 基于温压耦合电阻率约束的深部地温场预测方法及装置
PCT/CN2022/131137 WO2023202047A1 (zh) 2022-04-22 2022-11-10 基于温压耦合电阻率约束的深部地温场预测方法及装置
US18/032,343 US12158557B2 (en) 2022-04-22 2022-11-10 Method and device for predicting deep geothermal field based on temperature, pressure and resistivity coupling constraint

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210428404.5A CN114895364B (zh) 2022-04-22 2022-04-22 基于温压耦合电阻率约束的深部地温场预测方法及装置

Publications (2)

Publication Number Publication Date
CN114895364A CN114895364A (zh) 2022-08-12
CN114895364B true CN114895364B (zh) 2023-04-18

Family

ID=82717893

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210428404.5A Active CN114895364B (zh) 2022-04-22 2022-04-22 基于温压耦合电阻率约束的深部地温场预测方法及装置

Country Status (3)

Country Link
US (1) US12158557B2 (zh)
CN (1) CN114895364B (zh)
WO (1) WO2023202047A1 (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6502037B1 (en) * 1999-04-02 2002-12-31 Conoco Inc. Method for gravity and magnetic data inversion using vector and tensor data with seismic imaging and geopressure prediction for oil, gas and mineral exploration and production
CN106401574A (zh) * 2015-07-28 2017-02-15 中国石油化工股份有限公司 一种钻前高温地热井地层压力的预测方法
CN106814388A (zh) * 2016-12-27 2017-06-09 中国石油大学(北京) 砂泥岩储层地层压力的地震预测方法及装置
CN109667573A (zh) * 2018-12-12 2019-04-23 中国石油化工股份有限公司江汉油田分公司勘探开发研究院 三维页岩储层孔隙压力预测方法、装置和电子设备
CN113176618A (zh) * 2021-03-23 2021-07-27 中国地质大学(武汉) 一种深部地温场预测方法、设备及存储介质

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2563048A (en) * 2017-06-01 2018-12-05 Equinor Energy As Method of calculating temperature and porosity of geological structure
CN109931054A (zh) 2018-12-27 2019-06-25 西南石油大学 致密砂岩储层压力的预测方法
CN110566171A (zh) 2019-07-15 2019-12-13 西南石油大学 一种超高压致密裂缝性砂岩气藏出砂预测方法
CN110727031B (zh) * 2019-11-18 2021-04-13 科吉思石油技术咨询(北京)有限公司 一种基于三维叠前地震反演结果的地应力获取方法
US20210264262A1 (en) 2020-02-21 2021-08-26 Saudi Arabian Oil Company Physics-constrained deep learning joint inversion
CN112861065B (zh) * 2021-02-10 2022-07-08 湖北省地震局(中国地震局地震研究所) 一种自适应的震源区重力密度变化实时反演计算方法
CN113625364B (zh) 2021-08-18 2022-08-02 西南石油大学 一种基于双重校正的泥页岩地层孔隙压力的计算方法
CN113901681B (zh) 2021-09-22 2022-09-09 中国石油大学(华东) 一种全寿命周期页岩气储层双甜点三维可压性评估方法
CN114895365B (zh) * 2022-04-22 2023-03-28 中国地质大学(武汉) 高温高压岩心电阻率标定的深部地温场预测方法及装置

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6502037B1 (en) * 1999-04-02 2002-12-31 Conoco Inc. Method for gravity and magnetic data inversion using vector and tensor data with seismic imaging and geopressure prediction for oil, gas and mineral exploration and production
CN106401574A (zh) * 2015-07-28 2017-02-15 中国石油化工股份有限公司 一种钻前高温地热井地层压力的预测方法
CN106814388A (zh) * 2016-12-27 2017-06-09 中国石油大学(北京) 砂泥岩储层地层压力的地震预测方法及装置
CN109667573A (zh) * 2018-12-12 2019-04-23 中国石油化工股份有限公司江汉油田分公司勘探开发研究院 三维页岩储层孔隙压力预测方法、装置和电子设备
CN113176618A (zh) * 2021-03-23 2021-07-27 中国地质大学(武汉) 一种深部地温场预测方法、设备及存储介质

Also Published As

Publication number Publication date
CN114895364A (zh) 2022-08-12
US12158557B2 (en) 2024-12-03
WO2023202047A1 (zh) 2023-10-26
US20240264332A1 (en) 2024-08-08

Similar Documents

Publication Publication Date Title
He et al. Three-dimensional reservoir description from multiwell pressure data and prior information
CN113176618B (zh) 一种深部地温场预测方法、设备及存储介质
Cui et al. Fracture diagnosis in multiple-stage-stimulated horizontal well by temperature measurements with fast marching method
US10775531B2 (en) Big data point and vector model
WO2021146409A1 (en) Generation of a virtual three-dimensional model of a hydrocarbon reservoir
Ribeiro et al. Detecting fracture growth out of zone by use of temperature analysis
Gherabati et al. A large scale network model to obtain interwell formation characteristics
Hazlett et al. Optimal well placement in heterogeneous reservoirs through semianalytic modeling
Chen et al. Fracture inference and optimal well placement using a multiscale history matching in a HPHT tight gas reservoir, Tarim Basin, China
CN114895364B (zh) 基于温压耦合电阻率约束的深部地温场预测方法及装置
CN114895365B (zh) 高温高压岩心电阻率标定的深部地温场预测方法及装置
Liu et al. A fractal model for characterizing hydraulic properties of fractured rock mass under mining influence
Xu et al. The information content and integration of distributed-temperature-sensing data for near-wellbore-reservoir characterization
Alan et al. Prediction of production-inflow profile of a well producing single-phase flow of slightly compressible fluid from multilayer systems by temperature and/or pressure transient data
Putra et al. Thermal modeling and heat flow density interpretation of the onshore Northwest Java Basin, Indonesia
CN110062897B (zh) 使用自组织映射来进行的岩石物理场评估
CN116384267A (zh) 致密储层压裂水平井的最终可采储量的确定方法及设备
Xu et al. Reservoir rock and fluid characterization using distributed temperature sensing DTS systems data
CN118708906B (zh) 一种增强型地热系统采热性能的lstm预测方法及系统
Xu et al. Uncertainty reduction in reservoir geostatistical description using distributed temperature sensing (DTS) systems data
CN118223873B (zh) 一种用于油气田气驱注入的油层压力监测方法
JP2903502B2 (ja) 切羽前方地下水予測システム及び方法
Adeboye et al. PERMEABILITY ESTIMATION AND HYDRAULIC ZONE PORE STRUCTURES IDENTIFICATION USING CORE AND WELL LOGS DATA.
Shtabel et al. Numerical Simulation of the Time-domain Electromagnetic Field for Thaw Zone Detection
Zhang Interpretation of Downhole Temperature Measurements for Multistage Fracture Stimulation in Horizontal Wells

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
CB03 Change of inventor or designer information

Inventor after: Hu Xiangyun

Inventor after: Huang Guoshu

Inventor after: Zhou Wenlong

Inventor after: Cai Jianchao

Inventor after: Yang Jian

Inventor after: Ma Huolin

Inventor before: Huang Guoshu

Inventor before: Hu Xiangyun

Inventor before: Cai Jianchao

Inventor before: Yang Jian

Inventor before: Ma Huolin

Inventor before: Zhou Wenlong

CB03 Change of inventor or designer information