CN109031419A - A kind of method and system for picking up microseism first arrival - Google Patents

A kind of method and system for picking up microseism first arrival Download PDF

Info

Publication number
CN109031419A
CN109031419A CN201810841073.1A CN201810841073A CN109031419A CN 109031419 A CN109031419 A CN 109031419A CN 201810841073 A CN201810841073 A CN 201810841073A CN 109031419 A CN109031419 A CN 109031419A
Authority
CN
China
Prior art keywords
shearing wave
useful signal
microseism
coefficient
wave coefficient
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.)
Pending
Application number
CN201810841073.1A
Other languages
Chinese (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.)
Yangtze University
Original Assignee
Yangtze 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 Yangtze University filed Critical Yangtze University
Priority to CN201810841073.1A priority Critical patent/CN109031419A/en
Publication of CN109031419A publication Critical patent/CN109031419A/en
Pending legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • G01V1/306Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/62Physical property of subsurface
    • G01V2210/624Reservoir parameters

Landscapes

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

Abstract

The invention discloses a kind of method and system for picking up microseism first arrival, are applied to micro-seismic monitoring.The method include that the micro-seismic monitoring signal to acquisition carries out shearing wave conversion;Time frequency analysis is carried out to microseism data monitoring data, the dominant frequency of microseism useful signal is determined, determines the shearing wave coefficient in useful signal direction corresponding to durection component under different scale;The correlation for calculating the shearing wave coefficient of the useful signal corresponding direction determines the corresponding shearing wave coefficient in the useful signal direction and the useful signal direction according to maximum correlation;According to the corresponding shearing wave coefficient in the useful signal direction, seeks STA/LTA and pick up curve, determine microseism signal first arrival time.The accuracy of microseism first break pickup is greatly improved, and accelerates pickup velocity, and then ensures the real-time of first break picking.

Description

A kind of method and system for picking up microseism first arrival
Technical field
The present invention relates to geology monitoring field more particularly to a kind of method and system for picking up microseism first arrival
Background technique
When subsurface boring, digging the destructively bed rock stone structure such as mine, can all cause microseism, microseism be it is a kind of small-sizedly Shake activity.Microseism can all be monitored by generally carrying out pressure break to underground, to obtain underground state.And pick up microseism P wave first arrival It is of great significance to seismic source location, FRACTURE PREDICTION, seismic source rupture mechanism analysis, thus picks up microseism first arrival extremely have must It wants.
Currently, for the microseism signal monitored, it will usually after handling microseism signal by Wavelet transformation, pass through list One recognition mode first break picking, since microseism signal-to-noise ratio is low, and the bad this mode of wavelet decomposition directionality causes to pick up essence Exactness is lower, and based on after Shearlet conversion process microseism initial signal, it will usually high-order measurement Law is used, in single ruler Degree is picked up by kurtosis Functional Analysis, but the real-time processing speed of this mode is slower, it is difficult to ensure requirement of real-time.
Summary of the invention
The embodiment of the invention provides a kind of method and system for picking up microseism first arrival, can fast and accurately pick up micro- The P wave first arrival of earthquake.
The embodiment of the present invention in a first aspect, provide it is a kind of pick up microseism first arrival method, this method comprises:
S1, acquisition micro-seismic monitoring signal;
S2, shearing wave conversion is carried out to the micro-seismic monitoring signal, obtain the shearing wave of different scale and different directions Transformation coefficient;
S3, time frequency analysis is carried out to microseism data monitoring data, determines the dominant frequency of microseism useful signal, and according to micro- The dominant frequency of earthquake useful signal determines the shearing wave in useful signal direction corresponding to durection component under different scale respectively Coefficient;
S4, the correlation for calculating useful signal shearing wave coefficient in direction corresponding to direction vector under different scale Property, according to the maximum correlation of shearing wave coefficient corresponding to useful signal direction vector under different scale, determine institute State the corresponding shearing wave coefficient in useful signal direction and the useful signal direction;
S5, according to the corresponding shearing wave coefficient in the useful signal direction, seek STA/LTA and pick up curve, described in judgement Picking up the curve acquirement maximum value corresponding time is microseism signal first arrival time.
In the second aspect of the embodiment of the present invention, a kind of system for picking up microseism first arrival is provided, which includes:
Acquisition module: for acquiring micro-seismic monitoring signal;
Conversion module: for carrying out shearing wave conversion to the micro-seismic monitoring signal, different scale and non-Tongfang are obtained To shearing wave transformation coefficient;
Confirmation module: for carrying out time frequency analysis to microseism data monitoring data, the master of microseism useful signal is determined Frequently, and according to the dominant frequency of microseism useful signal, determine the useful signal under different scale corresponding to durection component respectively The shearing wave coefficient in direction;
Computing module: for calculating the shearing wave system in useful signal direction corresponding to direction vector under different scale Several correlation, according to the maximal correlation of shearing wave coefficient corresponding to useful signal direction vector under different scale Property, determine the corresponding shearing wave coefficient in the useful signal direction and the useful signal direction.
Determination module: for seeking STA/LTA and picking up song according to the corresponding shearing wave coefficient in the useful signal direction Line determines that the pickup curve acquirement maximum value corresponding time is microseism signal first arrival time.
As can be seen from the above technical solutions, the embodiment of the present invention has the advantage that
In the embodiment of the present invention, multiple dimensioned, multidirectional signal is obtained by Shearlet transformation, obtains useful signal pair After the Shearlet coefficient answered, useful signal first arrival time is picked up by STA/LTA algorithm.Allow to by Shearlet Transformation obtains multi-direction signal after decomposing, and ensures the extraction of useful signal, and then improve the accuracy picked up, and is based on STA/ LTA can greatly improve pickup velocity, and then ensure the real-time of first break picking.
Detailed description of the invention
It to describe the technical solutions in the embodiments of the present invention more clearly, below will be to embodiment or description of the prior art Needed in attached drawing be briefly described, it should be apparent that, the accompanying drawings in the following description is only of the invention some Embodiment for those of ordinary skill in the art without any creative labor, can also be according to these Attached drawing obtains other attached drawings.
Fig. 1 is the method one embodiment flow chart provided in an embodiment of the present invention for picking up microseism first arrival;
Fig. 2 is the system one embodiment structural schematic diagram provided in an embodiment of the present invention for picking up microseism first arrival.
Specific embodiment
The embodiment of the invention provides a kind of method and system for picking up microseism first arrival, at the beginning of picking up microseism P wave Extremely, it provides and picks up accuracy and real-time.
In order to make the invention's purpose, features and advantages of the invention more obvious and easy to understand, below in conjunction with the present invention Attached drawing in embodiment, technical scheme in the embodiment of the invention is clearly and completely described, it is clear that disclosed below Embodiment be only a part of the embodiment of the present invention, and not all embodiment.Based on the embodiments of the present invention, this field Those of ordinary skill's all other embodiment obtained without making creative work, belongs to protection of the present invention Range.
Referring to Fig. 1, one embodiment flow chart of the method for the pickup microseism first arrival provided in the embodiment of the present invention Include:
S101, acquisition micro-seismic monitoring signal;
The micro-seismic monitoring signal can be collected by monitoring instrument, especially by sensor collection microseism signal, It is simulated again by instrument and generates digital signal or analog signal.
S102, shearing wave conversion is carried out to the micro-seismic monitoring signal, obtain the shearing of different scale and different directions Wave conversion coefficient;
The shearing wave conversion, that is, Shearlet transformation, is generally used for multiple dimensioned signal decomposition, described in the present invention micro- It include the interference signals such as noise in seismic monitoring signal, it is decomposable multiple dimensioned out by Shearlet transformation, different directions Useful signal and noise can guarantee the accuracy of signal extraction in this way.
Specifically, carrying out Shearlet transformation according to formula (1) and formula (2)
Wherein, a ∈ R+, s ∈ R, t ∈ R2, f expression frequency;
Wherein, a,a∈R+, s ∈ R, t ∈ R2, Ma,sFor the matrix that stretches;S is shear parameters, a For scale parameter, t is translation parameters.
S103, time frequency analysis is carried out to microseism data monitoring data, determines the dominant frequency of microseism useful signal, and according to The dominant frequency of microseism useful signal determines the shearing in useful signal direction corresponding to durection component under different scale respectively Wave system number;
The microseism data monitoring data, that is, microseism signal set, history micro-seismic monitoring signal frequency division when passing through Analysis can determine the dominant frequency of useful signal in microseism, according to the dominant frequency of useful signal can determine current microseism signal not With shearing wave coefficient corresponding on scale and different directions.
Optionally, the size of the corresponding shearing wave coefficient is determined according to the frequency of useful signal height.
S104, the correlation for calculating useful signal shearing wave coefficient in direction corresponding to direction vector under different scale Property, according to the maximum correlation of shearing wave coefficient corresponding to useful signal direction vector under different scale, determine institute State the corresponding shearing wave coefficient in useful signal direction and the useful signal direction;
Optionally, step S104 is specifically included:
S1, the correlation that the shearing wave coefficient of corresponding direction described in direction vector under different scale is calculated according to formula (1):
Correlation (k)=U (j, k) * U (j+1, k) (1)
Wherein, U indicates that shearing wave coefficient, j indicate that scale, k indicate direction;
S2, corresponding direction k is maximized as useful signal direction according to Correlation (k), and obtains direction k Corresponding shearing wave coefficient;
S3, the shearing wave coefficient of direction k is added up, the shearing wave coefficient after adding up is the corresponding shearing of the useful signal Wave system number.
S105, according to the corresponding shearing wave coefficient in the useful signal direction, seek STA/LTA and pick up curve, determine institute Stating the pickup curve acquirement maximum value corresponding time is microseism signal first arrival time.
The STA/LTA that seeks picks up curve, can specifically be calculated according to the following formula:
Wherein, x (n) is signal time sequence, WSTAFor short time-window, WLTAWindow width when being long, n represent total sampling of x (n) Points, k represent current corresponding sampled point.
Optionally, when the amplitude for picking up curve reaches maximum value, current maximum corresponding time point is taken.
First arrival time can be determined according to the peak value that useful signal amplitude changes by picking up curve based on STA/LTA, be greatly reduced The high exponent arithmetic(al) bring time delay of traditional high-order measurement Law improves real-time.
It should be understood that the size of the serial number of each step is not meant that the order of the execution order in above-described embodiment, each process Execution sequence should be determined by its function and internal logic, the implementation process without coping with the embodiment of the present invention constitutes any limit It is fixed.
A kind of method for obtaining microseism first arrival is essentially described above, will be to a kind of acquisition microseism first arrival below System is described in detail.
Fig. 2 is one embodiment structural schematic diagram of the system provided in an embodiment of the present invention for obtaining microseism first arrival, should System includes:
Acquisition module 210: for acquiring micro-seismic monitoring signal;
Conversion module 220: for carrying out shearing wave conversion to the micro-seismic monitoring signal, different scale and difference are obtained The shearing wave transformation coefficient in direction;
Confirmation module 230: for carrying out time frequency analysis to microseism data monitoring data, microseism useful signal is determined Dominant frequency, and according to the dominant frequency of microseism useful signal, determine that useful signal durection component institute under different scale is right respectively Answer the shearing wave coefficient in direction;
Optionally, the determination useful signal described respectively direction corresponding to durection component under different scale is cut Cut wave system number specifically:
The size of the corresponding shearing wave coefficient is determined according to the frequency of useful signal height.
Computing module 240: for calculating the shearing in useful signal direction corresponding to direction vector under different scale The correlation of wave system number, according to the maximum phase of shearing wave coefficient corresponding to useful signal direction vector under different scale Guan Xing determines the corresponding shearing wave coefficient in the useful signal direction and the useful signal direction;
Optionally, the computing module 240 includes:
Computing unit: the shearing wave system for the corresponding direction according to direction vector under formula (1) calculating different scale Several correlations:
Correlation (k)=U (j, k) * U (j+1, k) (1)
Wherein, U indicates that shearing wave coefficient, j indicate that scale, k indicate direction;
Confirmation unit: for being maximized corresponding direction k as useful signal direction according to Correlation (k), and Obtain the corresponding shearing wave coefficient of direction k;
Summing elements: for the shearing wave coefficient of direction k to add up, the shearing wave coefficient after adding up is the useful signal Corresponding shearing wave coefficient.
Determination module 250: for seeking STA/LTA pickup according to the corresponding shearing wave coefficient in the useful signal direction Curve determines that the pickup curve acquirement maximum value corresponding time is microseism signal first arrival time.
The system provided in an embodiment of the present invention for obtaining microseism first arrival, energy multi-resolution decomposition micro-seismic monitoring signal, and After determining Shearlet coefficient, seeks STA/LTA pickup curve and determine first arrival time, not only can accurately determine first arrival time, but also It can ensure the real-time picked up.
It is apparent to those skilled in the art that for convenience and simplicity of description, the system of foregoing description, The specific work process of device and unit, can refer to corresponding processes in the foregoing method embodiment, and details are not described herein.
In the above-described embodiments, it all emphasizes particularly on different fields to the description of each embodiment, is not described in detail or remembers in some embodiment The part of load may refer to the associated description of other embodiments.
Those of ordinary skill in the art may be aware that each embodiment described in conjunction with the examples disclosed in this document Module, unit and/or method and step can be realized with the combination of electronic hardware or computer software and electronic hardware.This A little functions are implemented in hardware or software actually, the specific application and design constraint depending on technical solution.Specially Industry technical staff can use different methods to achieve the described function each specific application, but this realization is not It is considered as beyond the scope of this invention.
In several embodiments provided herein, it should be understood that disclosed system, device and method can be with It realizes by another way.For example, the apparatus embodiments described above are merely exemplary, for example, the unit It divides, only a kind of logical function partition, there may be another division manner in actual implementation, such as multiple units or components It can be combined or can be integrated into another system, or some features can be ignored or not executed.Another point, it is shown or The mutual coupling, direct-coupling or communication connection discussed can be through some interfaces, the indirect coupling of device or unit It closes or communicates to connect, can be electrical property, mechanical or other forms.
The unit as illustrated by the separation member may or may not be physically separated, aobvious as unit The component shown may or may not be physical unit, it can and it is in one place, or may be distributed over multiple In network unit.It can select some or all of unit therein according to the actual needs to realize the mesh of this embodiment scheme 's.
It, can also be in addition, the functional units in various embodiments of the present invention may be integrated into one processing unit It is that each unit physically exists alone, can also be integrated in one unit with two or more units.Above-mentioned integrated list Member both can take the form of hardware realization, can also realize in the form of software functional units.
The above, the above embodiments are merely illustrative of the technical solutions of the present invention, rather than its limitations;Although referring to before Stating embodiment, invention is explained in detail, those skilled in the art should understand that: it still can be to preceding Technical solution documented by each embodiment is stated to modify or equivalent replacement of some of the technical features;And these It modifies or replaces, the spirit and scope for technical solution of various embodiments of the present invention that it does not separate the essence of the corresponding technical solution.

Claims (4)

1. a kind of method for picking up microseism first arrival characterized by comprising
S1, acquisition micro-seismic monitoring signal;
S2, shearing wave conversion is carried out to the micro-seismic monitoring signal, obtain the shearing wave conversion of different scale and different directions Coefficient;
S3, time frequency analysis is carried out to microseism data monitoring data, determines the dominant frequency of microseism useful signal, and according to microseism The dominant frequency of useful signal determines the shearing wave system in useful signal direction corresponding to durection component under different scale respectively Number;
S4, the correlation for calculating useful signal shearing wave coefficient in direction corresponding to direction vector under different scale, root According to the maximum correlation of useful signal shearing wave coefficient corresponding to direction vector under different scale, determine described effective The corresponding shearing wave coefficient of sense and the useful signal direction;
S5, according to the corresponding shearing wave coefficient in the useful signal direction, seek STA/LTA and pick up curve, determine the pickup The curve acquirement maximum value corresponding time is microseism signal first arrival time.
2. the method as described in claim 1, which is characterized in that the difference determination useful signal is in different scale The shearing wave coefficient in direction corresponding to lower durection component specifically:
The size of the corresponding shearing wave coefficient is determined according to the frequency of useful signal height.
3. the method as described in claim 1, which is characterized in that the step S4 specifically:
S41, the correlation that the shearing wave coefficient of corresponding direction described in direction vector under different scale is calculated according to formula (1):
Correlation (k)=U (j, k) * U (j+1, k) (1)
Wherein, U indicates that shearing wave coefficient, j indicate that scale, k indicate direction;
S42, corresponding direction k is maximized as useful signal direction according to Correlation (k), and it is corresponding to obtain direction k Shearing wave coefficient;
S43, the shearing wave coefficient of direction k is added up, the shearing wave coefficient after adding up is the corresponding shearing wave of the useful signal Coefficient.
4. a kind of system for picking up microseism first arrival characterized by comprising
Acquisition module: for acquiring micro-seismic monitoring signal;
Conversion module: for carrying out shearing wave conversion to the micro-seismic monitoring signal, different scale and different directions are obtained Shearing wave transformation coefficient;
Confirmation module: for determining the dominant frequency of microseism useful signal to microseism data monitoring data progress time frequency analysis, and According to the dominant frequency of microseism useful signal, useful signal direction corresponding to durection component under different scale is determined respectively Shearing wave coefficient;
Computing module: for calculating the shearing wave coefficient of useful signal corresponding direction described in direction vector under different scale Correlation, according to the maximum correlation of shearing wave coefficient corresponding to useful signal direction vector under different scale, Determine the corresponding shearing wave coefficient in the useful signal direction and the useful signal direction;
Determination module: for seeking STA/LTA and picking up curve, sentence according to the corresponding shearing wave coefficient in the useful signal direction Determining the pickup curve acquirement maximum value corresponding time is microseism signal first arrival time.
CN201810841073.1A 2018-07-27 2018-07-27 A kind of method and system for picking up microseism first arrival Pending CN109031419A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810841073.1A CN109031419A (en) 2018-07-27 2018-07-27 A kind of method and system for picking up microseism first arrival

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810841073.1A CN109031419A (en) 2018-07-27 2018-07-27 A kind of method and system for picking up microseism first arrival

Publications (1)

Publication Number Publication Date
CN109031419A true CN109031419A (en) 2018-12-18

Family

ID=64645957

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810841073.1A Pending CN109031419A (en) 2018-07-27 2018-07-27 A kind of method and system for picking up microseism first arrival

Country Status (1)

Country Link
CN (1) CN109031419A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110133715A (en) * 2019-05-29 2019-08-16 长江大学 A kind of microseism seismic source location method based on the first arrival time difference and addition of waveforms

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140052377A1 (en) * 2012-08-17 2014-02-20 Schlumberger Technology Corporation System and method for performing reservoir stimulation operations
CN106646598A (en) * 2016-12-27 2017-05-10 吉林大学 FAST-AIC-algorithm micro-seismic signal collecting method
CN107991707A (en) * 2017-12-05 2018-05-04 吉林大学 A kind of borehole microseismic first break picking method based on kurtosis characteristic in shear let domains
CN109031416A (en) * 2018-06-22 2018-12-18 长江大学 The method of microseism P wave first break pickup

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140052377A1 (en) * 2012-08-17 2014-02-20 Schlumberger Technology Corporation System and method for performing reservoir stimulation operations
CN106646598A (en) * 2016-12-27 2017-05-10 吉林大学 FAST-AIC-algorithm micro-seismic signal collecting method
CN107991707A (en) * 2017-12-05 2018-05-04 吉林大学 A kind of borehole microseismic first break picking method based on kurtosis characteristic in shear let domains
CN109031416A (en) * 2018-06-22 2018-12-18 长江大学 The method of microseism P wave first break pickup

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
宋维琪: "《微地震监测新技术与新方法》", 30 June 2014 *
秦晅,等: "基于时窗能量比与互信息量的微地震初至拾取方法", 《物探与化探》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110133715A (en) * 2019-05-29 2019-08-16 长江大学 A kind of microseism seismic source location method based on the first arrival time difference and addition of waveforms

Similar Documents

Publication Publication Date Title
CN111164462B (en) Artificial source surface wave exploration method, surface wave exploration device and terminal equipment
CN112329855B (en) Underdetermined working mode parameter identification method and detection method based on self-adaptive dictionary
CN102288843B (en) Power quality disturbance signal detection method
CN110132403A (en) A kind of vacuum pump vibration signal noise-reduction method based on EEMD and wavelet threshold
CN109443719B (en) Drill bit vibration signal online virtual test method and system thereof
CN106886044B (en) A kind of microseism first break pickup method based on shearing wave and Akaike's Information Criterion
CN110110627B (en) Real-time target detection method for computing resource limitation platform deployment
CN104808243B (en) A kind of pre-stack seismic Bayes inversion method and device
CN103742131A (en) Method for extracting time difference in real time for logging-while-drilling acoustic underground signal acquisition and processing system
CN105212922A (en) The method and system that R wave of electrocardiosignal detects automatically are realized towards FPGA
CN112180429B (en) Unfavorable geological structure detection system and method by utilizing tunnel blasting vibration inversion
CN107966287B (en) Weak fault feature extraction method for self-adaptive electromechanical equipment
CN105259410A (en) Under-sampling waveform frequency estimation method and device under strong noise interference
CN109766798A (en) Gesture data processing method, server and awareness apparatus based on experience small echo
CN109031419A (en) A kind of method and system for picking up microseism first arrival
CN104391325A (en) Discontinuous heterogeneous geologic body detection method and device
CN117214950B (en) Multiple wave suppression method, device, equipment and storage medium
WO2016205679A1 (en) System and method for event detection using streaming signals
CN112711032A (en) Radar target detection method and system based on graph data and GCN
CN103605880B (en) Closely spaced mode damping ratio precisely-diagnosing method
Li et al. Noise reduction method of microseismic signal of water inrush in tunnel based on variational mode method
CN112578439B (en) Seismic inversion method based on space constraint
CN106896404B (en) The recognition methods of thin reservoir and device
CN102394844B (en) FPGA (Field Programmable Gate Array)-based spike potential signal parallel detection device and method
CN105301655A (en) Method and device for eliminating linear noise of common imaging point gather

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
RJ01 Rejection of invention patent application after publication

Application publication date: 20181218

RJ01 Rejection of invention patent application after publication