CN102916922B - Adaptive search Doppler compensation method for underwater sound OFDM - Google Patents

Adaptive search Doppler compensation method for underwater sound OFDM Download PDF

Info

Publication number
CN102916922B
CN102916922B CN201210388724.9A CN201210388724A CN102916922B CN 102916922 B CN102916922 B CN 102916922B CN 201210388724 A CN201210388724 A CN 201210388724A CN 102916922 B CN102916922 B CN 102916922B
Authority
CN
China
Prior art keywords
doppler
mse
data
factor
signal
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201210388724.9A
Other languages
Chinese (zh)
Other versions
CN102916922A (en
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.)
Nanhai innovation and development base of Sanya Harbin Engineering University
Original Assignee
Harbin Engineering 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 Harbin Engineering University filed Critical Harbin Engineering University
Priority to CN201210388724.9A priority Critical patent/CN102916922B/en
Publication of CN102916922A publication Critical patent/CN102916922A/en
Application granted granted Critical
Publication of CN102916922B publication Critical patent/CN102916922B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Abstract

The invention aims to provide an adaptive search Doppler compensation method for an underwater sound OFDM (Orthogonal Frequency Division Multiplexing). The method comprises the following steps: adopting CW (Continuous Wave) signals as training sequences to roughly measure Doppler frequency deviation factors; utilizing high-power DFT to compensate the Doppler frequency deviation; taking the Doppler frequency deviation factors measured by the CW signals as original values; taking mean-square error of compensated and balanced data as a price function; and unceasingly regulating the size of the factor to perform search until the conditions are met. According to the method provided by the invention, the defect that the Doppler frequency deviation factors measured by a current block Doppler estimation or single-frequency signal frequency measurement method have error, as well as the defect that the divergence of a data star map, namely the increasing of data mean-square error, is caused by directly utilizing the factors to compensate Doppler frequency deviation, is avoided.

Description

Underwater sound OFDM adaptable search Doppler Compensation Method
Technical field
What the present invention relates to is a kind of compensation method of high-speed underwater sound communication.
Background technology
Sound wave is the information carrier of carried out remote transmission unique under water, along with deepening continuously of ocean development, to the efficiency of transmission of mobile communication under water and reliability requirement also more and more higher.Orthogonal frequency division multiplexi has the advantage such as anti-multipath fading, anti-bursty interference, and is suitable for high speed data transfer.OFDM symbol is actually multiple sub-carrier signal sum through ovennodulation, utilizes orthogonality to distinguish between each subcarrier.If but the interference existed between subcarrier, this orthogonality will be destroyed, cause systematic function to decline.Cause the reason of inter-sub-carrier interference to have multiple, one of them is exactly that receiving terminal and transmitting terminal relative motion cause Doppler shift, destroys orthogonality between subcarrier, causes demodulation mistake.So carrying out Doppler effect correction for OFDM System for Underwater Acoustic is the key improving systematic function.
Existing OFDM System for Underwater Acoustic Doppler effect correction technology, many employing block Doppler estimate or simple signal frequency measuring method estimating Doppler frequency offset factor, are then compensated Doppler by Variable sampling method.It is respectively insert a known signal in the front and back of Frame that block Doppler estimates, at receiving terminal by carrying out with local signal the width that relevant treatment records receiving data frames, compare with transmission data frame length and record the Doppler shift factor, its certainty of measurement, by the restriction of the ambiguity figure of known signal, namely requires that the time-bandwidth product of signal is enough large.Simple signal frequency measuring method is by measuring the CW signal frequency received, and compare with the CW signal frequency sent and record the Doppler shift factor, its certainty of measurement by the restriction of CW signal pulse length, and is vulnerable to noise and multi-path jamming.
Summary of the invention
The object of the present invention is to provide the underwater sound OFDM adaptable search Doppler Compensation Method not being vulnerable to noise and multi-path jamming.
The object of the present invention is achieved like this:
Underwater sound OFDM adaptable search Doppler Compensation Method of the present invention, is characterized in that:
(1) CW signal is adopted to carry out the bigness scale of the Doppler shift factor as training sequence: definition Doppler shift factor mu d=v/c, wherein v represents the speed of related movement of transmitting terminal and receiving terminal, and c is the speed of sound wave at water transmission, and the CW signal frequency of transmission is f t, the CW signal frequency of reception is f r, first filtering is carried out to the CW signal received, then Fourier transform is done to it, then record CW signal frequency by curve-parabola-fitting method and be the then Doppler shift factor of bigness scale
(2) high power DFT is utilized to compensate Doppler shift: the OFDM symbol length of transmission is that N length N does not comprise that circulation is forward and backward sews length, because signal is subject to the impact of Doppler shift, make the time domain waveform of OFDM symbol that compression or expansion occur, namely the sampled point of a complete ofdm signal waveform becomes N ' point, does 2 after intercepting at N ' to the ofdm signal received npoint DFT, recalculates the position of data subcarrier according to the Doppler shift factor recorded;
(3) the Doppler shift factor mu gone out by CW Signal estimation, as initial condition, is searched for using the MSE of the reception data after equilibrium as cost function, until the MSE calculated can find minimum:
Before intercepting OFDM time-domain symbol, the frequency offset factor μ of bigness scale is utilized to recalculate OFDM symbol length, that is:
N ′ = [ N 1 + μ ] ,
Wherein [] represents round, high power DFT Doppler effect correction is carried out at N ' after again intercepted to OFDM symbol, OFDM symbol adopts Comb Pilot structure, Comb Pilot interpolation is utilized to carry out channel estimating, equalized data subcarriers, because the frequency offset factor of bigness scale exists error, Doppler effect correction with balanced after the mean square error of data comparatively large, so will search procedure be carried out.The mean square error now calculated is stored:
MSE(0)=MSE μ
Then to Doppler factor with step size mu fcarry out new Doppler factor search, obtain new Doppler factor:
μ k=μ+kμ f,k=±1,±2,...,±K
By μ ksubstitution formula calculate the OFDM symbol length made new advances, and again intercept reception data, Doppler effect correction afterwards and channel equalization process repeat above-mentioned steps, and the data mean square error calculated are stored, that is:
MSE ( k ) = MSE μ k ,
Make k=± 1, ± 2 ..., ± K, increases μ kabsolute value search for, and compare the MSE value calculated, if can find the minimum of MSE, then the Doppler shift factor of its correspondence is estimated correct, jumps out circulation, exports data, completes channel disturbance and suppresses; If the minimum of MSE can not be found, then repeat said process, until find the minimum of MSE then to jump out circulation, export data, complete channel disturbance and suppress.
Advantage of the present invention is: existing piece of Doppler estimates or simple signal frequency measuring method, owing to being subject to the restriction of signal itself, and noise and multi-path jamming, all there is certain error in the Doppler shift factor recorded.If directly utilize its factor recorded to carry out Doppler shift compensation, dispersing of data constellation figure can be caused, i.e. the increase of data mean square error.The present invention is using the Doppler factor of CW Signal estimation as initial value, with the mean square error of the data after compensation, equilibrium for cost function, continuous adjustment Doppler shift factor size is searched for, until the data mean square error calculated is enough little, then data have converged in each phase place of transmission planisphere, can stop search.
Accompanying drawing explanation
Fig. 1 is OFDM water sound communication signal frame structure of the present invention;
Fig. 2 is high power DFT method Doppler compensate of frequency deviation of the present invention;
Fig. 3 is adaptable search Doppler effect correction flow chart of the present invention;
Fig. 4 is the MSE of data after different frequency offset factor of the present invention compensates;
Fig. 5 is the BER of data after different frequency offset factor of the present invention compensates;
Fig. 6 is Doppler shift factor bigness scale compensation result of the present invention;
Fig. 7 is the Doppler shift factor of the present invention search compensation result.
Embodiment
Below in conjunction with accompanying drawing citing, the present invention is described in more detail:
Composition graphs 1 ~ 7, adopts CW signal to carry out the bigness scale of the Doppler shift factor as training sequence.High power DFT is utilized to compensate Doppler shift.The Doppler shift factor recorded using CW signal is as initial value, and the mean square error of the data to compensate, after equilibrium is cost function, and continuous Dynamic gene size is searched for, until satisfy condition.
Modelled signal frame structure as shown in Figure 1.Time synchronizing signal LFM is used for carrying out frame synchronization and timing synchronization, can find frequency deviation estimated signal CW by relevant treatment.Measure the CW signal frequency of reception again, compared with transmitting terminal CW signal frequency, the rough estimate of the Doppler shift factor can be obtained.
Doppler effect correction takes high power DFT method, does not comprise circulation forward and backwardly sew length if the OFDM symbol length sent is N(), the OFDM symbol length after affecting by Doppler shift is N '.N ' the point data of intercepting is done 2 npoint DFT, by increasing the method that DFT counts, increasing frequency resolution, correctly extracting data message after recalculating data subcarrier position.High power DFT method Doppler compensate of frequency deviation as shown in Figure 2.
The Doppler shift factor mu gone out by CW Signal estimation is as initial value, using the mean square error (MSE) of the reception data after high power DFT Doppler effect correction, channel equalization as cost function, the size of the continuous adjustment Doppler shift factor is searched for, until the MSE calculated stops search after satisfying condition.Flow chart as shown in Figure 3.
1, CW signal Doppler shift factor bigness scale
Adopt the signal frame structure shown in Fig. 1.Time synchronizing signal LFM is used for carrying out frame synchronization and timing synchronization, can find frequency deviation estimated signal CW by relevant treatment.Definition Doppler shift factor mu d=v/c, wherein v represents the speed of related movement of transmitting terminal and receiving terminal, and c is the speed of sound wave at water transmission.The CW signal frequency sent is f t, the CW signal frequency of reception is f r.First filtering is carried out to the CW signal received, then Fourier transform is done to it, then record CW signal frequency by curve-parabola-fitting method and be the then Doppler shift factor of bigness scale due to what record with the f of reality rthere is certain deviation, then also there is certain error in the Doppler shift factor of bigness scale.
2, high power DFT Doppler shift compensates
High power DFT Doppler shift compensates flow chart as shown in Figure 2.Do not comprise circulation forward and backwardly sew length if the OFDM symbol length sent is N().Because signal is subject to the impact of Doppler shift, make the time domain waveform of OFDM symbol that compression or expansion occur, namely the sampled point of a complete ofdm signal waveform becomes N ' point.2 are done at N ' after intercepted to the ofdm signal received npoint DFT, being then equivalent to carrying out zero padding after OFDM symbol, adding frequency resolution.Recalculate the position of data subcarrier according to the Doppler shift factor recorded, can data message be extracted.
3, adaptable search Doppler effect correction
As previously described, because CW signal is vulnerable to the impact of multi-path jamming and noise jamming, and frequency-measurement accuracy is limited to pulsewidth length to a certain extent, therefore can error be there is in its Doppler shift factor mu recorded, if directly utilize μ to carry out Doppler effect correction to reception data, then the planisphere of data there will be Divergent Phenomenon, is unfavorable for correct demodulation.If but correctly estimated Doppler shift factor mu, the reception data constellation figure after compensation can well be converged in send planisphere each phase place on, the error rate is corresponding reduction also.For QPSK mapping mode, then the phase place sending data is amplitude is 1, lays respectively in the middle of four quadrants.If R k, k=0,1 ... n-1 is the data point being distributed in same quadrant, and quantity is n, and lower footnote k is subcarrier ordinal number.In order to weigh the degree of divergence of data constellation figure, the mean square error MSE (dB) defining the data of each quadrant is:
MSE = 10 lg ( 1 n Σ k = 0 n - 1 ( R k - R ‾ ) ( R k - R ‾ ) * ) - - - ( 1 )
Wherein represent the mean value of the data subcarrier in a quadrant, namely symbol * represents complex conjugate operation; The inverse exponent arithmetic that it is the end that lg represents with 10.In order to the mean square error of verification msg and the proportional relation of the error rate, under underwater sound multipath channel, ofdm system is emulated.If receiving terminal and transmitting terminal counter motion, Doppler shift factor mu d=-1.7 × 10 -4.After different frequency offset factor compensates, as shown in Figure 4, the error rate as shown in Figure 5 for the mean square error of data.As can be seen from Fig. 4, Fig. 5, along with Doppler shift factor mu dchange, MSE and the BER variation tendency of OFDM symbol data is substantially identical, all at μ dclose to reaching minimum during theoretical value.In addition, the minimum point of MSE and the BER of data is not quite affected by noise, and the frequency offset factor that namely minimum is corresponding under low signal-to-noise ratio is still near theoretical value.
Based on above discussion, consider that the Doppler shift factor mu gone out by CW Signal estimation is as initial condition, constantly search for as cost function, until the MSE calculated stops search after satisfying condition using the mean square error (MSE) of the reception data after equilibrium.Flow chart as shown in Figure 3.
Wherein time synchronized is identical with aforesaid method with frequency deviation bigness scale process.Before intercepting OFDM time-domain symbol, the frequency offset factor μ of bigness scale be utilized to recalculate OFDM symbol length, that is:
N ′ = [ N 1 + μ ] - - - ( 2 )
Wherein [] represents round.High power DFT Doppler effect correction can be carried out at N ' after again intercepted to OFDM symbol.OFDM symbol adopts Comb Pilot structure, therefore Comb Pilot interpolation can be utilized to carry out channel estimating, equalized data subcarriers.Because the frequency offset factor of bigness scale exists error, Doppler effect correction with balanced after the mean square error of data comparatively large, so will search procedure be carried out.The mean square error now calculated is stored:
MSE(0)=MSE μ (3)
Then to Doppler factor with step size mu fcarry out new Doppler factor search, obtain new Doppler factor:
μ k=μ+kμ f,k=±1,±2,...,±K (4)
By μ ksubstitution formula (2) calculates the OFDM symbol length made new advances, and again intercepts reception data.Doppler effect correction afterwards and channel equalization process same as above, and the data mean square error calculated to be stored, that is:
MSE ( k ) = MSE μ k - - - ( 5 )
Make k=± 1, ± 2 ..., ± K, constantly increases μ kabsolute value search for, and compare the MSE value calculated, if the minimum of MSE can be found, then think that the Doppler shift factor of its correspondence is estimated correct, can circulation be jumped out, export data, complete channel disturbance and suppress; If do not met, then repeat said process, until find the minimum of MSE.
Below self-adaptive search algorithm is emulated.Still receiving terminal and transmitting terminal counter motion is supposed, Doppler shift factor mu d=-1.667 × 10 -4.OFDM symbol subcarrier adopts QPSK modulation system, and Comb Pilot carries out channel estimating, and signal to noise ratio is 10dB.First utilize the CW signal of frame structure front portion to carry out Doppler shift bigness scale, obtain frequency offset factor and estimate μ=-2.2 × 10 -4.Compensate with it and receive data and obtain data constellation figure as shown in Figure 6 after carrying out channel equalization.The mean value calculating four quadrant mean square errors after data point being divided by quadrant is-0.692dB.
Data constellation figure disperses comparatively serious as can be seen from Figure 6, illustrates that the frequency offset factor of bigness scale is inaccurate.Below by step size mu f=0.5 × 10 -4search for, result as shown in Figure 7.As seen from Figure 7, MSE corresponding during k=+1 is minimum, i.e. μ k=-1.7 × 10 -4for the final Doppler shift factor is estimated.

Claims (1)

1. underwater sound OFDM adaptable search Doppler Compensation Method, is characterized in that:
(1) CW signal is adopted to carry out the bigness scale of the Doppler shift factor as training sequence: definition Doppler shift factor mu d=v/c, wherein v represents the speed of related movement of transmitting terminal and receiving terminal, and c is the speed of sound wave at water transmission, and the CW signal frequency of transmission is f t, the CW signal frequency of reception is f r, first filtering is carried out to the CW signal received, then Fourier transform is done to it, then record CW signal frequency by curve-parabola-fitting method and be the then Doppler shift factor of bigness scale
(2) high power DFT is utilized to compensate Doppler shift: the OFDM symbol length of transmission is that N length N does not comprise that circulation is forward and backward sews length, because signal is subject to the impact of Doppler shift, make the time domain waveform of OFDM symbol that compression or expansion occur, namely the sampled point of a complete ofdm signal waveform becomes N ' point, does 2 after intercepting at N ' to the ofdm signal received npoint DFT, recalculates the position of data subcarrier according to the Doppler shift factor recorded;
(3) the Doppler shift factor mu gone out by CW Signal estimation, as initial condition, is searched for using the MSE of the reception data after equilibrium as cost function, until the MSE calculated can find minimum:
Before intercepting OFDM time-domain symbol, the frequency offset factor μ of bigness scale is utilized to recalculate OFDM symbol length, that is:
N ′ = [ N 1 + μ ] ,
Wherein [] represents round, high power DFT Doppler effect correction is carried out at N ' after again intercepted to OFDM symbol, OFDM symbol adopts Comb Pilot structure, Comb Pilot interpolation is utilized to carry out channel estimating, equalized data subcarriers, because the frequency offset factor of bigness scale exists error, Doppler effect correction with balanced after the mean square error of data comparatively large, so will search procedure be carried out; The mean square error now calculated is stored:
MSE(0)=MSE μ
Then to Doppler factor with step size mu fcarry out new Doppler factor search, obtain new Doppler factor:
μ k=μ+kμ f,k=±1,±2,...,±K
By μ ksubstitution formula calculate the OFDM symbol length made new advances, and reception data are intercepted again, Doppler effect correction afterwards and channel equalization process repeats step (2), and the data mean square error calculated is stored, that is:
MSE ( k ) = MSE μ k ,
Make k=± 1, ± 2 ..., ± K, increases μ kabsolute value search for, and compare the MSE value calculated, if can find the minimum of MSE, then the Doppler shift factor of its correspondence is estimated correct, jumps out circulation, exports data, completes channel disturbance and suppresses; If the minimum of MSE can not be found, then repeat said process, until find the minimum of MSE then to jump out circulation, export data, complete channel disturbance and suppress.
CN201210388724.9A 2012-10-15 2012-10-15 Adaptive search Doppler compensation method for underwater sound OFDM Active CN102916922B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201210388724.9A CN102916922B (en) 2012-10-15 2012-10-15 Adaptive search Doppler compensation method for underwater sound OFDM

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201210388724.9A CN102916922B (en) 2012-10-15 2012-10-15 Adaptive search Doppler compensation method for underwater sound OFDM

Publications (2)

Publication Number Publication Date
CN102916922A CN102916922A (en) 2013-02-06
CN102916922B true CN102916922B (en) 2014-12-17

Family

ID=47615157

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201210388724.9A Active CN102916922B (en) 2012-10-15 2012-10-15 Adaptive search Doppler compensation method for underwater sound OFDM

Country Status (1)

Country Link
CN (1) CN102916922B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107896204A (en) * 2017-10-11 2018-04-10 中国船舶重工集团公司第七〇五研究所 A kind of OFDM underwater sound Modem time-frequency two-dimensionals search Doppler compensator and compensation method based on FPGA

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103354538B (en) * 2013-07-20 2016-10-05 华南理工大学 A kind of method that Doppler effect correction is carried out to the receipt signal in underwater sound communication
CN106027116B (en) * 2016-07-07 2018-08-31 哈尔滨工程大学 A kind of mobile underwater sound communication Doppler coefficient method of estimation based on chirp signals
CN106100692A (en) * 2016-08-29 2016-11-09 东南大学 MIMO OFDM underwater sound communication system doppler spread method of estimation
CN106330251B (en) * 2016-08-29 2019-08-13 东南大学 Underwater sound communication system doppler spread estimation method based on zero correlation band sequence
CN108566354A (en) * 2018-04-03 2018-09-21 哈尔滨工程大学 DPFFT time-varying broadband Doppler Compensation Method in underwater sound OFDM

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102664840A (en) * 2012-04-26 2012-09-12 哈尔滨工程大学 Underwater sound OFDM (orthogonal frequency division multiplexing) Doppler estimation method based on cyclic prefixes
CN102724147A (en) * 2012-06-27 2012-10-10 哈尔滨工程大学 Channel estimation method for underwater sound orthogonal frequency division multiplexing

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102664840A (en) * 2012-04-26 2012-09-12 哈尔滨工程大学 Underwater sound OFDM (orthogonal frequency division multiplexing) Doppler estimation method based on cyclic prefixes
CN102724147A (en) * 2012-06-27 2012-10-10 哈尔滨工程大学 Channel estimation method for underwater sound orthogonal frequency division multiplexing

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107896204A (en) * 2017-10-11 2018-04-10 中国船舶重工集团公司第七〇五研究所 A kind of OFDM underwater sound Modem time-frequency two-dimensionals search Doppler compensator and compensation method based on FPGA

Also Published As

Publication number Publication date
CN102916922A (en) 2013-02-06

Similar Documents

Publication Publication Date Title
CN102916922B (en) Adaptive search Doppler compensation method for underwater sound OFDM
US7940848B2 (en) System having an OFDM channel estimator
US9020050B2 (en) Accounting for inter-carrier interference in determining a response of an OFDM communication channel
US6421401B1 (en) Method and apparatus for achieving and maintaining symbol synchronization particularly in an OFDM system
FI105963B (en) Procedure for forming a training period
JP2005287043A (en) Receiver
CN107231176B (en) OFDM-MFSK underwater acoustic communication broadband Doppler estimation and compensation method
CN102387115B (en) OFDM pilot scheme design and channel estimation method
JP2012028922A (en) Ofdm communication receiving device
CN107547143B (en) OFDM-MFSK underwater acoustic communication broadband Doppler estimation and compensation method with known subcarrier frequency
CN102404268A (en) Method for estimating and compensating doppler frequency offset in Rician channels in high-speed mobile environment
US8279743B2 (en) Method for interference estimation for orthogonal pilot patterns
CN107257324B (en) Time-frequency joint synchronization method and device in OFDM system
CN103312643B (en) Balancer, receiving system and equalization methods
CN102377726B (en) Timing synchronization method of OFDM (Orthogonal Frequency Division Multiplexing) system
CN101536440A (en) Subblock-wise frequency domain equalizer
CN102215184B (en) Method and system for estimating uplink timing error
CN103078819B (en) Fine symbol timing synchronization method and device thereof
CN101895487B (en) Confidence-based method and device for suppressing noises in channel estimation results
CN102957641B (en) Channel estimation method of TDS (Total Dissolved Solids)-OFDM (Orthogonal Frequency Division Multiplexing) system
CN111585740A (en) Transmission signal synchronization processing method, system, storage medium, program, and terminal
EP2704387B1 (en) SFO estimation technique for MIMO-OFDM frequency synchronization
TW200915751A (en) Orthogonal frequency division multiplexing system and its channel estimating appliance and method
KR100992369B1 (en) Apparatus for estimating channel of ofdm system
KR101063072B1 (en) Integer Frequency Error Estimation System and Method in WiBro System

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20201126

Address after: Area A129, 4th floor, building 4, Baitai Industrial Park, Yazhou Bay science and Technology City, Yazhou District, Sanya City, Hainan Province, 572024

Patentee after: Nanhai innovation and development base of Sanya Harbin Engineering University

Address before: 150001 Heilongjiang, Nangang District, Nantong street,, Harbin Engineering University, Department of Intellectual Property Office

Patentee before: HARBIN ENGINEERING University