• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    A New Method of Frequency Fluctuation Estimation and IPS Processing Results Based on the Downlink Signal of Tianwen-1

    2023-09-03 01:37:12YuChenLiuDeQingKongDongHaoLiuXinYingZhuYanSuHongBoZhangandZanWang

    Yu-Chen Liu ,De-Qing Kong,? ,Dong-Hao Liu ,Xin-Ying Zhu ,Yan Su ,Hong-Bo Zhang ,and Zan Wang

    1 Key Laboratory of Lunar and Deep Space Exploration,National Astronomical Observatories,Chinese Academy of Sciences,Beijing 100101,China;kdq@bao.ac.cn

    2 University of Chinese Academy of Sciences,Beijing 100049,China

    Abstract The radio-occultation observations taken by Tianwen-1 are aiming to study the properties of solar wind.A new method of frequency fluctuation (FF) estimation is presented for processing the down-link signals of Tianwen-1 during the occultation period to study the properties of the coronal plasma at the heliocentric distances of 4.48–19 R⊙.Because of low S/N as well as the phase fluctuation phenomena caused by solar activity,a Kalman based on polynomial prediction methods is proposed to avoid the phase locked loop loss lock.A new detrend method based on multi-level iteration correction is proposed to estimate Doppler shift to get more accurate power density spectra of FF in the low frequency region.The data analyze procedure is used to get the properties of the solar corona during the occultation.The method was finally verified at the point when the solar offset is 5.7 R⊙,frequency tracking was successfully performed on data with a carrier-to-noise ratio of about 28 dBHz.The density spectra obtained by the improved method are basically the same when the frequency is greater than 2 mHz,the uncertainty in the result of the rms of the FF obtained by removing the trend term with different order polynomials is less than 3.3%.The data without eliminating interference show a large error for different detrending orders,which justifies the need for an improved approach.Finally,the frequency fluctuation results combined with the information on intensity fluctuation obtained by the new method are compared with the results of the integrated Space Weather Analysis system and theoretical formula,which verifies that the processing results in this paper have a certain degree of credibility.

    Key words: methods: data analysis– Sun: corona– Sun: heliosphere– (Sun:) solar wind– Physical Data and Processes

    1.Introduction

    Interplanetary scintillation (IPS) refers to electromagnetic wave signals that are affected by the activity of the outer coronal plasma layer of the Sun that they pass through,thus exhibiting rapid changes in parameters such as angle of arrival,amplitude,and phase (Xu &Song 2019).

    Early IPS research mainly used natural radio sources(Hewish 1955;Hewish et al.1964)to invert the characteristics of the Sun.However,natural radio sources can only provide information on amplitude scintillation,the use of electromagnetic waves emitted by deep space detectors to conduct IPS experiments can further explore the structure of the corona by using the frequency fluctuation(FF),group delay,and Faraday rotation of the signal.In the late 1960s,Goldstein studied the relationship between spectral broadening and overall electron density in the coronal region using data from Mariner-4,the first spacecraft to successfully fly by Mars,and later Mariner-6(Goldstein 1969).In terms of Doppler measurements,Berman and Wackley used Pioneer 10,11 and Helios 1 data in 1975 to show that Doppler scintillation is a good way to study solar wind disturbances (Berman &Wackley 1976),Berman further showed in 1981 that the mean fluctuation of the Doppler data is proportional to the total electron density (Berman et al.1982).Research by Woo (1988) and Woo &Armstrong (1992)showed that Doppler scintillation is not only affected by heliocentric distance,but the dynamic solar wind also causes instantaneous scintillation in Doppler data.By comparing experiments with plasma measurements,Woo et al.(1985)found that interplanetary shocks are generally responsible for this instantaneous scintillation during periods of high solar activity,while Doppler scintillation is sometimes caused by corotating particle streams during periods of low solar activity.

    Scintillation research during the radio occultation period has always been a common concern of the communication field and the astronomical community.Through the study of interplanetary scintillation,more information on the physical parameters of the solar wind can be obtained.On the other hand,it is of great scientific significance to study and understand the Sunterrestrial space environment to provide a guarantee for the stable operation of deep space probes.

    At the beginning of this century,a series of IPS studies were also carried out in China (Zhu et al.2012;Liu et al.2021),which were all based on ground observation settings to observe natural radio sources.With the successful launch of the Tianwen-1 deep space probe and the advancement of more deep space missions,it is urgent to establish a systematic data processing method for IPS research using deep space probes.A few months ago,scientists used VLBI methods,combined with Mars Express and Tianwen-1 data,to probe the solar corona(Ma et al.2022).Similar efforts have been carried out by researchers in the past decade,and their methods have been applied to such as Messenger (Wexler et al.2019a),Mars Express (Bocanegra-Bahamón et al.2018;Ma et al.2021),Akatsuki (Wexler et al.2020a;Chiba et al.2022) and other missions.The inversion of solar parameters has achieved good results.On the basis of them,this paper proposes an improved data processing method for calculating more accurate phase scintillation information based on the characteristics of the Tianwen-1 detector signal.

    In Section 2,two schemes based on median filtering to improve the accuracy of frequency estimation in deep space exploration in highly dynamic environments are proposed,which have different purposes but similar principles.In Section 2.2,when the Phase-Locked-Loop (PLL) may lose lock under high dynamic noise,a judgment of interference of PLL based on Multi-Level Iteration Correction (MLIC) is proposed,so that the PLL can be assisted by Kalman filter to track Doppler frequency when the loop is at the edge of out-oflock,thereby eliminating the tracking error of phase-locked loops due to high dynamics and eliminates the wild value.In Section 2.3,based on the previous outlier judgment scheme,a new method is proposed to fit the long-period slow trend of frequency that results in more accurate frequency residuals for analysis of solar activity.In Section 3,the above two schemes are verified with the real data from Tianwen-1,and the accuracy of the results is compared.

    2.Method of Frequency Fluctuation Estimation

    2.1.The Characteristics of the Downlink Data

    The downlink data received from Tianwen-1 are transmitted by medium gain antenna in theX-band,during superior conjunction of the Mars over solar offset (SO) range from 4.5 to 19R⊙(solar radii),which defined as the distance from the proximate point to the heliocenter.In order to obtain the phase scintillation (Doppler scintillation) information of the signal,the one-way Doppler measurement mode is adopted(Zheng 2013),which means after the satellite transmits the beacon signal of the nominal frequency,the ground radio telescope completes the signal reception and frequency measurement.The single-station single-frequency measurement mode was adopted in this experiment.However,the weak signal transmitted from a medium gain antenna and the strong interference from the inner solar wind bring difficulties in frequency measurement.Figure 1 shows the spectrogram of the Tianwen-1 downlink signal received by the Wuqing 70 m radio telescope (WRT70) on 2021 October 12.

    Figure 1.Power density spectra for Tianwen-1 at WRT70 at 08:02:09 when R⊙=5.7.The sampling rate is 6.25 MHz with 8 bit digitized,integration time is 32 s,spectral resolution is 5.96 Hz.It can be seen that due to the Doppler shift,the frequency point of the signal is not at zero.

    Costas loop is used to obtain the instantaneous frequency(Zheng 2013;Buccino et al.2018),and after deducting the long-time varying trend caused by the relative motion between the satellite and Earth,we obtain the residual frequency which provides the information on solar activity.In view of the characteristics of low signal-to-noise ratio (S/N) and large scintillation of downlink data,this section will propose some new methods suitable for processing Frequency Fluctuation data for a single station.

    We observe density fluctuations using electromagnetic signals traversing the Sun’s corona,which are mainly affected by the proximate point on the line of sight (LOS).The signal frequency observed at the receiver system,fobs,has the following relationship with the original transmitted frequencyf0(P?tzold et al.2016;Wexler et al.2019b)

    where λ=c/f0is the transmitted carrier frequency,cis the speed of light,neis the electron number density function with respect to the LOS integration pathsand timet,andreis the classical electron radius which equals to 2.82×10?15m.

    The second term in Equation (1),?f0Vrel/c,is the Doppler frequency shift,which comes from the relative velocity,Vrel,between the detector and the Earth receiving station,which usually manifests as a long-period slow motion of the signal(Ma et al.2021).The method of eliminating the long period will affect subsequent frequency fluctuations spectral power estimation,especially causing fitting errors in the lowfrequency part (Wexler et al.2020a).In the past,there are mainly two processing methods to remove the gradual change caused by relative motion,which are direct polynomial fitting and the ephemeris method.The literature Wexler et al.(2020a)proposes using a third-order polynomial to fit the data every 4000 s to eliminate long periods.Ma et al.(2021)use a higherorder polynomial for fitting,and compare the method of using the satellite ephemeris to eliminate the period.It is found that the errors caused by the two methods are concentrated below 3 mHz.However,for the data of this study,neither simple polynomial fitting nor ephemeris can provide a solution precise enough.

    Because the local Doppler data sometimes has a large scintillation,it will lead to the polynomial long-period fitting affected by the local jump,and the error of eliminating the period will be brought about.For Mars exploration,the computational complexity of ephemeris is very high,and the accuracy of calculating relative motion through ephemeris is affected by the long distance,long delay,and complex relative motion in deep space (probe autonomous orbit adjustment,Mars revolution,Earth Rotation),many effects lead to high computational complexity,and this also means that the Doppler scintillation data cannot be processed immediately and accurately when the downlink data obtained.After removing the slow motion of the signal,we can obtain the instantaneous signal changes,so as to further study the influence of the solar corona on the signal.

    Fluctuation varianceis defined to infer the number density fluctuation information based on the numerical integration of FF spectral density over frequency limits [νlo,νup].(Wexler et al.2019b)

    where { }F ...represents the Fourier Transform.The lower frequency limit of the FF spectral integration νlois affected by the length of observation and the detrending method.In previous work (Ma et al.2021),it was found that the effect of the detrending term on the results above 3 mHz was negligible.The upper limit νuprequired to below the Nyquist frequency and noise floor (Wexler et al.2020a).

    The digital PLL (Molera Calvés 2012;Yuldashev 2013;Buccino et al.2018) is used to estimate the instantaneous frequency and phase.When phase-locking is successful,the Doppler frequency estimated by the PLL has a strong longperiod trend with time.For the phase disturbance record of the Cassini satellite in 2000 (Morabito et al.2003),the frequency shift is on a second scale and the long-time regularity presented is not caused by the disturbance of the solar wind.This is mainly caused by the influence of the relative motion of the earth.In order to obtain the information on the solar activity,that is,the third information of Equation (1),it is necessary to remove the slow change caused by the relative motion of the aircraft,which is a quasi-sine wave shape on the long-period scale.After the long-period slow change is eliminated,the residual frequency reflecting the influence of the solar wind can be obtained.

    To this end,we proposed a long-period fitting method based on multi-level iteration correction and denoising to eliminate the long-period error of the signal.

    2.2.PLL Tracking Based on Adaptive Filtering

    In the field of deep space exploration,especially during the occultation period,solar activity has a great influence on Doppler measurements,so that the effect of solar activity can cause instability of the PLL.The PLL-based closed-loop parameter estimator is a classical frequency phase estimation method(Gardner 2005;Molera Calvés 2012;Buccino et al.2018;Deng et al.2021).Due to the harsh conditions of the down-link narrow-band signal of deep space detection,it is easy to have a low S/N of the received signal and a dynamic frequency variation beyond the allowable range of the PLL,resulting in serious problems such as the inability to track,loss of lock,and large parameter estimation error (Paik &Asmar 2011).

    In order to solve the problem of frequency jump and loss of lock of PLL under solar activity,the following frequency estimation and tracking structure is proposed in Figure 2.

    Figure 2.Flow chart of the frequency estimation algorithm.

    The data processing flow is divided into two working modes,Frequency Tracking Mode and Frequency Correction Mode,switching depends on whether the PLL is at the edge of out-oflock.Under normal conditions,enter the Frequency Tracking Mode.In this mode,the initial frequency offset due to the Doppler effect is first estimated by the Fast Fourier Transform(FFT),and then the estimated value of the initial frequency offset is used as a priori condition to assist the locking of the PLL and get the estimated Doppler frequency.At the same time,the estimated frequency offset obtained after the PLL is used as feedback to assist the subsequent frequency tracking.Through the interference detection of the output results from the PLL which will discuss in Section 2.2.2,if the solar wind causes relatively large frequency scintillation to break the loop of PLL,the Frequency Correction Mode is entered,and the Doppler frequency is predicted and the instantaneous frequency is corrected by the Kalman filter.

    2.2.1.Frequency Tracking Mode based on PLL

    The down-link narrow-band signal transmitted by the deep space probe received by the ground receiving station,the signal after down-converted with sampling ratefsares(n),which can be expressed as

    wheref(n)=δf(n)+fc+fdand φ(n) are the instantaneous received frequency and phase changing in time,andfcis the frequency at which the spacecraft transmits signals,fdis the Doppler frequency,a(n)is the modulated information of signal.

    The received signal multiplied by the cosine signal output by the numerically controlled oscillator (NCO) can be expressed as (Yuldashev 2013)

    wherefest(n) and(n) are the estimation of instantaneous frequency and phase.

    After low-pass filtering,it can be obtained,where Δf=f(n)?fest,and Δφ(n)=φ(n) ?(n)

    The rough estimation of initial Doppler frequencyfinican be calculated by FFT(Molera Calvés 2012),with the resolution as ΔfFFT=fs/NFFT,whereNFFTis the number of FFT counts points.The accuracy of the initial frequency estimation affects the capture time required before the phase-locked loop enters the tracking state.

    2.2.2.Interference Detection

    The movement of the deep space detector can be described as follows without considering the extreme situations such as EDL (Entry,Descent,and Landing) (Paik &Asmar 2011;Niu et al.2014),the Doppler frequencyfd(t) can be represented as(Qi et al.2020)

    Since the S/N of the return signal during superior conjunction of the Mars is extremely low,in order to improve the S/N of the signal and simplify the subsequent calculations,we integrate the Doppler data.Without losing its generality,everyjdata have been integrated at the current timei.

    Then perform median filtering onfsum(i)

    wherekrepresents the current moment.In order to determine whether there are jittering disturbing points in the Doppler,we make a difference between the results after median filtering to get the difference of amplitude

    However,the above operation is based on the assumption that the polynomial fitting results(t)of the time intervalttot?Tis accurate at the momentt.In fact,Doppler estimationfdduring the upper conjunction period is greatly affected by the instantaneous activity of the solar wind,and Figure 3 shows the resulting jumps in frequency estimation that inevitably affect the accuracy of fitting trends directly using polynomials.In Figure 4(a),the error in estimating the trend term caused by direct polynomial fitting (the green line in Figure 4(a)) can be seen visually.At the same time,this is bound to affect the determination of interference points discussed above.

    Figure 3.The comparison of Δ (k)and the abnormal points in fd.The shaded areas in red represent the intervals detected as points of interference.

    Figure 4.Frequency residuals extracting process(a)Comparison of traditional detrending method and the Kalman filter to predict Doppler shift.(b)Prediction performance of the Medfilt-based Kalman filter.(c) Comparison after different detrend methods.The upper one corresponds to the time domain diagram obtained after the direct polynomial detrend method(set )with p=5.(d)This is the plot of FF spectraldensity.The green line and yellow line in the figure are the spectral density deduced from)of order 5.The lower figure is the time domain diagram using the new approach(set and ,while red and blue line are calculated byand .

    In order to make the trend estimation results of polynomial fittingfit the real long-term change trend,an improvement scheme of multi-level iteration correction is proposed in the following.For illustrative purposes,let us assume that the discrete numerical coordinates corresponding to the time fromt?Ttotarek=1 tok=M.In order to eliminate the jump interference point,a suitable threshold is required.We adopted a multi-level iteration method to set the threshold.To eliminate interference points,the process of themthtime is as follows:

    (1)Find average.To eliminate larger amplitude values,the average is calculated as

    Fqis the set of frequency pointkq,iof the interference signal detected by the q-level interference detection module.And

    especially,defineF0=φ whenq=0.

    (2)Calculate variance.In the second step,find the variance of the amplitude difference

    (3)Determine interference via threshold.This article decides to set the threshold to the following form

    When Δfmed(k)≥Qm,the correspondingkcan be calculated.Assume that there areL?s value that meets the above conditions,to represent as.

    (4)Delete interference and their neighbor.For our purpose,the interference points,which are caused by the effect of loss of lock of PLL,do not occur in isolation.Therefore,a neighborhood of the interference points should be chose to delete.

    2.2.3.Frequency Correction Mode Based on Kalman Filter

    We introduce Kalman filtering into the PLL to denote the state equation and observation equation as follows

    where the observation value z(k) is taken to be(k)that calculated in Section 2.2.2,H(k),V(k) are the observation vector,measurement matrix and observation noise at timek.Transmission matrix,measurement matrix and state vector are respectively

    The new Doppler frequency estimation can be given by the following algorithm (Niu et al.2014;Qi et al.2020)

    with the initial condition

    Replace the frequency value of thek-th cycle in Equation(8)with the Kalman filter value as follows.

    Then we can get the result of frequency correctionfcorrec(tk),which is an estimation offest(tk)based on the result of the previous cycle of PLL,or frequency correction after suppressing Doppler frequency jump caused by instability of PLL.The total result of frequency tracking that achieves the purpose of suppressing out-of-lock effects of PLL on Doppler frequency estimation can be expressed as

    For the convenience of the following article,the traditional tracking results without out-of-lock effects suppression can be listed,where we can define the range of the observation time setTobs.

    I think he is not in his right senses! said the Princess, and walked on, butwhen she had gone a little way, she stopped again. One must encourage art, said she, I am the Emperor s daughter. Tell him he shall, as on yesterday,have ten kisses from me, and may take the rest from the ladies of the court.

    2.3.Long-period Detrending Based on Multi-level Iteration Correction

    In the previous section,this article introduced a signal processing method to suppress frequency jump due to loss of stabilization of PLL at a low S/N and obtained a set of frequency estimation resultsG with wild value eliminated.

    After removing the interference caused by the phase-locked loop,the next step is to eliminate the long-period trend of the time series.As it had been discussed in Equation (2),the effective information of solar activity is not related to the longperiod variation of the Doppler frequency due to the movement of the satellite.The task in extracting solar activity parameters is to efficiently obtain the residuals after removing the trend term.The traditional method of removing trend terms is to remove the direct polynomial fitting introduced in previous Section 2.2.2,and there is also a method of removing the residuals using ephemeris (Ma et al.2021).In general,in the case of near-Earth satellites,ephemeris fitting is a method that can be considered,but for remote-distance,multi-planetary motion coupled deep space exploration,this method not only cannot process the latest data in real-time but also brings excessive calculations.

    Since the S/N of the return signal during superior conjunction of the Mars is extremely low,in order to improve the S/N of the signal and simplify the subsequent calculations,we integrate the Doppler data as before,where integrate everyknumdata

    Then a median filter was performed on the integrated Doppler datafsum(i) as in Equation (12).

    Because the operation at this time is not limited to the information at the current moment as in Section 2,in Equation (16) becomes

    whereNis the largest time point in the observation setTobs.

    Referring to the operation of Equation (17) filtering the interference point,we get Δfmed.It is important to note that although the formula is identical,the previous Δfmedis based on the result of the current momentkderived only from the information at the previous moment,and the Δfmedin this section is the result based on the global information.So now the judgment will be more accurate.

    The larger value of Δfmedcorresponds to the frequency value of jumping after being affected by coupling of solar wind and loss of lock of PLL.As introduced in Section 2.2.2,these jump values will affect the accuracy of the subsequent longterm fitting.

    Similarly,the frequency residuals deduced by traditional method can be expressed as

    3.Result and Test

    An extensive series of interplanetary scintillation (IPS)observations was made by Tianwen-1 down-link signals during the superior solar conjunction of Mars over solar offset range 4.48–19R⊙in 2021.

    To test the effectiveness of the proposed method,we choose the downlink data on 2021 October 12 from 07:59:00 ~10:00:00 (UTC) as the experimental data when the SO value is 5.7R⊙,which is in the middle corona region(from 1.5 to 6R⊙(Wexler et al.2022)).The signal is observed atX-band and received by the Wuqing 70 m radio telescope(WRT70) that sampled at the rate of 6.25MHz with 8 bit digitized,and the phase locking bandwidth is set to be 20 Hz.Figure 1 is the power density spectra for Tianwen-1 at 08:02 on October 12.

    3.1.Frequency Fluctuation Data

    In order to visually illustrate the previous processing method,we selected a data processing result of about 33 s which starts at 8:07:38 to draw Figures 4(a)–(c).Figure 4(a)is the comparative verification of the long period elimination method based on multi-level iteration correction proposed in Sections 2.2 and 2.3.The blue line in Figure 4(a)is the Doppler frequency estimationfdmade only by PLL without correction.The green line in the figure is the result of the traditional direct use of polynomial fitting,where the order of fitting is taken to bep=5.It can be clearly seen that the traditional direct fitting will cause the distortion of the fitting trend due to the jitter of the instantaneous frequency.The red line in the figure is the result of(k∣k)kin Equation (34).Even in the case of severe jitter,the Doppler trend can be well estimated.

    Figure 4(b) can well show the comparison between the frequency tracking result of the classic PLL method(setO)and the result processed by the algorithm in this paper(setG).The red line in this figure is the frequency tracking result of the original PLL.It can be seen that due to the close distance from the Sun (SO=5.7R⊙),the frequency estimation is affected by the activity of the solar corona,and the jitter occur at the same time.The blue line is the advanced Tracking method discussed in Section 2.2.It can be clearly seen that the improved processing method greatly improves the stability of frequency tracking.

    Figure 4(c) is a time domain plot of frequency residuals obtained using two different data processing methods,and its time coordinates are the same as the previous two figures.The figure above is the frequency residuals obtained by subtracting the directly fitting trend term represented by the green line,in Figure 4(a),from the instantaneous frequency tracking result that has not been corrected by the Kalman filter,that is the setwithp=5.It can be observed that the slow trend term is not completely removed when the frequency jump is large.And the frequency jumps generated in the PLL are not entirely affected by the signal as it passes through the solar corona.The frequency jumping in the upper Figure 4(c) for up to several seconds is a loss of lock due to the properties of the loop when tracking a highly dynamic target (Gardner 2005).The lower time domain plot is obtained by using the Kalman filter to compensate for the frequency jump,that is the setwithp=5.At this time,the method of detrending is also more refined.And before the trend term is fitted,there is a disturbance point of frequency jump after the median filtering,which is the method represented by the red line in Figure 4(a).

    Figure 4(d) is a comparison of frequency fluctuation power spectra plotted with data of about 3000 s.The process of plotting the power spectrum can be found in the work of previous (Imamura et al.2005;Wexler 2020b).In this figure,we not only add the comparison of two different detrend schemes (the green and red line) add,but also make a comparison of the fitting results of the scheme under different orders(the red and blue line).It can be seen that compared with the traditional long-period fitting scheme,the result of the improved method is not sensitive to the fitting order.In order to accurately estimate the electron spectral density (Ma et al.2021),the trend fitting scheme proposed is very necessary.For the improved approach,the spectra are almost the same when ν ≥2 mHz.In terms of accuracy,it is 1 mHz higher than the work of previous (Ma et al.2021).

    The large perturbation of the Doppler data in Figures 4(c)and(d)is due to the coupling of solar atmospheric activity with the PLL out-of-lock effect.

    These disturbances affect the elimination of long-term trends in the data,resulting in subsequent calculations of the error of the electron density spectrum in the low-frequency range.

    Once the long-time trend term is removed and the frequency residuals are obtained,we can calculate the FF spectral density.The methods taken to calculate the frequency fluctuation spectra are described by Imamura et al.(2005).Chiba et al.(2022)study the frequency fluctuation spectra by taking 2000 s sub-intervals per day with a resolution of 1.31 s using Akatasuki data detecting the heliocentric distance of 1.5–8.9R⊙.Bruno et al.(2009)point out that the time resolution of frequency fluctuation to study the FF spectra is just required to be smaller than the local ion gyro-period with enough long time duration compared with the period (Bruno &Carbone 2013).Taking into account the effective observation time and the amount of calculation afterward,3000 s time segment and 0.1 s time resolution is used.(See Figure 4(d)).

    where λ is the radio transmitter wavelength in meter which isX-band here,ris the SO value in solar radii.The fluctuation variance σFMthat based on different methods,which shows in Figure 4(d),are listed in Table 1 (Wexler et al.2020a).

    Table 1Information Deduced from the FF Spectral Density on 2021 October 12 from 08:00 to 08:50 (UTC)

    It can be seen that the data processed with the new detrending method is not sensitive to the order of the fit,while the traditional method is the opposite.The new method uses two different fitted orders to obtain σFMwith an uncertainty of less than 3.3% of the electron number density compared to 8.8% in a traditional way.

    The electron number densityNecan be calculated from the empirical model(Wexler et al.2019a)under the assumption of average case of slow wind circumstance

    Then,we can calculate the theoretical plasma density on 2021 October 12,which is 2.5716×1010m?3.During periods of calmer solar activity,an interpretation of FF in the inner part of the middle corona for givenneis proposed by Wexler et al.(2022)

    A term incorporates information of solar wind outflow speed is added for larger SO

    Figure 5 shows these two curves and the rms frequency fluctuation results σFMobtained using different methods.

    Figure 5.Frequency measure fluctuation results.Areas shaded in pink represent the extent of the middle coronal.Equation(44)and a correction from Equation(45)for the information on solar wind speed surrounding the shaded area of light blue.The two dots in red are the result of detrending by different orders using the improved method,and the maximum uncertainty is chosen.Similarly,the blue dots represent the results obtained by data with interference.

    The figure once again verifies that for the improved method,the detrending of different orders has little effect on the results,which almost coincide.While for the traditional method,the difference can be seen.The results obtained by improving the method are consistent with the theoretical formula.

    The results with interference are much larger than the theoretical value.But Wexler et al.(2020a)also suggest that if there were the presence of coronal mass ejection (CME),the rms frequency fluctuation measure σFMwould be much larger than the theoretical value predicted by Equation (45).We will verify in the next section that the solar wind outflow speed at that time meets the conditions for low speed,which excludes the possibility of CME.

    3.2.Comparative Analysis of Intensity Scintillation Data

    In order to verify the premise of the formula in the previous section,the case of the quiet solar wind,we combine the information reveals by intensity scintillation to get the inversion of the solar wind properties.We apply a similar data reduction and processing procedure for single station single frequency(SSSF)given by Liu et al.(2021)to get the intensity power spectral (see Figure 6(a)),the integration time is 1ms with a length of ≈100 s.

    Figure 6.The Interplanetary plasma properties for the time of our observations (2021 October 12 at 06:57 UTC).(a) Intensity scintillation spectrum and its corresponding model inversion.(b)Frequency scintillation spectrum obtained by the new detrend approach in this paper.(c)(Right)Electron density content times the square of the radial distance from the Sun[N(cm?3)].(Left)Solar wind radial speed on the elliptic plane in km s?1.It can be seen that the solar activity at that time was relatively calm.Simulations are available through the Integrated Space Weather Analysis System (iSWA,https://ccmc.gsfc.nasa.gov/iswa/).

    Under the weak-scattering assumption,the intensity scintillation spectrum can be expressed as follows (Scott et al.1983;Liu et al.2021;Chiba et al.2022)

    wherereis the classical electron radius and λ is the observing wavelength which is determined as we use the singleX-band for our experiment.Vis the component of solar wind velocity perpendicular to the LOS.

    To fit the theoretical intensity scintillation spectra,L-BFGS(Limited-memory Broyden–Fletcher–Goldfarb–Shanno) method is applied to obtain the global optimal (Chiba et al.2022).The accuracy of the measurement of solar wind speed and other three parameters will suffer selection effects during the fitting process of SSSF (Zhang 2007).Since the parameter AR andsimainly affect the slope between Fresnel knee and noise floor(Scott et al.1983),on the basis of previous research (Liu et al.2021;Chiba et al.2022),we added new optimization goals to the model fitting of SSSF

    whereWis the weight andKis the inclined rate between Fresnel knee and the noise floor.

    We compare the solar activity parameters based on the inversion of intensity (Figure 6(a)) and frequency scintillation(Figure 6(b)) on 2021 October 12 with the interplanetary observations of that day (Figure 6(c)) provide by NASA’s the Integrated Space Weather Analysis System(iSWA).According to iSWA’s,the speed of the solar wind at about 06:57 on 2021 October 12 is low.However,iSWA’s observations can only indicate that solar activity was relatively quiet at the interplanetary scale (1 au ≈215R⊙),which could not provide the speed of the solar wind in the middle corona region.The activity of the solar middle corona was also relatively calm at that time is the only information that can be inferred.The inversion of solar wind velocity from intensity scintillation is the radial component of the outflow velocity,approximately at 80 km s?1,which consists of iSWA’s model and empirical values at such SO(Wexler et al.2020a),and better reflects the calm state of the solar activity.Therefore,the conditions at that time make the empirical Equations (43),(44) and (45) valid.This result illustrates the effectiveness of the treatment method proposed in this paper.

    4.Conclusion and Prospect

    In this paper,a systematic Doppler frequency measurement method for Mars exploration missions and its application to the estimation of solar activity region parameters are proposed.The main conclusions are:

    1.A scheme of frequency tracking and correction is proposed in the presence of high dynamic Doppler noise(such as Sun noise)in the deep space communication scene.In the actual test when SO is less than 6R⊙,X-band,rapid changes in frequency were successfully tracked by PLL.

    2.In order to measure the electron density fluctuation spectrum in the LOS direction by using the Doppler scintillation data of spacecraft radio,a more effective long-period elimination scheme is proposed,which is verified to correct the error of the frequency scintillation spectrum in the low frequency band,the uncertainty of estimating σFMimproves from 8.8% to 3.3%.This increases the accuracy trusted in the lower frequency limit of power spectrum integration νloto 2 mHz.

    3.Using the inversion method of intensity scintillation combined with NASA’s observations,it was determined that the Sun was in a relatively calm state around 8:00 a.m.on 2021 October 12.However,iSWAs data can only provide information at the interplanetary scale (1 au≈215R⊙).In order to obtain the parameters of the solar corona,the paper combines with intensity scintillation data to reflect the results of the quiet solar wind,verifying the relatively calm state of the solar corona during observation.Thus,it can be verified that the results obtained by the improved method are consistent with the theoretical values.

    This paper verifies the advantages of the improved method for calculating FF,and verifies that the results are consistent with the theoretical values,and that in the middle coronal,the frequency is affected byandr?3.2.The second of these represents the effects of the low-speed solar wind.The use of intensity scintillation data can just reflect the information on solar wind speed,thereby improving the credibility and accuracy of the fitting.In the next step,we will use the improvements proposed in this paper to make an overall statistical analysis of the data during the Tianwen-1 occultation period.

    Acknowledgments

    This work was funded by the Astronomical Joint Fund of the National Natural Science Foundation of China and Chinese Academy of Sciences (Grant Nos.U1831114,11941002,and 12073048).

    日韩亚洲欧美综合| 精品久久久久久电影网 | 中国国产av一级| 成人特级av手机在线观看| 亚洲伊人久久精品综合 | 最近视频中文字幕2019在线8| 欧美一区二区国产精品久久精品| 爱豆传媒免费全集在线观看| 免费播放大片免费观看视频在线观看 | 亚洲av成人精品一二三区| 永久免费av网站大全| 日本-黄色视频高清免费观看| 国产老妇女一区| 亚洲经典国产精华液单| 亚洲色图av天堂| 国产精品1区2区在线观看.| 欧美激情久久久久久爽电影| 精品久久久久久久人妻蜜臀av| 午夜久久久久精精品| 亚洲天堂国产精品一区在线| 内地一区二区视频在线| 蜜桃久久精品国产亚洲av| 久久精品国产鲁丝片午夜精品| 免费黄网站久久成人精品| a级毛色黄片| 色视频www国产| 国产中年淑女户外野战色| 99热全是精品| 99久国产av精品| 亚洲国产精品久久男人天堂| 麻豆精品久久久久久蜜桃| 禁无遮挡网站| 熟女人妻精品中文字幕| 99热6这里只有精品| 非洲黑人性xxxx精品又粗又长| 秋霞在线观看毛片| 51国产日韩欧美| 精品人妻熟女av久视频| 精品人妻熟女av久视频| 色视频www国产| 亚洲国产成人一精品久久久| 特大巨黑吊av在线直播| 亚洲国产精品国产精品| 久久久a久久爽久久v久久| 国产欧美日韩精品一区二区| 麻豆一二三区av精品| 久久精品国产亚洲av涩爱| 国产黄片美女视频| 三级男女做爰猛烈吃奶摸视频| 色综合色国产| 亚洲国产精品合色在线| 乱人视频在线观看| 精品午夜福利在线看| 日韩人妻高清精品专区| 在线观看66精品国产| 国产综合懂色| 亚洲国产精品国产精品| 青春草视频在线免费观看| 麻豆一二三区av精品| 国产精品.久久久| 在线a可以看的网站| ponron亚洲| 亚洲av中文av极速乱| 亚洲av中文av极速乱| 亚洲精品乱码久久久久久按摩| 国产亚洲5aaaaa淫片| 亚洲精品乱码久久久久久按摩| 亚洲av中文av极速乱| 日韩制服骚丝袜av| 99在线人妻在线中文字幕| 亚洲精品国产成人久久av| 天美传媒精品一区二区| 99久久无色码亚洲精品果冻| 日韩高清综合在线| 久久国内精品自在自线图片| 51国产日韩欧美| 特大巨黑吊av在线直播| 九九爱精品视频在线观看| 欧美3d第一页| 高清毛片免费看| 国产精品麻豆人妻色哟哟久久 | 国产色婷婷99| 99久久精品一区二区三区| 国产伦一二天堂av在线观看| 亚洲综合精品二区| 一级毛片我不卡| 建设人人有责人人尽责人人享有的 | 亚洲av二区三区四区| 我的老师免费观看完整版| 日产精品乱码卡一卡2卡三| 亚洲国产精品sss在线观看| 国产一区二区在线观看日韩| 少妇的逼好多水| 一级毛片久久久久久久久女| 女人十人毛片免费观看3o分钟| 亚洲欧美清纯卡通| 亚洲四区av| 超碰97精品在线观看| 亚洲图色成人| 99热精品在线国产| 伦理电影大哥的女人| 国产午夜福利久久久久久| 免费观看a级毛片全部| 亚洲国产欧美在线一区| 舔av片在线| 亚洲欧美清纯卡通| 一级黄片播放器| 欧美高清性xxxxhd video| 国产日韩欧美在线精品| 伊人久久精品亚洲午夜| 国产女主播在线喷水免费视频网站 | 成年免费大片在线观看| 尤物成人国产欧美一区二区三区| 欧美变态另类bdsm刘玥| 精品国产一区二区三区久久久樱花 | 免费在线观看成人毛片| 99热全是精品| 精品人妻一区二区三区麻豆| 91精品一卡2卡3卡4卡| 国产老妇女一区| av播播在线观看一区| 99在线人妻在线中文字幕| 日日撸夜夜添| 国产视频首页在线观看| videossex国产| 特大巨黑吊av在线直播| 一级毛片aaaaaa免费看小| 国产伦在线观看视频一区| 你懂的网址亚洲精品在线观看 | 国内精品宾馆在线| 欧美3d第一页| 男女那种视频在线观看| 少妇熟女aⅴ在线视频| 身体一侧抽搐| 国产 一区精品| 床上黄色一级片| 搡老妇女老女人老熟妇| 插阴视频在线观看视频| 麻豆乱淫一区二区| 又黄又爽又刺激的免费视频.| 伊人久久精品亚洲午夜| 日本黄色视频三级网站网址| 久久99精品国语久久久| 国产精品一及| 欧美性感艳星| 午夜老司机福利剧场| 18禁动态无遮挡网站| 日韩欧美在线乱码| h日本视频在线播放| 亚洲精品乱码久久久久久按摩| 久久久成人免费电影| 我要看日韩黄色一级片| 免费观看的影片在线观看| 男人狂女人下面高潮的视频| av在线老鸭窝| 国内精品一区二区在线观看| 亚洲激情五月婷婷啪啪| 亚洲伊人久久精品综合 | 老司机影院成人| 成人亚洲精品av一区二区| 亚洲在线观看片| .国产精品久久| 人人妻人人澡欧美一区二区| 国产高清视频在线观看网站| 久久精品久久精品一区二区三区| 乱系列少妇在线播放| 色尼玛亚洲综合影院| 国产精品人妻久久久久久| 欧美三级亚洲精品| 在线观看一区二区三区| 国产成人精品婷婷| 日韩欧美精品免费久久| 精品久久久久久成人av| 国产成人午夜福利电影在线观看| 国语自产精品视频在线第100页| 只有这里有精品99| 国产精品日韩av在线免费观看| 亚洲国产色片| 乱码一卡2卡4卡精品| 国产午夜精品论理片| 边亲边吃奶的免费视频| 亚洲成人精品中文字幕电影| 丰满乱子伦码专区| 麻豆成人午夜福利视频| 特大巨黑吊av在线直播| 欧美日本视频| 天天躁夜夜躁狠狠久久av| 亚洲一级一片aⅴ在线观看| 狠狠狠狠99中文字幕| 干丝袜人妻中文字幕| 熟女电影av网| 免费看a级黄色片| 国产高清视频在线观看网站| 亚洲av免费高清在线观看| 日韩视频在线欧美| 日本-黄色视频高清免费观看| 亚洲国产精品专区欧美| 老司机影院成人| 97超视频在线观看视频| 淫秽高清视频在线观看| 韩国av在线不卡| 黄色配什么色好看| 日本wwww免费看| 国产精品熟女久久久久浪| 免费观看精品视频网站| 午夜福利在线在线| 免费看美女性在线毛片视频| 久久人人爽人人片av| 亚洲最大成人中文| 日韩三级伦理在线观看| 乱系列少妇在线播放| 久久久久久久久久成人| 日韩制服骚丝袜av| 亚洲欧美清纯卡通| 精品久久久久久久人妻蜜臀av| 国产精品福利在线免费观看| 国产精品爽爽va在线观看网站| 午夜福利成人在线免费观看| 全区人妻精品视频| 国产成人精品婷婷| 非洲黑人性xxxx精品又粗又长| kizo精华| 久久欧美精品欧美久久欧美| 久久久成人免费电影| 欧美变态另类bdsm刘玥| 亚洲精品aⅴ在线观看| 日韩大片免费观看网站 | videossex国产| 国产91av在线免费观看| 麻豆成人av视频| 亚洲乱码一区二区免费版| 成人漫画全彩无遮挡| 亚洲一区高清亚洲精品| 不卡视频在线观看欧美| 免费看a级黄色片| 乱人视频在线观看| 麻豆成人午夜福利视频| 国产伦精品一区二区三区四那| 神马国产精品三级电影在线观看| 国产视频首页在线观看| 精品无人区乱码1区二区| 麻豆成人av视频| 国产精品福利在线免费观看| 精品一区二区免费观看| 日日摸夜夜添夜夜爱| 亚洲欧美精品综合久久99| 欧美+日韩+精品| 亚洲国产精品sss在线观看| 视频中文字幕在线观看| 国产亚洲一区二区精品| 午夜视频国产福利| 自拍偷自拍亚洲精品老妇| 男女国产视频网站| 亚洲欧美中文字幕日韩二区| 成人一区二区视频在线观看| 午夜福利在线观看免费完整高清在| 国产精品日韩av在线免费观看| 人妻系列 视频| 精品少妇黑人巨大在线播放 | 欧美成人精品欧美一级黄| 国产精品国产高清国产av| 日本色播在线视频| av卡一久久| 色尼玛亚洲综合影院| 美女脱内裤让男人舔精品视频| 久久这里只有精品中国| 97热精品久久久久久| 国产综合懂色| 久久精品人妻少妇| 噜噜噜噜噜久久久久久91| 久久精品国产亚洲网站| 色综合色国产| 日韩大片免费观看网站 | kizo精华| www日本黄色视频网| 最后的刺客免费高清国语| 亚洲天堂国产精品一区在线| 两个人的视频大全免费| 一级毛片电影观看 | 欧美激情久久久久久爽电影| 亚洲美女视频黄频| 春色校园在线视频观看| 波多野结衣巨乳人妻| 亚洲精品久久久久久婷婷小说 | 老司机福利观看| 精品不卡国产一区二区三区| 国产成年人精品一区二区| 免费黄色在线免费观看| 国产午夜精品一二区理论片| 中文欧美无线码| 日韩av在线大香蕉| 国产色爽女视频免费观看| 中文字幕av成人在线电影| 国产精品三级大全| 国产亚洲一区二区精品| 麻豆成人av视频| 国内揄拍国产精品人妻在线| 好男人在线观看高清免费视频| 看片在线看免费视频| 日韩人妻高清精品专区| 一区二区三区四区激情视频| 尤物成人国产欧美一区二区三区| 亚洲av日韩在线播放| 亚洲欧洲日产国产| 亚洲人与动物交配视频| 国产精品99久久久久久久久| 男女下面进入的视频免费午夜| 亚洲综合精品二区| 国产精品不卡视频一区二区| 亚洲国产精品sss在线观看| 能在线免费看毛片的网站| 日韩视频在线欧美| 91av网一区二区| 亚洲精品影视一区二区三区av| 一级爰片在线观看| 国产乱人视频| 99国产精品一区二区蜜桃av| 国产综合懂色| 亚洲欧美精品自产自拍| 97超视频在线观看视频| 国产美女午夜福利| 国产探花极品一区二区| 国产精品久久久久久久电影| 99在线人妻在线中文字幕| 免费看日本二区| 欧美区成人在线视频| 色综合站精品国产| 一级av片app| 一级黄色大片毛片| 亚洲18禁久久av| 麻豆精品久久久久久蜜桃| www.av在线官网国产| 亚洲精品乱久久久久久| 国产黄色视频一区二区在线观看 | 国产亚洲午夜精品一区二区久久 | 3wmmmm亚洲av在线观看| 最近手机中文字幕大全| 成人美女网站在线观看视频| 国产真实伦视频高清在线观看| 丰满少妇做爰视频| 看片在线看免费视频| 深夜a级毛片| 亚洲欧美清纯卡通| 97超碰精品成人国产| 哪个播放器可以免费观看大片| 99久久精品国产国产毛片| 国产亚洲精品av在线| 久久婷婷人人爽人人干人人爱| 三级毛片av免费| 色5月婷婷丁香| 老司机福利观看| 国产精品久久久久久久电影| 欧美日韩精品成人综合77777| 亚洲欧美日韩高清专用| 国产精品久久久久久精品电影小说 | 国产色婷婷99| 三级男女做爰猛烈吃奶摸视频| 久久这里有精品视频免费| 国产美女午夜福利| 欧美极品一区二区三区四区| 欧美区成人在线视频| 欧美精品国产亚洲| 国语对白做爰xxxⅹ性视频网站| 女的被弄到高潮叫床怎么办| 尾随美女入室| 亚洲av中文av极速乱| 日韩中字成人| 久久久久久久久大av| 亚洲综合精品二区| 免费一级毛片在线播放高清视频| 免费看a级黄色片| 女人十人毛片免费观看3o分钟| 欧美丝袜亚洲另类| 中文字幕制服av| 亚洲天堂国产精品一区在线| 久久久久久久久久久免费av| 日产精品乱码卡一卡2卡三| 免费大片18禁| 男人和女人高潮做爰伦理| 99久久精品一区二区三区| 麻豆久久精品国产亚洲av| 两个人视频免费观看高清| 免费黄网站久久成人精品| 夜夜看夜夜爽夜夜摸| 天堂√8在线中文| 熟女人妻精品中文字幕| 亚洲成av人片在线播放无| 久久婷婷人人爽人人干人人爱| 国产真实乱freesex| 97超碰精品成人国产| 亚洲精品亚洲一区二区| 成人午夜高清在线视频| 免费看av在线观看网站| 欧美丝袜亚洲另类| 久久精品熟女亚洲av麻豆精品 | 草草在线视频免费看| 亚洲aⅴ乱码一区二区在线播放| a级毛色黄片| 午夜视频国产福利| 久久99精品国语久久久| 狂野欧美激情性xxxx在线观看| 人妻少妇偷人精品九色| 国产精品麻豆人妻色哟哟久久 | 日本免费一区二区三区高清不卡| 亚洲精品国产成人久久av| 久久精品国产亚洲网站| 在线播放无遮挡| 91狼人影院| 国产一区二区三区av在线| 美女内射精品一级片tv| 国产成人一区二区在线| 观看美女的网站| 国产午夜精品论理片| 免费看美女性在线毛片视频| 久久鲁丝午夜福利片| 亚洲av熟女| 国语自产精品视频在线第100页| 一个人观看的视频www高清免费观看| 成年版毛片免费区| 一本久久精品| 22中文网久久字幕| 国产美女午夜福利| 欧美日韩在线观看h| 国产精品久久久久久久久免| 亚洲美女视频黄频| 国产精品麻豆人妻色哟哟久久 | 午夜亚洲福利在线播放| 伦理电影大哥的女人| 好男人视频免费观看在线| 美女脱内裤让男人舔精品视频| 亚洲精品国产av成人精品| 国产一级毛片七仙女欲春2| 亚洲欧美日韩高清专用| 亚洲精品久久久久久婷婷小说 | 午夜福利高清视频| 99热这里只有是精品在线观看| 美女大奶头视频| 一卡2卡三卡四卡精品乱码亚洲| 91aial.com中文字幕在线观看| 欧美极品一区二区三区四区| 免费不卡的大黄色大毛片视频在线观看 | 丰满少妇做爰视频| 特级一级黄色大片| 亚洲18禁久久av| 如何舔出高潮| 国产亚洲av嫩草精品影院| 婷婷色av中文字幕| 内射极品少妇av片p| 汤姆久久久久久久影院中文字幕 | 久久精品91蜜桃| 亚洲欧洲日产国产| 亚洲欧美精品综合久久99| av在线蜜桃| 国产高清有码在线观看视频| 男女那种视频在线观看| 村上凉子中文字幕在线| 久久午夜福利片| 赤兔流量卡办理| 色吧在线观看| 国产一级毛片在线| 亚洲精品aⅴ在线观看| 亚洲av一区综合| 高清在线视频一区二区三区 | 日本午夜av视频| 久久精品国产亚洲av天美| 爱豆传媒免费全集在线观看| 欧美高清成人免费视频www| 一夜夜www| 夫妻性生交免费视频一级片| 国产淫片久久久久久久久| 久久久精品欧美日韩精品| 最近最新中文字幕免费大全7| 国产亚洲一区二区精品| 亚洲中文字幕一区二区三区有码在线看| 嫩草影院精品99| 真实男女啪啪啪动态图| www日本黄色视频网| 国产中年淑女户外野战色| 你懂的网址亚洲精品在线观看 | 日本一本二区三区精品| 国产极品天堂在线| 国产老妇伦熟女老妇高清| 女人久久www免费人成看片 | 久久久久久久久久成人| 亚洲欧洲日产国产| 中文字幕久久专区| 日本wwww免费看| 99热全是精品| 久久精品人妻少妇| 18禁动态无遮挡网站| 色吧在线观看| 91精品一卡2卡3卡4卡| 毛片一级片免费看久久久久| 成人毛片a级毛片在线播放| 久久久久性生活片| 2021天堂中文幕一二区在线观| 一卡2卡三卡四卡精品乱码亚洲| 国产成年人精品一区二区| 欧美日本视频| 精品久久久久久久末码| 亚洲激情五月婷婷啪啪| 国产v大片淫在线免费观看| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 乱人视频在线观看| 一级黄色大片毛片| 少妇人妻精品综合一区二区| 亚洲欧洲国产日韩| 久久这里有精品视频免费| 亚洲欧美日韩卡通动漫| 1000部很黄的大片| 中文字幕制服av| АⅤ资源中文在线天堂| 亚洲,欧美,日韩| 最后的刺客免费高清国语| 神马国产精品三级电影在线观看| 高清毛片免费看| 97在线视频观看| 免费看光身美女| 日本欧美国产在线视频| 精品人妻视频免费看| 国产又黄又爽又无遮挡在线| 欧美一区二区亚洲| 精品酒店卫生间| av女优亚洲男人天堂| 精华霜和精华液先用哪个| 国产色婷婷99| 久久6这里有精品| 国产免费福利视频在线观看| 美女高潮的动态| 久久久久久久久久黄片| 日韩在线高清观看一区二区三区| 卡戴珊不雅视频在线播放| 欧美精品一区二区大全| 久久久久久久国产电影| 成年女人永久免费观看视频| 亚洲成色77777| 午夜激情欧美在线| 亚洲av.av天堂| 国产亚洲精品av在线| 日韩欧美精品免费久久| av在线观看视频网站免费| 乱码一卡2卡4卡精品| 青春草亚洲视频在线观看| .国产精品久久| 久久精品国产亚洲网站| 国产精品国产三级国产av玫瑰| 国产伦一二天堂av在线观看| 亚洲欧美清纯卡通| 亚洲成人精品中文字幕电影| 亚洲人成网站高清观看| 日日干狠狠操夜夜爽| 日本黄大片高清| 一级av片app| 高清av免费在线| 网址你懂的国产日韩在线| 在线免费观看不下载黄p国产| 久久婷婷人人爽人人干人人爱| 亚洲精华国产精华液的使用体验| 丰满少妇做爰视频| 国产一区二区三区av在线| 亚洲三级黄色毛片| 亚洲av电影在线观看一区二区三区 | 能在线免费观看的黄片| 国产精品av视频在线免费观看| 日韩一区二区视频免费看| 亚洲精品影视一区二区三区av| 国产高潮美女av| 九九爱精品视频在线观看| 噜噜噜噜噜久久久久久91| 插逼视频在线观看| 婷婷色av中文字幕| 精品久久国产蜜桃| .国产精品久久| 成年免费大片在线观看| 亚洲av免费在线观看| 国内精品一区二区在线观看| 成人av在线播放网站| 亚洲美女视频黄频| 色5月婷婷丁香| 99在线视频只有这里精品首页| 26uuu在线亚洲综合色| 免费一级毛片在线播放高清视频| 亚洲国产最新在线播放| 亚洲乱码一区二区免费版| 国产色婷婷99| 美女内射精品一级片tv| 国产午夜精品久久久久久一区二区三区| 亚洲在线自拍视频| 一区二区三区四区激情视频| 亚洲中文字幕一区二区三区有码在线看| 亚洲av电影不卡..在线观看| 欧美性猛交黑人性爽| 欧美日韩综合久久久久久| 干丝袜人妻中文字幕| 精品午夜福利在线看| 国产精品.久久久| 麻豆精品久久久久久蜜桃| 亚州av有码| 日韩国内少妇激情av| 午夜爱爱视频在线播放| 免费av毛片视频| 国产伦理片在线播放av一区| 日韩高清综合在线| 中文欧美无线码| 国产免费又黄又爽又色| 卡戴珊不雅视频在线播放| 国产三级在线视频| 可以在线观看毛片的网站| 国模一区二区三区四区视频| 亚洲精品日韩av片在线观看| 欧美97在线视频| 精品久久国产蜜桃| 伦理电影大哥的女人| 国内精品宾馆在线| 国产三级在线视频|