Automatic identification method of azimuth of planetary seismometer based on passive source

Abstract

The present invention relates to a method for automatically identifying the azimuth of a planetary seismometer based on a passive source. The method includes obtaining optical images of the same area on the planetary surface at different times and solving them to obtain the longitude and latitude of a meteorite impact crater; determining the azimuth from the meteorite impact crater to the planetary seismometer according to the longitude and latitude of the meteorite impact crater and the longitude and latitude of the planetary seismometer; obtaining the three-component seismic signal generated by the meteorite impact crater, extracting the horizontal longitudinal wave signal after bandpass filtering; then, based on the unknown azimuth of the planetary seismometer and the solved azimuth from the meteorite impact crater to the planetary seismometer, constructing an objective function by performing lateral and radial rotation on the horizontal longitudinal wave signal and solving it to obtain the azimuth of the planetary seismometer. The present invention can realize the automatic identification of the azimuth of a planetary seismometer without adding additional designs, and has the advantages of simplicity, efficiency, and accuracy.

Description

Technical Field

[0001] The present invention relates to the field of planetary detection technology, and more particularly to an automatic azimuth identification method of a planetary seismograph based on a passive source. Background Art

[0002] Planetary seismometers are important tools for detecting the internal structure of planets and are of great significance to human understanding of the evolution of planets. At present, the deployment of planetary seismometers is usually carried out unmanned, that is, through landers, due to the high cost, technical difficulty and high safety risks of manual deployment.

[0003] However, there is a key problem with the unmanned deployment method, that is, the azimuth of the planetary seismometer cannot be accurately predicted, which makes it impossible to accurately invert the internal structure of the planet. In response to this problem, Chinese researchers have proposed a series of solutions. Among them, one solution is to install a fixed vibration excitation device on the support leg of the lander to realize the automatic identification of the azimuth of the planetary seismometer (publication number: CN116164699A), but this method has problems such as increased equipment size, power consumption and weight, complex design and high failure rate.

[0004] Another solution is based on the principle of planetary sundials, which realizes automatic identification of azimuth by installing a style and a dial face on the surface of a planetary seismograph and taking photos with an optical camera (Publication No.: CN116481487A). However, this method is greatly affected by weather, its accuracy is affected by its own structure and manufacturing materials, and inaccurate scales on the style or dial face will also affect the accuracy of the results.

[0005] In order to solve the above problems, the present invention proposes a method for automatic identification of the azimuth of a planetary seismograph based on a passive source. Summary of the invention

[0006] The passive source-based planetary seismograph azimuth automatic identification method disclosed in the present invention aims to at least partially solve the above-mentioned technical problems and provide a simple, efficient and accurate solution for my country's planetary exploration.

[0007] In order to achieve the above object, the present invention adopts the following technical solution:

[0008] A method for automatically identifying the azimuth of a planetary seismograph based on a passive source, comprising:

[0009] Obtain optical images of the same area on the planetary surface at different times and calculate the longitude and latitude of the meteorite impact crater;

[0010] Determine the azimuth from the meteorite impact crater to the planetary seismometer based on the longitude and latitude of the meteorite impact crater and the longitude and latitude of the planetary seismometer;

[0011] Obtain the three-component seismic signal generated by the meteorite impact crater, and extract the horizontal longitudinal wave signal after bandpass filtering;

[0012] Based on the unknown azimuth of the planetary seismometer and the calculated azimuth from the meteorite impact crater to the planetary seismometer, the objective function is constructed by performing lateral and radial rotations on the horizontal longitudinal wave signal.

[0013] Solve the objective function to determine the azimuth of the planetary seismometer.

[0014] Preferably, the optical camera on the orbiter is used to continuously take pictures of the planetary surface to obtain optical images of the same area at different times, and the longitude and latitude of the meteorite impact crater are calculated by comparing the optical images.

[0015] Preferably, the azimuth from the meteorite impact crater to the planetary seismometer is determined according to the following formula based on the longitude and latitude of the meteorite impact crater and the longitude and latitude of the planetary seismometer;

[0016]

[0017] Where (latA, lonA) is the longitude and latitude of the meteorite impact crater, and (latB, lonB) is the longitude and latitude of the planetary seismometer.

[0018] Preferably, before bandpass filtering, the instrument response in the three-component seismic signal is deducted according to the following formula:

[0019] x(t) =y(t) / i(t)

[0020] Where x(t) is the original seismic signal after deducting the instrument response; y(t) is the seismic signal recorded by the planetary seismograph; and i(t) is the instrument response.

[0021] Preferably, the bandpass filtering includes a second-order bandpass filtering, and the transfer function of the second-order bandpass filter after Z transformation is:

[0022]

[0023] Where N is the number of zeros and poles, b is N is the coefficient of the forward path, b0, b1, b2 are the serial numbers of the forward path coefficients, a N is the coefficient of the feedback path, and a0, a1, a2 are the serial numbers of the feedback path coefficients.

[0024] Preferably, extracting the horizontal longitudinal wave signal includes:

[0025] The average value of the signal amplitude in the short-time window and the average value of the signal amplitude in the long-time window are obtained at each moment. When the ratio of the average values ​​exceeds the preset threshold, the longitudinal wave signal at the corresponding moment is extracted.

[0026] Preferably, the average value of the signal amplitude in the short-time window is:

[0027]

[0028] The average value of the signal amplitude in the long time window is:

[0029]

[0030] F(i) is the influencing factor at time i, S(i) is the average value of the signal amplitude in the short-time window, L(i) is the average value of the signal amplitude in the long-time window, Q is the total number of sample points of the signal in the short-time window, and M is the total number of sample points of the signal in the long-time window.

[0031] Preferably, the impact factor at time i is calculated as:

[0032]

[0033] Where Y(i) is the amplitude of the signal arriving at time i.

[0034] Preferably, the objective function expression is:

[0035]

[0036] Where, E1(φ) is the transverse longitudinal wave energy, E2(φ) is the radial longitudinal wave energy, A is the ratio of the transverse longitudinal wave energy to the radial longitudinal wave energy, X p is the longitudinal wave signal extracted from the horizontal X-axis direction, Y p is the longitudinal wave signal extracted from the horizontal Y-axis direction, θ c is the azimuth from the meteorite impact crater to the planetary seismograph, φ is the unknown azimuth of the planetary seismograph, S represents the number of P-wave signals, and j represents the sequence number of the current P-wave signal.

[0037] Preferably, when solving the objective function, scanning is performed between 0°-180° at intervals of 1°, and the φ value corresponding to the minimum ratio of the sum of transverse longitudinal wave energy to the sum of radial longitudinal wave energy is the planetary seismograph azimuth corresponding to the current meteorite impact event.

[0038] Preferably, the planetary seismograph azimuths corresponding to multiple meteorite impact events are obtained, and the average value is calculated as the planetary seismograph azimuth, so as to reduce the influence of a single measurement error on the result and improve the reliability of azimuth identification.

[0039] Compared with the prior art, the passive source-based planetary seismograph azimuth automatic identification method disclosed in the present application can more accurately determine the azimuth of a planetary seismograph by positioning the longitude and latitude of a meteorite impact crater in combination with three-component seismic signal processing; and is not affected by the topography and structure of the planetary surface, and is suitable for seismograph azimuth identification on various planetary surfaces.

[0040] Compared with traditional methods, there is no need to add additional structural design, which reduces the equipment size, power consumption and weight, and reduces deployment costs and safety risks. At the same time, it is not affected by weather and can perform azimuth identification under various meteorological conditions, which improves the applicability of the solution.

[0041] This application is significantly innovative, practical and economical, and provides new ideas and means for the development of planetary earthquake detection technology. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying creative work.

[0043] Figure 1 This is a flow chart of the automatic azimuth identification method of a planetary seismograph based on a passive source of the present invention;

[0044] Figure 2 This is a longitudinal wave waveform diagram extracted from the three-component meteorite impact data after bandpass filtering in the present invention. DETAILED DESCRIPTION

[0045] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0046] When a meteorite hits the surface of a planet, a planetary seismograph will collect three-channel seismic wave signals. By analyzing these signals, the azimuth of the planetary seismograph can be accurately identified. Based on this principle, the present invention proposes an automatic azimuth identification method for a planetary seismograph based on a passive source.

[0047] Specifically, after the planetary seismometer is released to the surface of the planet by the lander, the optical camera is used to take pictures of the planetary surface, and then the position of the meteorite impact crater is calculated. Further, based on the seismic signal generated by the meteorite impact collected by the planetary seismometer, the azimuth of the planetary seismometer on the planetary surface is calculated.

[0048] In an practicable embodiment, referring to Figure 1 , the automatic identification method of the azimuth of a planetary seismometer based on passive sources is implemented by the following steps:

[0049] 1) Obtain optical images of the same area on the planetary surface at different times and calculate the longitude and latitude of the meteorite impact crater;

[0050] In this embodiment, the optical camera on the orbiter continuously takes pictures of the planetary surface, and then returns the images to the ground station computer. By regularly comparing the optical images of the same planetary surface area, the longitude and latitude of the meteorite impact crater are calculated.

[0051] The orbiter refers to the Lunar Orbiter, also known as the "Lunar Orbiter", which is NASA's lunar probe and an artificial lunar satellite that studies the lunar environment and surface structure.

[0052] 2) Determine the azimuth from the meteorite impact crater to the planetary seismometer based on the longitude and latitude of the meteorite impact crater and the longitude and latitude of the planetary seismometer;

[0053] Assuming that the longitude and latitude of the meteorite impact crater are (latA, lonA) and the longitude and latitude of the planetary seismometer are (latB, lonB), the azimuth angle from the meteorite impact crater to the planetary seismometer is calculated as follows:

[0054]

[0055] 3) Obtain the three-component seismic signal generated by the meteorite impact crater, and extract the horizontal longitudinal wave signal after bandpass filtering;

[0056] In one embodiment, it is preferred to deduct the instrument response from the three-component ground motion signal;

[0057] In this embodiment, the planetary seismograph transmits the collected three-component seismic signal generated by the meteorite impact to the ground station computer, and then the ground station computer deducts the instrument response from the signal to obtain the three-component original seismic signal, which is expressed as:

[0058] x(t) = y(t) / i(t)

[0059] Where x(t) is the original seismic signal after deducting the instrument response; y(t) is the seismic signal recorded by the planetary seismograph; and i(t) is the instrument response.

[0060] In one embodiment, the original earthquake signal is bandpass filtered, including determining the bandpass filter frequency band according to the epicenter distance, wherein the energy of distant earthquakes is mainly concentrated between 3600s and 2s, the energy of nearby earthquakes is mainly concentrated between 30s and 10Hz, and the energy of local earthquakes is mainly concentrated between 5s and 30Hz.

[0061] In one embodiment, second-order bandpass filtering is performed, and a suitable filtering frequency band is selected according to the epicenter distance, and then a corresponding bandpass filter is designed, wherein the normalized expression of the transfer function of the second-order bandpass filter is as follows:

[0062]

[0063] in, , , , , f s is the sampling frequency, f2 is the upper cutoff frequency of the bandpass filter, and f1 is the lower cutoff frequency of the bandpass filter.

[0064] Further, Substituting into the above formula, we get the system function H(z) after the second-order bandpass filter Z transform, which is expressed as:

[0065]

[0066] Where N is the number of zeros and poles, the numerator and denominator in the formula correspond to the zeros and poles of the filter, respectively. N is the coefficient of the forward path, b0, b1, b2 are the serial numbers of the forward path coefficients, a N is the coefficient of the feedback path, and a0, a1, a2 are the serial numbers of the feedback path coefficients.

[0067] In one embodiment, longitudinal wave extraction is performed on the three-component signals after filtering, and the extraction steps are: obtaining the average value of the signal amplitude in the short-time window at each moment and the average value of the signal amplitude in the long-time window, and when the ratio of the average values ​​exceeds a preset threshold, extracting the longitudinal wave signal at the corresponding moment.

[0068] When a longitudinal wave arrives, there will be a sudden change in the ratio of the signal amplitude in the short-time window at the current moment to the average value of the signal amplitude in the long-time window. When the ratio exceeds the set threshold, it is considered to be a valid signal.

[0069] The expression of the average signal amplitude S(i) in the short-time window at time i is:

[0070]

[0071] Where x(i) represents the signal in the short-time window, and Q represents the total number of sample points of the signal in the short-time window;

[0072] The expression of the average signal amplitude L(i) in the long time window at time i is:

[0073]

[0074] Where h(i) represents the signal in the long time window, and M represents the total number of signal samples in the long time window.

[0075] In a preferred embodiment, when calculating the average value, the influencing factor is taken into account to improve the accuracy of the average value. The influencing factor determined at time i is:

[0076]

[0077] Among them, Y(i) is the amplitude of the signal arriving at time i.

[0078] And substitute the impact factor into the following expression:

[0079]

[0080]

[0081] Furthermore, a recursive processing method is used to determine whether the ratio of S(i) / L(i) exceeds a set threshold. If so, it is considered that a longitudinal wave signal has arrived and the longitudinal wave signal is extracted.

[0082] It should be noted that, since the short window and the long window overlap in time, and the long window needs to contain more historical information than the short window, when calculating the mean of the long window, the influence of the signal Q time units ago needs to be considered to ensure that the mean calculation of the long window is based on the data after a complete short window period, thereby reflecting the average characteristics of the signal over a longer period of time, rather than just based on the most recent data.

[0083] In this application, the extracted longitudinal wave time series in the horizontal X-axis direction, the horizontal Y-axis direction, and the vertical Z-axis direction are respectively recorded as Xp(n), Yp(n), and Zp(n).

[0084] 4) Based on the unknown azimuth of the planetary seismometer and the calculated azimuth from the meteorite impact crater to the planetary seismometer, the objective function is constructed by rotating the horizontal longitudinal wave signal;

[0085] Assuming that the azimuth angle of the planetary seismometer on the planetary surface is φ, and the angle between the line connecting the meteorite impact crater and the planetary seismometer and the north pole of the instrument is θa, then φ and θa satisfy the following formula:

[0086] θa=θc-φ

[0087] Then the horizontal longitudinal wave signal is rotated transversely and radially. The converted transverse longitudinal wave time series is Ts(n), and the radial longitudinal wave time series is Rs(n). The conversion formula is as follows:

[0088] Ts(n)= -Xp(n)*sin(θc-φ)+Yp(n)*cos(θc-φ)

[0089] Rs(n)= Xp(n)*cos(θc-φ)+Yp(n)*sin(θc-φ)

[0090] The ratio A of the transverse longitudinal wave energy and E1(φ) to the radial longitudinal wave energy and E2(φ) is further calculated, and the formula is as follows:

[0091]

[0092]

[0093]

[0094] That is, the objective function expression constructed in this application is:

[0095]

[0096] Where, E1(φ) is the transverse longitudinal wave energy, E2(φ) is the radial longitudinal wave energy, A is the ratio of the transverse longitudinal wave energy to the radial longitudinal wave energy, X p is the longitudinal wave signal extracted from the horizontal X-axis direction, Y p is the longitudinal wave signal extracted from the horizontal Y-axis direction, θ c is the azimuth from the meteorite impact crater to the planetary seismograph, φ is the unknown azimuth of the planetary seismograph, S represents the number of P-wave signals, and j represents the sequence number of the current P-wave signal.

[0097] 5) Solve the objective function to determine the azimuth of the planetary seismometer.

[0098] In this embodiment, when solving the objective function, scanning is performed between 0°-180° at intervals of 1°, with the value corresponding to the minimum ratio of the sum of the transverse longitudinal wave energy and the sum of the radial longitudinal wave energy being the value, and φ is the azimuth of the planetary seismometer corresponding to the current meteorite impact event.

[0099] In a preferred embodiment, multiple meteorite impact events may occur over time. Therefore, the present application calculates each meteorite impact event in sequence according to the above steps, respectively obtains the azimuth corresponding to each event, and then calculates the average azimuth obtained for the meteorite impact events occurring within a fixed time as the final azimuth of the planetary seismometer on the planetary surface.

[0100] An example is as follows:

[0101] Taking a meteorite impact event as an example, the calculation method of the azimuth of the planetary seismometer is described as follows.

[0102] The optical camera on the orbiter continuously takes pictures of the planet's surface, and then returns the images to the ground station computer. By regularly comparing the optical images of the same area on the planet's surface, the longitude and latitude of the meteorite impact crater of the event can be calculated.

[0103] Assuming that the longitude and latitude of the meteorite impact crater returned to the ground station computer are (-3.85°, 152.2°), the location of the planetary seismometer on the planetary surface is known, and assuming that the longitude and latitude are (29.58°, 113.41°), then according to the above formula, it can be calculated that the azimuth from the meteorite impact crater to the planetary seismometer is 314°.

[0104] According to the above steps, the three-component data of the meteorite impact event collected by the planetary seismometer are band-pass filtered at (0.05-1) Hz. Then, the longitudinal waves are extracted from the filtered three-component data respectively. The waveform of the extracted data is shown in the figure below: Figure 2 As shown, the extracted horizontal and vertical longitudinal wave time series are expressed as:

[0105] Xp(n)=(Xp1, Xp2,…, Xp39700),n=1,2,…,39700;

[0106] Yp(n)=(Yp1, Yp2,…, Yp39700),n=1,2,…,39700;

[0107] Zp(n)=(Zp1, Zp2,…, Zp39700),n=1,2,…,39700.

[0108] According to the known azimuth from the meteorite impact crater to the planetary seismometer, combined with the unknown azimuth of the planetary seismometer on the planetary surface, the extracted horizontal longitudinal wave signals are rotated laterally and radially respectively, and then the objective function is constructed to minimize the ratio of the sum of the transverse longitudinal wave energy to the sum of the radial longitudinal wave energy. The obtained azimuth is used as the azimuth of the planetary seismometer on the planetary surface of the meteorite impact event.

[0109] Compared with the active source method (CN116164699A), the present invention can realize the automatic identification of the azimuth of the planetary seismometer without adding additional design, thereby reducing the complexity of the overall design and the failure rate;

[0110] At the same time, the active source (CN116164699A) calculates the azimuth of the planetary seismometer through the azimuth of multiple support legs on the lander. The distance between the support legs on the lander is relatively close, which will lead to large errors in the results. The present method calculates the azimuth of the planetary seismometer based on the azimuth between the meteorite impact crater and the planetary seismometer, so the error of the calculated azimuth will be smaller and the accuracy will be higher.

[0111] Compared with the planetary sundial method (CN116481487A), this method is not affected by external conditions and its own structure and manufacturing materials, so the calculated azimuth is more accurate.

[0112] In this specification, each embodiment is described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the embodiments can be referred to each other. For the device disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple, and the relevant parts can be referred to the method part.

[0113] The above description of the disclosed embodiments enables one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but rather to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for automatic identification of the azimuth of a planetary seismograph based on a passive source, characterized in that: Obtain optical images of the same area on the planetary surface at different times and calculate the longitude and latitude of the meteorite impact crater; Determine the azimuth from the meteorite impact crater to the planetary seismometer based on the longitude and latitude of the meteorite impact crater and the longitude and latitude of the planetary seismometer; Obtain the three-component seismic signal generated by the meteorite impact crater, and extract the horizontal longitudinal wave signal after bandpass filtering; Based on the unknown azimuth of the planetary seismometer and the calculated azimuth from the meteorite impact crater to the planetary seismometer, the objective function is constructed by performing lateral and radial rotations on the horizontal longitudinal wave signal. Solve the objective function to determine the azimuth of the planetary seismometer.

2. The method for automatically identifying the azimuth of a planetary seismograph according to claim 1, characterized in that: The optical camera on the orbiter continuously takes pictures of the planet's surface to obtain optical images of the same area at different times. By comparing the optical images, the longitude and latitude of the meteorite impact crater are calculated.

3. The method for automatically identifying the azimuth of a planetary seismograph according to claim 1, characterized in that: According to the longitude and latitude of the meteorite impact crater and the longitude and latitude of the planetary seismometer, the azimuth from the meteorite impact crater to the planetary seismometer is determined according to the following formula; ； Where (latA, lonA) is the longitude and latitude of the meteorite impact crater, and (latB, lonB) is the longitude and latitude of the planetary seismometer.

4. The method for automatically identifying the azimuth of a planetary seismograph according to claim 1, characterized in that: Before bandpass filtering, the instrument response in the three-component ground motion signal is deducted according to the following formula: x(t) =y(t) / i(t) Where x(t) is the original seismic signal after deducting the instrument response; y(t) is the seismic signal recorded by the planetary seismograph; and i(t) is the instrument response.

5. The method for automatically identifying the azimuth of a planetary seismograph according to claim 1, characterized in that: Bandpass filtering includes second-order bandpass filtering. The transfer function of the second-order bandpass filter after Z transformation is: ； Where N is the number of zeros and poles, b is N is the coefficient of the forward path, b0, b1, b2 are the serial numbers of the forward path coefficients, a N is the coefficient of the feedback path, and a0, a1, a2 are the serial numbers of the feedback path coefficients.

6. The method for automatically identifying the azimuth of a planetary seismograph according to claim 1, characterized in that: Extraction of horizontal longitudinal wave signals, including: The average value of the signal amplitude in the short-time window and the average value of the signal amplitude in the long-time window are obtained at each moment. When the ratio of the average values ​​exceeds the preset threshold, the longitudinal wave signal at the corresponding moment is extracted.

7. The method for automatically identifying the azimuth of a planetary seismograph according to claim 6, characterized in that: The average value of the signal amplitude in the short-time window is: ； The average value of the signal amplitude in the long time window is: ； F(i) is the influencing factor at time i, S(i) is the average value of the signal amplitude in the short-time window, L(i) is the average value of the signal amplitude in the long-time window, Q is the total number of sample points of the signal in the short-time window, and M is the total number of sample points of the signal in the long-time window.

8. The method for automatically identifying the azimuth of a planetary seismograph according to claim 7, characterized in that: The impact factor at time i is calculated as: ； Where Y(i) is the amplitude of the signal arriving at time i.

9. The method for automatically identifying the azimuth of a planetary seismograph according to claim 1, characterized in that: The objective function expression is: ； Where, E1(φ) is the transverse longitudinal wave energy, E2(φ) is the radial longitudinal wave energy, A is the ratio of the transverse longitudinal wave energy to the radial longitudinal wave energy, X p is the longitudinal wave signal extracted from the horizontal X-axis direction, Y p is the longitudinal wave signal extracted from the horizontal Y-axis direction, θ c is the azimuth from the meteorite impact crater to the planetary seismograph, φ is the unknown azimuth of the planetary seismograph, S represents the number of P-wave signals, and j represents the sequence number of the current P-wave signal.

10. The method for automatically identifying the azimuth of a planetary seismograph according to claim 1, characterized in that: Solve the planetary seismograph azimuth corresponding to multiple meteorite impact events, and calculate the average value as the planetary seismograph azimuth.

Images