CN113266335A - 一种基于最小二乘算法的航空电磁数据系统辨识的优化方法 - Google Patents
一种基于最小二乘算法的航空电磁数据系统辨识的优化方法 Download PDFInfo
- Publication number
- CN113266335A CN113266335A CN202110294070.2A CN202110294070A CN113266335A CN 113266335 A CN113266335 A CN 113266335A CN 202110294070 A CN202110294070 A CN 202110294070A CN 113266335 A CN113266335 A CN 113266335A
- Authority
- CN
- China
- Prior art keywords
- data
- magnetotelluric
- response
- magnetotelluric response
- emission current
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 71
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 15
- 238000005457 optimization Methods 0.000 title claims abstract description 14
- 230000004044 response Effects 0.000 claims abstract description 161
- 230000002035 prolonged effect Effects 0.000 claims abstract description 9
- 239000004020 conductor Substances 0.000 claims abstract description 5
- 238000004364 calculation method Methods 0.000 claims description 14
- 238000005070 sampling Methods 0.000 claims description 9
- FGUUSXIOTUKUDN-IBGZPJMESA-N C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 Chemical compound C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 FGUUSXIOTUKUDN-IBGZPJMESA-N 0.000 claims description 5
- 230000005540 biological transmission Effects 0.000 claims description 5
- 238000004088 simulation Methods 0.000 claims description 4
- GNFTZDOKVXKIBK-UHFFFAOYSA-N 3-(2-methoxyethoxy)benzohydrazide Chemical compound COCCOC1=CC=CC(C(=O)NN)=C1 GNFTZDOKVXKIBK-UHFFFAOYSA-N 0.000 claims description 2
- 238000005516 engineering process Methods 0.000 abstract description 5
- 230000008859 change Effects 0.000 description 5
- 238000001514 detection method Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 3
- 230000001419 dependent effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 101150071746 Pbsn gene Proteins 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005674 electromagnetic induction Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 230000006698 induction Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
Images
Classifications
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B47/00—Survey of boreholes or wells
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B49/00—Testing 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
-
- 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/08—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with magnetic or electric fields produced or modified by objects or geological structures or by detecting devices
- G01V3/087—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with magnetic or electric fields produced or modified by objects or geological structures or by detecting devices the earth magnetic field being modified by the objects or geological structures
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/02—Agriculture; Fishing; Forestry; Mining
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Geology (AREA)
- Mining & Mineral Resources (AREA)
- Environmental & Geological Engineering (AREA)
- Geochemistry & Mineralogy (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Remote Sensing (AREA)
- Fluid Mechanics (AREA)
- Geophysics (AREA)
- Business, Economics & Management (AREA)
- Health & Medical Sciences (AREA)
- Tourism & Hospitality (AREA)
- Animal Husbandry (AREA)
- Economics (AREA)
- General Health & Medical Sciences (AREA)
- Human Resources & Organizations (AREA)
- Marketing (AREA)
- Primary Health Care (AREA)
- Strategic Management (AREA)
- Marine Sciences & Fisheries (AREA)
- General Business, Economics & Management (AREA)
- Agronomy & Crop Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Electromagnetism (AREA)
- Measuring Magnetic Variables (AREA)
Abstract
本发明涉及一种基于最小二乘算法的航空电磁数据系统辨识的优化方法。录入伪随机发射电流和伪随机发射电流激发地下良导体产生的航空电磁数据；分别截断伪随机发射电流和航空电磁数据的第一周期；对航空电磁数据进行系统辨识，得到有限长度的大地电磁响应数据；对辨识得到的大地电磁响应数据和时间序列同时取自然对数，采用最小二乘法拟合的方法，得到大地电磁响应的一阶拟合式的斜率和截距；延长一倍时间序列长度，代入步骤大地电磁响应的一阶拟合式，得到延长时间段的大地电磁响应数据；对其延长时间部分的，大地电磁响应的幅值进行求解。最终得到延长的大地电磁响应。与传统的系统辨识大地电磁响应技术相比，可以获取更准确的辨识结果。
Description
技术领域
本发明属于时间域航空电磁资料处理技术领域，特别涉及一种基于最小二乘算法的航空电磁数据系统辨识的优化方法。
背景技术
以伪随机序列(PRBS)为发射电流，应用系统辨识方法获取大地电磁响应的时间域航空电磁法(ATEM)是一种基于机载方式，建立在电磁感应原理上，高效的地球物理浅层探测方法。适用于大面积地质，地貌的勘探与调查，具有探测效率高，可较好获取浅层地质信息的优点。
航空电磁数据系统辨识方法是以大地为待辨识系统，由实际发射电流对航空电磁数据进行辨识，得到大地电磁响应数据。系统辨识方法的误差主要由两点原因引起，第一个原因为：第一周期发射电流开启时由0A下降至最小幅值，与其他周期发射电流在正负最大幅值间变化不同，破坏了系统辨识所需要的周期性条件，会使辨识结果存在较大误差。第二个原因为：航空电磁数据为实际发射电流与无限长度的大地电磁响应的卷积，而因可辨识长度的限制，只能辨识有限的长度的大地电磁响应数据，则未参与辨识的大地电磁响应带来的影响将计入辨识得到的大地电磁响应中，并因为卷积运算，会随着采集时间长度的增加，出现累积误差，对系统辨识结果精度出现较大影响。因此对系统辨识方法进行优化，获取高精度的大地电磁响应是航空电磁探测中的关键问题。去除第一周期数据的影响，再通过最小二乘的方法对辨识得到的大地电磁响应数据进行拟合，将未辨识的大地电磁响应对航空电磁数据的影响去除，最后通过迭代提高辨识精度，相对原有的系统辨识方法更加贴近实际情况，减少辨识误差，得到正确的大地电磁响应数据。
武欣等(2015)在“伪随机编码源电磁响应的精细辨识”一文中对辨识得到的大地电磁响应与仿真计算得到的大地电磁响应对比，对两者误差计算了误差函数，并只截取误差较小的部分进行分析。该方法浪费数据，且不能保证辨识精度。武欣等对于误差较大部分的辨识结果，认为对数据进行适当平滑即可提高精度，缺乏理论依据。武欣等的方法并未涉及通过最小二乘法对系统辨识方法进行优化。
根据长期的观察和资料，大多数应用系统辨识方法并未针对系统辨识中，因辨识长度有限带来的误差进行改进。如石琦等在“伪随机编码磁性源瞬变电磁发射技术及电磁响应分析”一文中分别对不含噪声的航空电磁数据和含噪声的电磁数据进行系统辨识，因含噪声的系统辨识存在更大的误差，则忽略了对不含噪声的系统辨识在尾端误差的分析与处理。本发明提出了一种基于最小二乘算法的航空电磁系统辨识优化方法，到目前为止，尚未发现关于应用最小二乘法对系统辨识得到的大地电磁响应进行拟合，并在航空电磁数据中去除未参与辨识的大地响应的影响，来提高系统辨识精度方法的报道。
发明内容
本发明所要解决的技术问题在于提供一种基于最小二乘算法的航空电磁数据系统辨识的优化方法，能够精准的获取大地电磁响应数据，解决原有系统辨识方法因辨识矩阵长度有限使得辨识结果存在误差的问题，并大幅度提高辨识结果精度。
本发明是这样实现的，
一种基于最小二乘算法的航空电磁数据系统辨识的优化方法，该方法包括：
步骤a、录入伪随机发射电流和伪随机发射电流激发地下良导体产生的航空电磁数据；
步骤b、分别截断伪随机发射电流和航空电磁数据的第一周期；并通过数据叠加的方法抑制噪声；
步骤c、对航空电磁数据进行系统辨识，得到有限长度的大地电磁响应数据；
步骤d、对辨识得到的大地电磁响应数据和时间序列同时取自然对数，采用最小二乘法拟合的方法，得到大地电磁响应的一阶拟合式的斜率和截距；
步骤e、延长一倍时间序列长度，代入步骤大地电磁响应的一阶拟合式，得到延长时间段的大地电磁响应数据；
步骤f、将延长时间段的大地电磁响应对应的航空电磁数据视为录入航空电磁数据中的偏差值，并去除；
步骤g、对去偏差后的航空电磁数据进行系统辨识，得到有限长度的大地电磁响应数据；
步骤h、判断去偏差前后的大地电磁响应差异是否达到目标精度，否，回到步骤c，是，则输出大地电磁响应数据。
进一步地，步骤a由r级反馈式线性移位寄存器生成幅值为±Amp变化的 m序列伪随机发射电流，当航空电磁系统发射机的发射频率为fc，航空电磁系统接收机的采样频率为fs，产生最大二进制序列长度为：
p周期的伪随机发射电流表示为i＝[0i(1)i(2)…i(pN)]；
数值模拟的航空电磁数据为正演计算得到的理想大地模型径向磁场的阶跃响应与差分发射电流的卷积，表示为
其中
为理想大地模型径向磁场的阶跃响应数据，
为差分伪随机发射电流数据，
为航空电磁数据。
进一步地，步骤b具体包括：
根据公式
需先将伪随机发射电流取差分，再将差分后的伪随机电流数据截掉第一周期，表示为
将伪随机发射电流生成的航空电磁数据截掉第一周期，表示为 d＝[d(1) d(2) … d(N) d(N+1) d(N+2) … d((p-1)N)]，实际的大地电磁响应数据表示为 g＝[g(1) g(2) g(3) …g(∞)]。
进一步地，将去掉第一周期数据后的差分伪随机发射电流数据和航空电磁数据叠加为单周期数据，分别表示为
差分伪随机发射电流数据序列长度与单周期伪随机序列长度相同，为 I(t)＝[I(1) I(2) … I(N)]，航空电磁数据表示为D(t)＝[D(1) D(2) … D(N)]。
进一步地，所述步骤c具体包括：航空电磁数据采用公式(4)表示，
D(t)＝-I(t)*g(t) (4)
对于公式(4)，基于Wiener-Hopf方程，将方程两端与I(t)互相关，得到：
CR(D(t),I(t))＝-AR(I(t))*g(t) (5)
其中CR(D(t),I(t))表示信号D(t)与I(t)的互相关，AR(I(t))表示信号I(t)的自相关；
采用循环相关的方法计算CR(D(t),I(t))和AR(I(t))，其计算方法为公式(6)：
其中Ic(t+τ)为长度为t的I(t)首尾相连的信号，此时大地电磁数据最长的可辨识的长度为单周期数据长度减一个不可辨识的码元长度，为，
式中Ng为最多的可辨识大地电磁响应的采样点数，系统辨识得到的有限长度的大地电磁响应表示为
辨识时间序列为
进一步地，步骤d包括：
一阶最小二乘法的拟合式为
y＝kx+b (8)
式中，k为拟合式的斜率，b为拟合式的截距，将系统辨识得到的大地电磁响应数据和时间序列数据同时取对数，并代入拟合式，得其近似表达式
log(G(t))＝k log(t)+b (9)
则拟合后的航空电磁数据近似为
G(t)＝tk·eb (10)。
进一步地，步骤e，延长一倍时间序列长度，带入大地电磁响应的一阶拟合式，得到延长时间段的大地电磁响应数据包括：
将时间长度t延长一倍，时间序列表示为
通过拟合后的航空电磁数据近似表达式(10)，拟合得到的大地电磁响为：
进一步地，步骤f包括：
由公式(4)可得，伪随机发射电流与大地电磁响应的卷积为航空电磁数据，航空电磁数据表示为，
D(t)＝-I(t)*g(t)＝-I(t)*G(t)＝-I(t)*Gp(t) (11)
上式对于拟合得到的大地电磁响应得到公式(12)，
将拟合得到的大地电磁响应Gp分为待辨识大地电磁响应
和偏差值
两部分，写为
式中第一项为去偏差的大地电磁响应与伪随机发射电流的卷积，第二项为偏差值与伪随机发射电流的卷积。
进一步地，步骤g，去偏差的航空电磁数据为:
将去偏差的航空电磁数据D'(t)代入公式(5)，可辨识得到第一次去偏差后辨识得到的大地电磁响应数据G'。
进一步地，步骤h，判断去偏差前后的大地电磁响应差异是否达到目标精度采用公式(15)得到：
若去偏差前后的大地电磁响应未达到目标精度，则将去偏差后的大地电磁响应重新应用最小二乘法进行拟合，并依据公式(14)，在第一次去偏差后的航空数据中重新去掉此次偏差值得影响；若去偏差前后的大地电磁响应达到目标精度，则输出基于最小二乘算法的航空电磁系统辨识优化方法辨识得到的高精度的大地电磁响应。
本发明与现有技术相比，有益效果在于：本发明与传统的系统辨识大地电磁响应技术相比，可以获取更准确的辨识结果。与现有通过硬件提高发射序列阶数，延长发射时间的技术相比，本发明更经济，简便，可用低阶伪随机序列达到高阶伪随机序列的探测效果。与通过舍弃或滤波晚期数据的技术相比，本发明更有效的利用数据，可获取更多的大地信息，优化了现有的系统辨识技术，辨识得到更准确大地电磁响应，提高了数据的使用率。。
附图说明
图1是基于最小二乘算法的航空电磁相关辨识优化方法流程图；
图2a是本发明实例提供的伪随机发射电流图；
图2b是本发明实例提供的航空电磁数据图；
图3a是本发明实例提供的系统辨识得到的大地电磁响应数据图(虚线条) 和实际大地电磁响应数据图(实线条)；
图3b是本发明实例提供的系统辨识得到的大地电磁响应数据与实际大地电磁响应数据的相对误差图；
图4是本发明实例提供的经最小二乘算法拟合后的大地电磁响应数据图；
图5a是本发明实例提供的第一次去偏差后系统辨识得到的大地电磁响应数据图(虚线条)和正演计算得到的理想大地电磁响应数据图(实线条)；
图5b是本发明实例提供的第一次去偏差后辨识得到的大地电磁响应数据与实际大地电磁响应数据的相对误差图；
图6a是本发明实例提供的十次迭代后辨识得到的大地电磁响应数据分别与实际大地电磁响应数据的相对误差图；
图6b是本发明实例提供的多次迭代去偏差前后的大地电磁响应差异的相对误差图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
参见图1所示，一种基于最小二乘算法的航空电磁系统辨识优化方法，该方法包括：
a、录入伪随机发射电流和伪随机发射电流激发地下良导体产生的航空电磁数据(即dB/dt)；
b、为保持数据的周期性，分别截断伪随机发射电流和航空电磁数据的第一周期；并通过数据叠加的方法抑制噪声；
c、对航空电磁数据进行系统辨识，得到有限长度的大地电磁响应数据；
d、对辨识得到的大地电磁响应数据和时间序列同时取自然对数，采用最小二乘法拟合的方法，得到其一阶拟合式的斜率和截距；
e、延长一倍时间序列长度，代入大地电磁响应的一阶拟合式，得到延长时间段的大地电磁响应数据；
f、将延长时间段大地电磁响应对应的航空电磁数据视为录入航空电磁数据中的偏差值，并去除；
g、对去偏差后的航空电磁数据进行系统辨识，得到有限长度的大地电磁响应数据；
h、判断去偏差前后的大地电磁响应差异是否达到目标精度，否，回到步骤c，是，则输出大地电磁响应数据；
步骤a，录入伪随机发射电流和伪随机发射电流激发地下良导体产生的航空电磁数据(即dB/dt)；
由r级反馈式线性移位寄存器生成幅值为±Amp变化的m序列伪随机发射电流，当航空电磁系统发射机的发射频率为fc，航空电磁系统接收机的采样频率为fs，可产生最大二进制序列长度为
p周期的伪随机发射电流表示为i＝[0 i(1) i(2) … i(pN)]。
数值模拟的航空电磁数据为正演计算得到的理想大地模型径向磁场的阶跃响应与差分发射电流的卷积，可表示为
其中
为理想大地模型径向磁场的阶跃响应数据，即本发明待辨识的大地电磁响应数据，
为差分伪随机发射电流数据，
为航空电磁数据。
步骤b，为保持数据的周期性，分别截断伪随机发射电流和航空电磁数据的第一周期；并通过数据叠加的方法抑制噪声；
伪随机发射电流和航空电磁数据的周期性是影响系统辨识精度的因素之一，第一周期发射电流由0A下降至-Amp，与后续周期发射电流在Amp与-Amp 之间交替变化的规律性不同，则截掉第一周期数据，由第二周期开始录入。根据公式(2)，需先将伪随机电流取差分，再将差分后的伪随机电流数据截掉第一周期，可表示为
将伪随机发射电流生成的航空电磁数据截掉第一周期，可表示为 d＝[d(1) d(2) … d(N) d(N+1) d(N+2) … d((p-1)N)]，实际的大地电磁响应数据可表示为 g＝[g(1) g(2) g(3) … g(∞)]。
考虑到实际测量中存在随机噪声等影响，将去掉第一周期数据后的多周期发射电流数据和航空电磁数据叠加为单周期数据，分别表示为
此时差分伪随机发射电流数据序列长度与单周期伪随机序列长度相同，为 I(t)＝[I(1) I(2) … I(N)]，航空电磁数据可表示为D(t)＝[D(1) D(2) … D(N)]。步骤c，伪随机发射电流对航空电磁数据进行系统辨识，得到有限长度的大地电磁响应数据。
伪随机发射电流对航空电磁数据进行系统辨识的基本原理基于 Wiener-Hopf方程。在不考虑噪声和其他因素的情况下，航空电磁探测可表示为一线性时不变系统，
D(t)＝-I(t)*g(t) (4)
对于公式(4)，基于Wiener-Hopf方程，将方程两端与I(t)互相关，得到：
CR(D(t),I(t))＝-AR(I(t))*g(t) (5)
其中CR(D(t),I(t))表示信号D(t)与I(t)的互相关，AR(I(t))表示信号I(t)的自相关。
采用循环相关的方法计算CR(D(t),I(t))和AR(I(t))，可以有效地避免线性相关移位计算过程中，因对序列的头部和尾部补零而产生的具有震荡特性的旁瓣对结果的影响。循环相关计算方法是将序列首位环状相连，避免补零，保持了数据的周期性，其计算方法为
其中Ic(t+τ)为长度为t的I(t)首尾相连的信号。此时大地电磁数据最长的可辨识的长度为单周期数据长度减一个不可辨识的码元长度，即表示为，
式中Ng为最多的可辨识大地电磁响应的采样点数。系统辨识得到的有限长度的大地电磁响应表示为
辨识时间序列可表示为
由于现实情况下实际的大地电磁响应数据是无限长度的，但计算过程中因最长可辨识长度的限制，无法辨识无限长度的大地电磁响应数据，则没有参与计算的大地电磁响应数据会算入辨识得到的大地电磁响应数据中，使数值结果偏小。因此需要从原航空电磁数据中去除未参与计算的大地电磁的带来的影响。
步骤d，对辨识得到的大地电磁响应数据和时间序列同时取自然对数，采用最小二乘法拟合的方法，得到其一阶拟合式的斜率和截距；
一阶最小二乘法的拟合式为
y＝kx+b (8)
式中，k为拟合式的斜率，b为拟合式的截距。将系统辨识得到的大地电磁响应数据和时间序列数据同时取对数，并代入拟合式，可得其近似表达式
log(G(t))＝k log(t)+b (9)
则拟合后的航空电磁数据可近似表示为
G(t)＝tk·eb (10)
步骤e，延长一倍时间序列长度，带入大地电磁响应的一阶拟合式，得到延长时间段的大地电磁响应数据；
将时间长度t延长一倍，此时时间序列可表示为
通过公式(10)，拟合得到的大地电磁响可表示为
并在后续部分参与运算。
步骤f，将大地电磁响应数据分为辨识部分和偏差值，并去除航空电磁数据中偏差值的影响；
由公式(4)可得，伪随机发射电流与大地电磁响应的卷积为航空电磁数据。航空电磁数据可表示为，
D(t)＝-I(t)*g(t)＝-I(t)*G(t)＝-I(t)*Gp(t) (11)
由上式可得，当航空数据一定时，因大地电磁响应辨识的长度不同，会影响辨识得到的大地电磁响应的结果。则上式对于拟合得到的大地电磁响应，还可以写为，
m，n分别为自变量和因变量。自变量m为I的长度，m∈(1,2Ng)，因变量 n为延长后的大地电磁响应Gp的长度，Gp∈(1,2Ng)。
将拟合得到的大地电磁响应Gp分为待辨识大地电磁响应
和偏差值
两部分，公式(12)可写为
式中第一项为去偏差的大地电磁响应与伪随机发射电流的卷积。第二项为偏差值与伪随机发射电流的卷积，需要从航空电磁数据中去除。
步骤g，伪随机发射电流对去偏差后的航空电磁数据进行系统辨识。
去偏差的航空电磁数据可表示为：
将去偏差的航空电磁数据D'(t)代入公式(5)，可辨识得到第一次去偏差后辨识得到的大地电磁响应数据G'。
步骤h，判断去偏差前后的大地电磁响应差异是否达到目标精度，否，回到步骤c，将去偏差后的航空电磁数据重新拟合，并辨识；
设置去偏差前后大地电磁响应差异的精度，可表示为，
若去偏差前后的大地电磁响应未达到目标精度，则将去偏差后的大地电磁响应重新应用最小二乘法进行拟合，并依据公式(14)，在第一次去偏差后的航空数据中重新去掉此次偏差值得影响。若去偏差前后的大地电磁响应达到目标精度，则输出基于最小二乘算法的航空电磁系统辨识优化方法辨识得到的高精度的大地电磁响应。
实施例
以一组4周期，9阶，幅值为5A的伪随机编码数据作为发射电流，通过电导率为0.01S/m的均匀大地模型，该理想大地模型径向磁场磁感应强度的变化率(dB/dt)为航空电磁数据为例。设航空电磁系统发射机的发射频率为200kHz，航空电磁系统接收机的采样频率为1MHz，录入的伪随机发射电流如图2a所示，录入航空电磁数据如图2b所示。
系统辨识需要保持数据的周期性，第一周期发射电流由0A下降至-5A，与后续周期电流由5A至-5A变化的规律性不同，则截掉第一周期数据，由第二周期开始录入。截掉的伪随机发射电流数据单周期序列长度为
对录入4周期的伪随机编码电流进行差分，并截掉其第一周期数据，可表示为
对录入航空电磁数据截掉第一周期数据，可表示为d＝(d(1) d(2) … d(7665))。为减少随机噪声的干扰，对现有三周期伪随机发射电流和航空电磁数据进行叠加，叠加后的单周期伪随机发射电流数据和航空电磁数据可表示为，
此时差分伪随机发射电流数据序列可表示为I(t)＝[I(1) I(2) … I(2555)]，航空电磁数据可表示为D(t)＝[D(1) D(2) … D(2555)]，实际的大地模型径向磁场的阶跃响应数据可表示为g＝[g(1) g(2) g(3) … g(∞)]。
航空电磁探测为线性时不变系统，将伪随机发射电流和航空电磁数据代入不考虑噪声和其他因素的Wiener-Hopf方程，并将将方程两端与I(t)互相关，得到：
CR(D(t),I(t))＝-AR(I(t))*g(t) (5)
若将公式(5)简化为A(t)＝CR(D(t),I(t))，B(t)＝AR(I(t))。公式(5)在离散系统可写为
Ng为可辨识大地电磁响应的采样点数，为2550点。A和B矩阵通过循环相关方法求得。根据公式(4)，将公式(3)在离散系统矩阵化，表示为为A＝BG，其中
式中，n1和n2分别为A(ni)和B(nj)序列中数据的最大值采样点的序列号。此时辨识得到的大地电磁响应是一个长度为Ng的列向量，G＝[G(1) G(2) G(3) … G(2550)]，其时间序列为t＝[t(1) t(2) t(3) … t(2550)]。系统辨识得到的大地电磁响应图如图 3a实线所示，因卷积运算的本质为累加计算，随着时间序列的增加会逐渐产生较大的累积误差，与实际的大地电磁响应数据相比，如图3a实线所示，在晚期有较大的误差。系统辨识得到的大地电磁响应数据与实际大地电磁响应数据的相对误差如图3b所示，在1.3ms时，其相对误差将达到100％。则需要对系统辨识进一步优化，减小辨识误差。
采用最小二乘法拟合辨识得到的大地电磁响应数据，将辨识得到的大地电磁响应数据和时间序列数据同时取对数，并代入一阶拟合式，计算可得斜率为 -2.89，截距为-50.65。则辨识得到的大地电磁响应数据和表示为
G(t)＝t-2.89·e-50.65 (10)
将时间序列t延长一倍，此时时间序列可表示为 T＝[T(1) T(2) … T(5100)]，通过公式(10)拟合得到延长后的大地电磁响应可表示为
如图4所示虚线左侧为辨识得到的大地电磁响应数据，虚线右侧为基于最小二乘拟合后的大地电磁响应数据。
航空电磁数据D可由延长后的大地电磁响应Gp与差分后的发射电流I表示。并将大地电磁响应Gp分为待辨识大地电磁响应
G1＝[G1(1) G1(2) G1(3) … G1(2555)]和偏差值G2＝[G2(1) G2(2) G2(3) … G2(2555)]两部分，可表示为
将大地电磁响应的偏差值感应的航空电磁数据从原航空电磁数据中去除，再经过系统辨识，可得第一次去偏差后辨识得到的大地电磁响应数据，如图5a 虚线所示，在晚期的误差与图3相比大幅度减小。第一次去偏差后系统辨识得到的大地电磁响应数据与实际大地电磁响应数据的相对误差如图5b所示，在 1.3ms时，其相对误差已下降至70％。可证明基于最小二乘算法的航空电磁系统辨识优化方法可以大幅度减小辨识误差。
每次迭代后辨识得到的大地电磁响应数据定义为Gn，Gn与实际的大地电磁响应的相对误差如图6a所示，经过多次迭代后，其去偏差效果逐渐减弱，在 1.3ms时辨识得到的大地电磁响应与实际大地电磁响应之间的相对误差逐渐收敛至15％。
实际情况下不会有实际大地电磁响应数据作为对比，则对比每次辨识前后的的大地电磁响应的相对误差作为精度，
如图6b所示，第一点即第一次迭代结果与第二次迭代结果的相对误差最大，可得第二次迭代效果做好，之后迭代结果去偏差效果较小，以经济的角度可以第二次结果作为最佳补偿结果。若以精度作为停止迭代标准，如当精度为 1％时，则第七次迭代后可输出所需结果。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
Claims (10)
1.一种基于最小二乘算法的航空电磁数据系统辨识的优化方法，其特征在于，该方法包括：
步骤a、录入伪随机发射电流和伪随机发射电流激发地下良导体产生的航空电磁数据；
步骤b、分别截断伪随机发射电流和航空电磁数据的第一周期；并通过数据叠加的方法抑制噪声；
步骤c、对航空电磁数据进行系统辨识，得到有限长度的大地电磁响应数据；
步骤d、对辨识得到的大地电磁响应数据和时间序列同时取自然对数，采用最小二乘法拟合的方法，得到大地电磁响应的一阶拟合式的斜率和截距；
步骤e、延长一倍时间序列长度，代入步骤大地电磁响应的一阶拟合式，得到延长时间段的大地电磁响应数据；
步骤f、将延长时间段的大地电磁响应对应的航空电磁数据视为录入航空电磁数据中的偏差值，并去除；
步骤g、对去偏差后的航空电磁数据进行系统辨识，得到有限长度的大地电磁响应数据；
步骤h、判断去偏差前后的大地电磁响应差异是否达到目标精度，否，回到步骤c，是，则输出大地电磁响应数据。
2.按照权利要求1所述的方法，其特征在于，步骤a由r级反馈式线性移位寄存器生成幅值为±Amp变化的m序列伪随机发射电流，当航空电磁系统发射机的发射频率为fc，航空电磁系统接收机的采样频率为fs，产生最大二进制序列长度为：
p周期的伪随机发射电流表示为i＝[0 i(1) i(2) … i(pN)]；
数值模拟的航空电磁数据为正演计算得到的理想大地模型径向磁场的阶跃响应与差分发射电流的卷积，表示为
其中
为理想大地模型径向磁场的阶跃响应数据，
为差分伪随机发射电流数据，
为航空电磁数据。
3.按照权利要求2所述的方法，其特征在于，步骤b具体包括：
根据公式
需先将伪随机发射电流取差分，再将差分后的伪随机电流数据截掉第一周期，表示为
将伪随机发射电流生成的航空电磁数据截掉第一周期，表示为d＝[d(1) d(2) … d(N) d(N+1) d(N+2) … d((p-1)N)]，实际的大地电磁响应数据表示为g＝[g(1) g(2) g(3) …g(∞)]。
4.按照权利要求3所述的方法，其特征在于，将去掉第一周期数据后的差分伪随机发射电流数据和航空电磁数据叠加为单周期数据，分别表示为
差分伪随机发射电流数据序列长度与单周期伪随机序列长度相同，为I(t)＝[I(1) I(2) … I(N)]，航空电磁数据表示为D(t)＝[D(1) D(2) … D(N)]。
5.按照权利要求1所述的方法，其特征在于，所述步骤c具体包括：航空电磁数据采用公式(4)表示，
D(t)＝-I(t)*g(t) (4)
对于公式(4)，基于Wiener-Hopf方程，将方程两端与I(t)互相关，得到：
CR(D(t),I(t))＝-AR(I(t))*g(t) (5)
其中CR(D(t),I(t))表示信号D(t)与I(t)的互相关，AR(I(t))表示信号I(t)的自相关；
采用循环相关的方法计算CR(D(t),I(t))和AR(I(t))，其计算方法为公式(6)：
其中Ic(t+τ)为长度为t的I(t)首尾相连的信号，此时大地电磁数据最长的可辨识的长度为单周期数据长度减一个不可辨识的码元长度，为，
式中Ng为最多的可辨识大地电磁响应的采样点数，系统辨识得到的有限长度的大地电磁响应表示为
辨识时间序列为
6.按照权利要求1所述的方法，其特征在于，步骤d包括：
一阶最小二乘法的拟合式为
y＝kx+b (8)
式中，k为拟合式的斜率，b为拟合式的截距，将系统辨识得到的大地电磁响应数据和时间序列数据同时取对数，并代入拟合式，得其近似表达式
log(G(t))＝klog(t)+b (9)
则拟合后的航空电磁数据近似为
G(t)＝tk·eb (10)。
7.按照权利要求6所述的方法，其特征在于，步骤e，延长一倍时间序列长度，带入大地电磁响应的一阶拟合式，得到延长时间段的大地电磁响应数据包括：
将时间长度t延长一倍，时间序列表示为
通过拟合后的航空电磁数据近似表达式(10)，拟合得到的大地电磁响为：
8.按照权利要求7所述的方法，其特征在于，步骤f包括：
由公式(4)可得，伪随机发射电流与大地电磁响应的卷积为航空电磁数据，航空电磁数据表示为，
D(t)＝-I(t)*g(t)＝-I(t)*G(t)＝-I(t)*Gp(t) (11)
上式对于拟合得到的大地电磁响应得到公式(12)，
将拟合得到的大地电磁响应Gp分为待辨识大地电磁响应
和偏差值
两部分，写为
式中第一项为去偏差的大地电磁响应与伪随机发射电流的卷积，第二项为偏差值与伪随机发射电流的卷积。
9.按照权利要求7所述的方法，其特征在于，步骤g，去偏差的航空电磁数据为:
将去偏差的航空电磁数据D'(t)代入公式(5)，可辨识得到第一次去偏差后辨识得到的大地电磁响应数据G'。
10.按照权利要求9所述的方法，其特征在于，步骤h，判断去偏差前后的大地电磁响应差异是否达到目标精度采用公式(15)得到：
若去偏差前后的大地电磁响应未达到目标精度，则将去偏差后的大地电磁响应重新应用最小二乘法进行拟合，并依据公式(14)，在第一次去偏差后的航空数据中重新去掉此次偏差值得影响；若去偏差前后的大地电磁响应达到目标精度，则输出基于最小二乘算法的航空电磁系统辨识优化方法辨识得到的高精度的大地电磁响应。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110294070.2A CN113266335B (zh) | 2021-03-19 | 2021-03-19 | 一种基于最小二乘算法的航空电磁数据系统辨识的优化方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110294070.2A CN113266335B (zh) | 2021-03-19 | 2021-03-19 | 一种基于最小二乘算法的航空电磁数据系统辨识的优化方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113266335A true CN113266335A (zh) | 2021-08-17 |
| CN113266335B CN113266335B (zh) | 2022-06-21 |
Family
ID=77228326
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110294070.2A Active CN113266335B (zh) | 2021-03-19 | 2021-03-19 | 一种基于最小二乘算法的航空电磁数据系统辨识的优化方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113266335B (zh) |
Citations (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| GB945365A (en) * | 1960-01-06 | 1963-12-23 | Barringer Research Ltd | Method and apparatus for the remote detection of conducting bodies |
| US20070021916A1 (en) * | 2003-06-10 | 2007-01-25 | Ohm Limited | Electromagnetic surveying for hydrocarbon reservoirs |
| WO2015026650A1 (en) * | 2013-08-20 | 2015-02-26 | Technoimaging, Llc | Systems and methods for remote electromagnetic exploration for mineral and energy resources using stationary long-range transmitters |
| WO2016010915A1 (en) * | 2014-07-12 | 2016-01-21 | Halliburton Energy Services, Inc. | Detecting defects in non-nested tubings and casings using calibrated data and time thresholds |
| CN105652325A (zh) * | 2016-01-05 | 2016-06-08 | 吉林大学 | 基于指数拟合-自适应卡尔曼的地空电磁数据去噪方法 |
| US20160195631A1 (en) * | 2015-01-07 | 2016-07-07 | The Regents Of The University Of California | System and method for groundwater detection and evaluation |
| US20160282498A1 (en) * | 2015-03-27 | 2016-09-29 | Cgg Services Sa | Apparatus and method for calculating earth's polarization properties from airborne time-domain electromagnetic data |
| US20170068014A1 (en) * | 2014-02-28 | 2017-03-09 | Cgg Services Sa | Systems and methods for a composite magnetic field sensor for airborne geophysical surveys |
| CN107121706A (zh) * | 2017-05-08 | 2017-09-01 | 厦门大学 | 基于波恩迭代法的航空瞬变电磁电导率三维反演方法 |
| CN110058317A (zh) * | 2019-05-10 | 2019-07-26 | 成都理工大学 | 航空瞬变电磁数据和航空大地电磁数据联合反演方法 |
| CN110398778A (zh) * | 2019-07-10 | 2019-11-01 | 吉林大学 | 一种基于等效采样的航空电磁浅层数据相关辨识方法 |
| US20210071640A1 (en) * | 2019-01-03 | 2021-03-11 | Lucomm Technologies, Inc. | Generative System |
-
2021
- 2021-03-19 CN CN202110294070.2A patent/CN113266335B/zh active Active
Patent Citations (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| GB945365A (en) * | 1960-01-06 | 1963-12-23 | Barringer Research Ltd | Method and apparatus for the remote detection of conducting bodies |
| US20070021916A1 (en) * | 2003-06-10 | 2007-01-25 | Ohm Limited | Electromagnetic surveying for hydrocarbon reservoirs |
| WO2015026650A1 (en) * | 2013-08-20 | 2015-02-26 | Technoimaging, Llc | Systems and methods for remote electromagnetic exploration for mineral and energy resources using stationary long-range transmitters |
| US20170068014A1 (en) * | 2014-02-28 | 2017-03-09 | Cgg Services Sa | Systems and methods for a composite magnetic field sensor for airborne geophysical surveys |
| WO2016010915A1 (en) * | 2014-07-12 | 2016-01-21 | Halliburton Energy Services, Inc. | Detecting defects in non-nested tubings and casings using calibrated data and time thresholds |
| US20160195631A1 (en) * | 2015-01-07 | 2016-07-07 | The Regents Of The University Of California | System and method for groundwater detection and evaluation |
| US20160282498A1 (en) * | 2015-03-27 | 2016-09-29 | Cgg Services Sa | Apparatus and method for calculating earth's polarization properties from airborne time-domain electromagnetic data |
| CN105652325A (zh) * | 2016-01-05 | 2016-06-08 | 吉林大学 | 基于指数拟合-自适应卡尔曼的地空电磁数据去噪方法 |
| CN107121706A (zh) * | 2017-05-08 | 2017-09-01 | 厦门大学 | 基于波恩迭代法的航空瞬变电磁电导率三维反演方法 |
| US20210071640A1 (en) * | 2019-01-03 | 2021-03-11 | Lucomm Technologies, Inc. | Generative System |
| CN110058317A (zh) * | 2019-05-10 | 2019-07-26 | 成都理工大学 | 航空瞬变电磁数据和航空大地电磁数据联合反演方法 |
| CN110398778A (zh) * | 2019-07-10 | 2019-11-01 | 吉林大学 | 一种基于等效采样的航空电磁浅层数据相关辨识方法 |
Non-Patent Citations (12)
Also Published As
| Publication number | Publication date |
|---|---|
| CN113266335B (zh) | 2022-06-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106990435B (zh) | 一种减弱依赖初至拾取精度的微震定位方法及装置 | |
| CN103941254A (zh) | 一种基于地质雷达的土壤物性类别识别方法和装置 | |
| CN107229033B (zh) | 基于高度维分段搜索的多目标到达时间差定位方法 | |
| CN105137498A (zh) | 一种基于特征融合的地下目标探测识别系统及方法 | |
| CN104502980A (zh) | 一种电磁大地冲激响应的辨识方法 | |
| Hua et al. | Optimal VMD-based signal denoising for laser radar via Hausdorff distance and wavelet transform | |
| Li | A new robust signal recognition approach based on holder cloud features under varying SNR environment | |
| Wanchun et al. | Location algorithms for moving target in non‐coherent distributed multiple‐input multiple‐output radar systems | |
| CN113866736B (zh) | 激光雷达回波信号高斯拐点选择分解方法 | |
| CN113266335B (zh) | 一种基于最小二乘算法的航空电磁数据系统辨识的优化方法 | |
| Li et al. | Identification and parameter estimation algorithm of radar signal subtle features | |
| CA2412995A1 (en) | Seismic survey system | |
| CN107918152B (zh) | 一种地震相干层析成像方法 | |
| CN115219991A (zh) | 一种基于希尔伯特变换的二相编码调制信号识别方法 | |
| Lin et al. | Signal generation and continuous tracking with signal attribute variations using software simulation | |
| Yang et al. | LPI radar signal detection based on the combination of FFT and segmented autocorrelation plus PAHT | |
| Soliman | The accuracy of the Gaussian tail and Dual Dirac model in jitter histogram and probability density functions | |
| Cheng et al. | A radar main lobe pulse correlation sorting method | |
| CN110726995B (zh) | 激光雷达高精度测距方法及系统 | |
| Li et al. | Pulse train detection algorithm based on edge enhancement | |
| Zhu et al. | A Multistep Method for Automatic Determination and Optimization of Microseismic P‐Phase Arrival Times in a Coal Mine | |
| CN110865351A (zh) | 一种高速机动目标参数估计方法 | |
| Aliyu et al. | Airwaves Detection and Elimination Using Fast Fourier Transform to Enhance Detection of Hydrocarbon | |
| Liu et al. | Adaptive Parameter Estimation for Compound Interrupted Sampling and Repeating Jamming | |
| CN111914802B (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 |