Analys av algoritmer för detektering av resonansfrekvenser i vibrationsmätningar på överhettartuber

Institutionen för systemteknik

Department of Electrical Engineering

Examensarbete

Analys av algoritmer för detektering av

resonansfrekvenser i vibrationsmätningar på

överhettartuber

Examensarbete utfört i Reglerteknik vid Tekniska högskolan i Linköping

av

Daniel Eriksson

LiTH-ISY-EX--10/4371--SE Linköping 2010

Department of Electrical Engineering Linköpings tekniska högskola

Linköpings universitet Linköpings universitet

Analys av algoritmer för detektering av

resonansfrekvenser i vibrationsmätningar på

överhettartuber

Examensarbete utfört i Reglerteknik

vid Tekniska högskolan i Linköping

av

Daniel Eriksson

LiTH-ISY-EX--10/4371--SE

Handledare: Peter Rosander

isy, Linköpings Universitet

Elisabet Blom

Qring Technology AB Examinator: Svante Gunnarsson

isy, Linköpings Universitet Linköping, 23 juni, 2010

Avdelning, Institution

Division, Department

Division of Automatic Control Department of Electrical Engineering Linköpings universitet

SE-581 83 Linköping, Sweden

Datum Date 2010-06-23 Språk Language  Svenska/Swedish  Engelska/English  ⊠ Rapporttyp Report category  Licentiatavhandling  Examensarbete  C-uppsats  D-uppsats  Övrig rapport  ⊠

URL för elektronisk version

http://www.control.isy.liu.se http://urn.kb.se/resolve?urn=urn:nbn:se:liu:diva-57716 ISBN — ISRN LiTH-ISY-EX--10/4371--SE

Serietitel och serienummer

Title of series, numbering

ISSN

—

Titel

Title

Analys av algoritmer för detektering av resonansfrekvenser i vibrationsmätningar på överhettartuber

Analysis of algorithms for detection of resonance frequencies in vibration measu-rements on super heater tubes

Författare

Author

Daniel Eriksson

Sammanfattning

Abstract

Combustion in thermal power plants emits particles which create coatings on the super heater tubes. The coatings insulate the tubes and impairs the efficiency of the heat transfer. Cleaning the tubes occurs while the power plant is running and a system that automatically calculates the amount of coatings could make the cleaning more needs-based.

The resonance frequencies of the super heater tubes are affected by the added mass of the coatings. A change in frequency corresponds to a change in mass. Vibration measurements have been made with strain gauges on the super heater tubes in Ryaverket’s power plant and in one of SAKAB’s power plants.

The purpose of this thesis work has been to analyse different methods to esti-mate resonance frequencies in the vibration measurements. ESPRIT, MUSIC, AR and a heuristic statistical method have been tested on generated signals. MUSIC and ESPRIT have given the best estimations and have thus been used to analyse the measurements.

Periodically some estimations of measurements from Ryaverket are following trends which indicates that they could be resonance frequencies. The rest of the estimations contain too large variations. The estimations made of the measure-ments from SAKAB have been made during a shorter time period but shows clearer trends which make them probable resonance frequencies.

To automatically trace resonance frequencies in the estimations, even though they contains large variations, a target tracking algorithm has been implemented. The algorithm finds estimations that follows expected trends between the cleaning periods. Tests shows that the target tracking algorithm finds probable resonance frequencies in the estimations but that it is hard to reach a conclusion if they contain large variations. Better measurements could give estimations with smaller variations.

An idea is presented where MUSIC or ESPRIT together with a target tracking algorithm could be used to calculate the amount of coatings on the super heater tubes.

Nyckelord

Abstract

Combustion in thermal power plants emits particles which create coatings on the super heater tubes. The coatings insulate the tubes and impairs the efficiency of the heat transfer. Cleaning the tubes occurs while the power plant is running and a system that automatically calculates the amount of coatings could make the cleaning more needs-based.

The resonance frequencies of the super heater tubes are affected by the added mass of the coatings. A change in frequency corresponds to a change in mass. Vibration measurements have been made with strain gauges on the super heater tubes in Ryaverket’s power plant and in one of SAKAB’s power plants.

The purpose of this thesis work has been to analyse different methods to esti-mate resonance frequencies in the vibration measurements. ESPRIT, MUSIC, AR and a heuristic statistical method have been tested on generated signals. MUSIC and ESPRIT have given the best estimations and have thus been used to analyse the measurements.

Periodically some estimations of measurements from Ryaverket are following trends which indicates that they could be resonance frequencies. The rest of the estimations contain too large variations. The estimations made of the measure-ments from SAKAB have been made during a shorter time period but shows clearer trends which make them probable resonance frequencies.

To automatically trace resonance frequencies in the estimations, even though they contains large variations, a target tracking algorithm has been implemented. The algorithm finds estimations that follows expected trends between the cleaning periods. Tests shows that the target tracking algorithm finds probable resonance frequencies in the estimations but that it is hard to reach a conclusion if they contain large variations. Better measurements could give estimations with smaller variations.

An idea is presented where MUSIC or ESPRIT together with a target tracking algorithm could be used to calculate the amount of coatings on the super heater tubes.

Sammanfattning

På överhettartuber i värmepannor bildas beläggningar på grund av sot och partik-lar från förbränningen. Beläggningarna isolerar överhettartuberna vilket försämrar värmepannans verkningsgrad. Sotning av tuberna sker under drift och ett system som automatiskt kan beräkna hur mycket beläggningar som finns på tuberna skulle kunna göra sotningen mer behovsstyrd.

Resonansfrekvenser hos överhettartuberna påverkas av masspåslaget som blir när beläggningar bildas. En förändring i frekvens kan då översättas till en föränd-ring i massa. Vibrationsmätningar har gjorts med töjningsgivare som är monterade på överhettartuber i en av Ryaverkets pannor och i SAKAB:s panna.

I detta examensarbete har syftet varit att analysera olika metoder för att skatta resonansfrekvenser i genomförda mätningar. Algoritmerna MUSIC, ESPRIT och AR samt en heuristisk statistisk metod har testats på genererade signaler. MUSIC och ESPRIT har givit bäst skattningar och har därefter använts för att analysera mätningarna.

Periodvis följer vissa skattningar av mätningar från Ryaverket trender mellan sotningarna vilket indikerar att det skulle kunna vara resonansfrekvenser. Annars innehåller skattningarna för mycket variationer. Skattningarna av mätningar från SAKAB:s panna är gjorda under en kortare tidsperiod men visar en tydligare trend som troliggör att det är resonansfrekvenser som detekterats.

För att automatiskt hitta och följa resonansfrekvenser i skattningarna, även om dessa innehåller stora variationer, har en målföljningsalgoritm implementerats. Algoritmen hittar skattningar som följer en förväntad trend mellan sotningarna. Tester visar att algoritmen hittar troliga resonansfrekvenser i skattningarna men att det är svårt att kunna dra några slutsatser om skattningarna varierar för mycket. Bättre signaler skulle kunna minska variationerna hos skattningarna.

Ett förslag presenteras hur MUSIC eller ESPRIT tillsammans med en målfölj-ningsalgoritm skulle kunna användas för att beräkna masspåslag på överhettartu-berna.

Tack

Till att börja med vill jag tacka Elisabet Blom på Qring Technology AB och hennes medarbetare Claes Fredö och Lars Gabrielsson för förtroendet att få ge-nomföra detta examensarbete. Det har varit ett intressant och lärorikt arbete och jag uppskattar allt stöd som jag har fått.

Jag vill också tacka min examinator Svante Gunnarsson och min handledare Peter Rosander på ISY som har hjälpt mig under arbetets gång. Jag vill också tacka Niklas Wahlström som jag har suttit tillsammans med för alla givande diskussioner vi har haft.

Slutligen vill jag tacka Ylva Jung som hela tiden har stöttat mig och funnits vid min sida.

Daniel Eriksson, Linköping juni 2010

Innehåll

1 Introduktion 3

1.1 Parter . . . 3

1.1.1 Värmeforsk . . . 3

1.1.2 Qring Technology AB . . . 4

1.1.3 Borås energi och miljö AB . . . 4

1.1.4 SAKAB . . . 4 1.2 Bakgrund . . . 4 1.2.1 Fluidbäddpannor . . . 4 1.2.2 Mekaniska vibrationer . . . 5 1.2.3 Beläggningar . . . 6 1.2.4 Sotning . . . 6

1.2.5 Detektering av beläggningar på överhettartuber . . . 7

1.3 Mål . . . 9

1.4 Upplägg . . . 10

2 Skattningsmetoder 11 2.1 Fouriertransform och periodogram . . . 11

2.1.1 Welchs metod . . . 12

2.2 Fönstring . . . 12

2.3 Koherens . . . 12

2.4 Heuristisk metod för gradientanalys av frekvensinnehåll . . . 13

2.4.1 Manuell analysmetod . . . 14

2.4.2 Automatiserad statistisk metod . . . 14

2.5 MUSIC-algoritmen . . . 15 2.5.1 Parameterval för MUSIC . . . 17 2.6 ESPRIT-algoritmen . . . 17 2.6.1 Parameterval för ESPRIT . . . 18 2.7 AR-modellering . . . 18 2.7.1 Parameterval för AR . . . 19

3 Utvärdering av skattningsmetoder på simulerade data 21 3.1 Skattning av simulerade data med AR, MUSIC och ESPRIT . . . 21

3.1.1 Skattning av sinuskomponenter med varierande frekvens . . 22

3.1.2 Skattning av sinuskomponenter med olika amplitud . . . . 24 ix

3.1.3 Skattning av sinuskomponenter med varierande amplitud

och frekvens . . . 25

3.1.4 Skattning av en dämpad sinus . . . 27

3.2 Gradientanalys av frekvensinnehåll med automatiserad statistisk metod . . . 30

3.3 Diskussion . . . 31

4 Analys av mätdata 35 4.1 Datainsamling . . . 35

4.1.1 Givarnas placering . . . 35

4.1.2 Tider hos genomförda mätserier . . . 36

4.2 Knackprov . . . 39

4.3 Signalkvalitet . . . 40

4.3.1 Oönskat signalinnehåll . . . 40

4.3.2 Koherensanalys mellan givarna . . . 43

5 Resultat från skattningar 45 5.1 Analysstruktur . . . 45

5.2 Analys av signaler från Ryaverkets panna . . . 47

5.2.1 MUSIC och ESPRIT . . . 47

5.2.2 Periodogram . . . 50

5.3 Analys av signaler från SAKAB:s panna . . . 53

5.4 Diskussion . . . 56

6 Följning av resonansfrekvenser i skattningar från MUSIC och ESPRIT 57 6.1 Målföljning . . . 57

6.1.1 Kalmanfilter . . . 58

6.1.2 Associera data . . . 59

6.1.3 Initiering och terminering . . . 59

6.1.4 Algoritm . . . 60

6.2 Följning av resonansfrekvenser . . . 61

6.2.1 Anpassning av modell för resonansfrekvenser . . . 61

6.2.2 Parameterval . . . 63

6.2.3 Resultat med målföljningsalgoritmen . . . 65

6.3 Diskussion . . . 68

7 Slutsatser och fortsatt arbete 71 7.1 Slutsatser . . . 71

7.2 Användande av skattnings- och målföljningsalgoritmer . . . 72

7.2.1 MUSIC och ESPRIT . . . 72

7.2.2 Målföljningsalgoritm . . . 73

7.2.3 Användning av resonansfrekvenser för att förbättra sotningsru-tiner . . . 74

7.3 Fortsatt arbete . . . 74

Innehåll xi

A Implementerad MATLAB-kod 79

A.1 Skattning av resonansfrekvenser med MUSIC eller ESPRIT . . . . 79 A.1.1 Funktioner . . . 83 A.2 Målföljningsalgoritm . . . 84 A.2.1 Klasser . . . 85

Beteckningar

N (0, 1) Normalfördelning med väntevärde 0 och varians 1. AT _{Transponat till A.}

A−1

Invers till A. A∗

Komplexkonjugat och transponat av A. {zk}Nk=1 Mängden z1, z2, . . . , zN.

q−n Förskjutningsoperator q−n

z(t) = z(t − nT ). ˆ

z Skattning av z.

z(t + n) Värde för en diskret variabel, n sampel efter tidpunkt t.

ˆ

z(t|t − 1) Skattning av z vid tidpunkt t givet information om z fram till och med tidpunkt t − 1.

Ordlista

Egenvärde

Egenvärdena λ beräknas för en matris A som lösningarna till det (A − Iλ) = 0.

FFT

Fast Fourier Transform. Algoritm för att beräkna den diskreta fouriertransformen för en tidsdiskret signal.

Frekvenskomponent

I rapporten används begreppet frekvenskomponent i huvudsak för en svängning i en uppmätt signal med en viss frekvens. Metoderna som används i rapporten

benämns detektera frekvenskomponenter då det inte är klart att det är just reso-nansfrekvenser som metoderna detekterar.

Konfektionsmodell

Kallas även för svartlådemodell. Samlingsnamn för modeller som används för sy-stem med ett bestämt antal in- och utsignaler där information om fysikaliska samband som beskriver systemet saknas.

Minstakvadratanpassning

Minstakvadratanpassning av en parameter θ till en modell y = f (x, θ) utifrån mätdata (x1, y1), (x2, y2), . . ., (xn, yn) ges av lösningen till

θ = arg min θ n X i=1 (yi− f(xi, θ)) 2 (1) Nollrum

Nollrummet för en matris A är alla vektorer ¯x som uppfyller A¯x = ¯0 där ¯0 är nollvektorn.

Notchfilter

Ett filter designat för att filtrera bort signalinnehållet vid en bestämd frekvens.

Ortogonal

Två vektorer ¯x1och ¯x2, som båda har dimensionen n × 1 sägs vara ortogonala om ¯

xT

1x¯2= 0.

Rang

Rangen för en matris A är detsamma som antalet linjärt oberoende kolumner (eller rader) i A.

Störning

Ett oönskat signalinnehåll som kan ha känt eller okänt ursprung. Störningar går inte att påverka eller välja.

Kapitel 1

Introduktion

Detta examensarbete har gjorts på uppdrag från företaget Qring Technology AB som en del i ett projekt för forskningssamarbetsorganisationen Värmeforsk. Pro-jektet undersöker möjligheterna att behovsstyra sotning i fastbränsleeldade pannor genom att utnyttja tubernas mekaniska egenskaper.

En presentation av de olika parterna som ingår i detta arbete finns i avsnitt 1.1. I rapporten har fokus varit på att analysera mätningar av överhettartubernas vibrationer och därför finns i avsnitt 1.2 en kort beskrivning av den panntyp som finns i Ryaverket där mätningarna utförts. Vidare i samma avsnitt förklaras vad beläggning är och hur de påverkar pannans verkningsgrad. Målen och syftet med detta examensarbete beskrivs i avsnitt 1.3 och slutligen ges en sammanfattning av rapportens upplägg.

1.1

Parter

I detta avsnitt presenteras de företag och organisationer som berörts i detta arbete.

1.1.1

Värmeforsk

Värmeforsk är sedan 1968 ett gemensamt forskningsorgan för företag inom skogs-och tillverkningsindustri, myndigheter skogs-och energiproducenter. Enligt [18] bedriver organisationen i huvudsak forskning inom bränslebaserad kraft- och värmepro-duktion. Ett mål med forskningen inom Värmeforsk är att resultat ska ge kon-kreta lösningar som ska kunna tillämpas i industrin inom fem år. Organisationen fungerar också som grund för initiering av nya projekt, nätverksbildande inom medlemsföretagen och för spridning av kunskap.

Aspekter som ny lagstiftning, ökad konkurrens och EU-direktiv gör att Vär-meforsk idag fungerar som en neutral part där medlemmarna kan samarbeta och gemensamt vara med och påverka forskningsprogrammen där resultaten är offent-liga.

Detta examensarbete är en del i ett projekt som genomförs av Qring Technology AB och finansieras av Värmeforsk och handlar om att undersöka om

bers mekaniska egenskaper kan användas för att få mer effektiva sotningsrutiner i värmeverk och värmekraftverk.

1.1.2

Qring Technology AB

Qring Technology AB är ett företag som bildades 2007. Qring har koncentrerat sig på specialistkunskap inom framför allt ljud och vibrationer, se [15]. Företaget består idag av tre specialister med överlappande kompetenser. Företaget arbetar i huvudsak inom energi- och offshorebranscherna.

1.1.3

Borås energi och miljö AB

Signalerna som har analyserats i denna rapport kommer i huvudsak från en av pannorna i Ryaverket i Borås som ägs av Borås energi och miljö AB, se [16]. Rya-verket har funnits sedan 60-talet och producerar fjärrvärme, el och fjärrkyla. Den panna från vilken data har hämtats vid Ryaverket är av typen fluidbäddpanna och används mest till att förbränna avfall. Sotning av överhettartuberna sker normalt en gång om dagen och genomförs med hjälp av ångsotning, se avsnitt 1.2.4.

1.1.4

SAKAB

De första mätningarna i det projekt som är en bakgrund till detta examensarbete gjordes i SAKAB:s panna WTE11_{i Normstorp. Pannan används för att förbränna} kemiska föroreningar som är svåra att bryta ned.

Pannan är en rotationsugn som används för att elda farligt avfall. Avfallet kommer till pannan med emballaget kvar. Detta innebär att explosioner inträffar då, till exempel, en tät burk med färg värms upp. Förbränningen blir momentant mycket kraftig i samband med en sådan explosion. Detta resulterar i att överhettar-tuberna smutsas ned snabbt och behöver därför vanligtvis sotas en gång var 30:e minut. Sotningen som genomförs i WTE1 är av typen slagsotning, se avsnitt 1.2.4.

1.2

Bakgrund

Som bakgrund till problemet beskrivs här vilken typ av värmepanna som analy-serats hos Ryaverket i Borås, en kort förklaring om mekaniska vibrationer, några förklaringar till hur och varför beläggningar uppkommer på överhettartuberna och hur sotning av dessa går till. Slutligen beskrivs idén som ligger till grund för pro-jektet och detta examensarbete.

1.2.1

Fluidbäddpannor

I en fluidbäddpanna blandas bränslet av luftströmmar som blåser underifrån där bränslet ligger vilket gör att den lyfter från bädden. På det viset blandas bränslet och syret vilket ger en effektivare förbränning, se [19]. Dessa pannor har fördelen

1

1.2 Bakgrund 5

att man kan använda många olika typer av bränslen som andra pannor inte klarar av. Den heta gasen från förbränningen passerar en mellandel för att upprätthålla kravet i direktivet EG 2000/76 med 850◦

C i minst två sekunder som går att läsa i [1]. Denna regel finns för att oxideringen av kolföreningar skall vara så fullständig som möjligt. Genom en fullständig förbränning minskar risken för återbildningar av dioxiner och furaner.

Nästa steg i pannan är överhettartuberna där ångan i tuberna överhettas. Överhettartuberna är i praktiken en tvärsströmstubvärmeväxlare, det vill säga gaserna från förbränningen som värmer upp ångan i överhettartuberna passerar tuberna vinkelrätt.

Rökgaserna kyls ned ytterligare i en ekonomiserdel som är en fortsättning av överhettartuberna men där tuberna innehåller vatten istället för ånga. Denna del används för att förvärma, till exempel, matarvattnet till pannan.

Ryaverkets pannor är också försedda med rökgaskondensering som förbättrar verkningsgraden ytterligare.

1.2.2

Mekaniska vibrationer

För att förklara mekaniska vibrationer så betraktas här ett rör som exempel. En överhettartub är format som ett rör, vilket gör det lättare att förstå de mekaniska vibrationerna som ska mätas och orsaken till deras beteende. Mekaniska vibratio-ner, eller svängningar, delas enligt [17] upp i tre karakteristiska typer utifrån hur de utbreder sig i röret:

• Longitudinell svängning • Transversell svängning • Skjuvsvängning

Longitudinell svängning är vibrationer som sker genom att ett rör komprimeras och expanderas längs med rörets riktning. Detta är också bland annat sättet som ljudvågor förflyttar sig i luften när man pratar. Transversell svängning i ett rör är istället vibrationer som rör sig vinkelrätt mot utbredningsriktningen. Svängning-arna hos en gitarrsträng är transversella. Skjuvsvängningar sker längs med ytan på röret vinkelrätt mot rörets riktning. Ett exempel på detta är ringsvängningar, se Figur 1.1.

Longitudinell svängning Transversell svängning Skjuvsvängning

x x x

Figur 1.1.Exempel på långitudinell, transversell och skjuvsvängning. x-axeln visar rö-rets längdriktning.

Ett systems naturliga frekvens fn är den frekvens som systemet svänger med

när externa krafter slutat verka på systemet. Även här är en gitarrsträng ett bra exempel. Genom att slå på strängen utsätts den för en påtvingad svängning fp.

När fingret släpper strängen kommer den att svänga med den frekvens som den är stämd för, fn.

Om ett system, till exempel en gunga, tvingas svänga med en viss frekvens kallas det för en påtvingad svängning fp, se [3]. Sammanfaller den påtvingade

svängningen med den naturliga frekvensen fn kommer svängningen att förstärkas.

Detta kallas resonans. På det viset går det att med liten kraft få en väldigt stor reaktion hos ett system. Om gungan tillåts svänga fram och tillbaka av sig själv sker det med dess naturliga frekvens. Händer inget med gungan kommer den fort-farande gunga med dess naturliga frekvens samtidigt som svängningens amplitud sakta avtar. Att svängningen avtar kallas för dämpning.

Genom att putta på gungan och försöka få den att svänga med en annan fre-kvens fp kommer gungan att göra det så länge som den tvingas svänga med den

frekvensen. Släpps gungan återgår den att svänga med sin naturliga frekvens fn.

Puttas gungan med samma frekvens som den naturliga frekvensen kommer gungan istället att gunga högre och högre. En liten påtvingad kraft förstärks av systemet vilket gör att naturlig frekvens också kallas för resonansfrekvens fr. Således

kom-mer i rapporten resonansfrekvens att beskriva då påtvingad frekvens och naturlig frekvens sammanfaller.

1.2.3

Beläggningar

Beläggningar, eller påslag, på överhettartuberna och andra ytor skapas av stoft och aska som frigörs vid förbränning i pannan. I Figur 1.2 är två foton tagna inuti Ryaverkets panna som visar hur beläggningar på överhettartuberna kan se ut. Typen av bränsle som används påverkar bildandet och egenskaperna hos dessa beläggningar. Förhållandena, såsom temperatur och rökgasflöde, i värmepannan orsakar olika typer av processer som leder till att beläggningar byggs upp. Ett par exempel på dessa processer är impaktion, kondensation och termoforces som beskrivs i [5].

Impaktion leder till att flygande partiklar i pannan fastnar på överhettartu-berna. Kondensation uppstår då gasformiga partiklar kondenserar och stelnar mot den relativt kallare överhettartuben. Ibland kan den kondenserade vätskan leda till att fasta partiklar limmas fast på tuben.

Termoforces åsyftar den kinetiska energin som partiklarna får på grund av temperaturgradienten i pannan. Vid de kallare ytorna i pannan har gasen lägre kinetisk energi än där det är varmare, vilket resulterar i att partiklar fastnar och beläggningar byggs upp på överhettartuberna. Olika stora partiklar i gasen påverkas olika mycket av de olika processerna.

1.2.4

Sotning

I [5] beskrivs olika typer av sotningsmetoder som används för att avlägsna belägg-ningar på (och inuti) överhettartuberna. Sammantaget för de flesta metoderna är

1.2 Bakgrund 7

att de utsätter pannan och överhettartuberna för stora mekaniska påfrestningar. Ett exempel på en sotningsmetod är användning av ånga för att avlägsna belägg-ningarna, så kallad ångsotning. Denna metod kan dock leda till erosionsskador eller sprickbildningar på överhettartuberna. Andra typer av sotningsmetoder använder slag, explosioner eller olika typer av metallkulor för att slå bort beläggningarna. Det kan dock nämnas att metallkulor företrädelsevis används för rengöring inuti överhettartuberna. Även här kan metoderna orsaka oönskade sprickbildningar och annat slitage på överhettartuberna.

När sotning genomförs bör det finnas beläggningar på överhettartuberna för att inte sotning ska ske direkt mot tuberna. Finns det inga beläggningar, till exempel om sotning sker för ofta, kommer sotningen att utsätta överhettartuberna för onormalt slitage vilket leder till att de kan skadas eller går sönder. Är det för långa perioder mellan sotningarna resulterar beläggningarna i ett sämre värmeutbyte. Detta är ineffektivt eftersom värme som produceras i pannan då inte utnyttjas till produktion av ånga utan istället försvinner med rökgasen.

Sotning sker idag oftast med jämna tidsintervall då det antagits att beläggning-ar hbeläggning-ar hunnit bildas på överhettbeläggning-artuberna, eller när operatörerna, som bland annat ansvarar för rengöringen av pannorna, anser att det behövs. Båda metoderna har brister, då ingen tar hänsyn till den aktuella storleken på beläggningarna. Ibland behöver pannan rengöras oftare och ibland mer sällan beroende på hur snabbt beläggningarna bildas. Det händer att sotningen inte avlägsnar allt slagg och sot som sitter på tuberna, och då fortsätter dessa att byggas på under tiden tills nästa sotning sker. Det vore därför intressant om det går att detektera sotpåslag som skulle kunna användas för att anpassa sotningen. Sotningen skulle därmed kunna bli mer behovsstyrd, vilket skulle effektivisera hela rengöringsprocessen och minska risken för slitage.

1.2.5

Detektering av beläggningar på överhettartuber

Den idé som ligger till grund för det projekt som detta examensarbete är en del av, går att läsa om i [4]. Det handlar om att mäta vibrationerna hos överhettartuber-na för att aöverhettartuber-nalysera dess resoöverhettartuber-nansfrekvenser och mäta hur resoöverhettartuber-nansfrekvenseröverhettartuber-na förändras över tiden på grund av beläggningarna som bildas på tuben. På grund av sotet som transporteras i pannan och dess egenskaper bör beläggningarna som bildas på överhettartuberna resultera i att tubernas resonansfrekvenser avtar. I [5] beskrivs hur knackningsprov har genomförts på ett rör i normal rumstempera-tur där resonansfrekvenserna analyserats före och efter att ett lager gips, som ska motsvara verkliga beläggningar, applicerats på röret. Under ett knackprov utsätts en struktur för slag, till exempel med en hammare, för att excitera dess resonans-frekvenser. Resultatet från testet var att det gick att se att resonansfrekvenserna hos röret tydligt minskade när det fanns beläggningar på röret. Genom att ana-lysera frekvensförändringen hos överhettartubernas resonansfrekvenser går det då att bestämma beläggningarnas storlek vilket beskrivs nedan.

Under det fortsatta arbetet, som beskrivs djupare i [5], har olika typer av givare analyserats och hur och var dessa ska monteras för att bäst kunna mäta resonansfrekvenserna på överhettartuberna under drift. De höga temperaturerna

Figur 1.2. Två foton tagna inifrån Ryaverkets panna efter sex månaders drift. Den övre bilden visar ett fotografi av sensorerna som använts för att mäta vibrationerna. Mätutrustningen är skyddad för att klara av påfrestningarna i pannan. Denna bild går att jämföra med Figur 4.1 som visar hur det såg ut när den precis hade monterats sex månader tidigare. Den undre bilden visar överhettartuberna. I båda bilderna syns belägg-ningarna som har bildats utanpå tuberna och som även hänger som små tappar under. Även om pannan sotas med jämna mellanrum byggs det kontinuerligt på beläggningar som inte går bort.

1.3 Mål 9

i värmepannan förkortar livslängden hos utrustningen. Partiklar från bränslet ger frätnings- och blästringsskador, vilket kräver att givare och kablar är ordentligt skyddade för att kunna fungera över tid i en sådan krävande miljö. Knackprov går inte att genomföra under drift, utan istället är det vibrationer som åsamkas av pannan som mäts. Mätutrustning som klarar de extrema förhållandena har framtagits och mätningar visar att de monterade givarna detekterar vibrationer hos tuberna. Därmed bör även information om resonansfrekvenserna kunna finnas med i dessa mätdata.

En detekterad frekvensförändring kan översättas till en massförändring hos överhettartuben utifrån en förenklad relation nedan enligt [4].

ff = C

r k mm

(1.1) ff är resonansfrekvensen, mm är massan, k är styvheten och C är en

proportio-nalitetskonstant. Styvheten antas vara konstant och då kan förhållandet mellan förändringen i frekvens δff och förändringen i massa δmmskrivas som

δmm= C2k 1 (ff + δff)2− 1 ff2 ! (1.2) eller δff = C √ k  ₁ √ mm+ δmm− 1 √_m m  (1.3) Både i (1.2) och (1.3) så kan en avtagande frekvens, δff < 0 ses som en massökning

δmm och tvärtom.

Syftet med examensarbetet är att följa små frekvensskiften hos resonansfre-kvenserna som sedan kan användas för att beräkna massförändring hos överhettar-tuberna. Svårigheten är att detektera och följa förändringarna i frekvens δff.

1.3

Mål

I detta examensarbete ska olika metoder för att skatta frekvensinnehåll utvär-deras på simulerade signaler. Metoderna ska sedan användas för att analysera mätdata från överhettartuber för att detektera förändring hos tubernas resonans-frekvenser. Utmaningen är att detektera, sortera samt värdera förändringen hos resonansfrekvenserna på bästa möjliga sätt. Med hjälp av frekvensförändringarna kan sotningsrutinerna för överhettartuberna förbättras.

De huvudsakliga målen med detta arbete kan sammanfattas i tre punkter: 1. Ett förarbete ska genomföras där olika signalmodelleringsmetoder skall

un-dersökas och analyseras för att se hur resonansfrekvenser kan identifieras samt hur brus och dynamik kan hanteras.

2. Utifrån genomförda mätningar från överhettartuber i värmepannor skall det undersökas om resonansfrekvenser kan identifieras hos dessa mätningar.

3. En algoritm där följning av resonansfrekvenser genomförs i realtid skall, om möjligt, implementeras i Matlab och en robusthetsanalys skall i så fall ge-nomföras. En diskussion ska finnas i rapporten hur styrparametrar för sot-ning kan bestämmas där det går att säga när sotsot-ning ska ske.

1.4

Upplägg

Hänvisningar i texten till ekvationer markeras med parentes, till exempel (1.1) är fortlöpande numrering genom kapitlet, där första siffran visar kapitel och den and-ra beskriver ekvationens position. Referenser till annan litteand-ratur markeand-ras med hakparentes (till exempel [1]) där siffran motsvarar ordningen i litteraturförteck-ningen.

Innehållet i rapporten är uppdelat i följande kapitel:

Kapitel 1 berättar om bakgrunden till examensarbetet: Parter, enkla be-skrivningar om vibrationer, värmepannor och beläggningar. Slutligen definieras examensarbetets mål.

Kapitel 2 förklarar teorin bakom de metoder och algoritmer som används i denna rapport.

Kapitel 3innehåller tester på simulerade data där olika egenskaper och pa-rameterval hos metoderna och algoritmerna i Kapitel 2 utreds.

Kapitel 4går igenom hur signaler har inhämtats och innehåller analyser om signalinnehåll och störningar hos signalerna.

Kapitel 5 innehåller en del av de skattningar som gjorts med de olika me-toderna och algoritmerna på verkliga signaler från Ryaverkets panna. Resultatet från de olika metoderna diskuteras.

Kapitel 6 beskriver en målföljningsalgoritm som används för att hitta reso-nansfrekvenser i skattningarna gjorda i Kapitel 5. Slutligen testas algoritmen på tidigare gjorda skattningar och resultatet diskuteras.

Kapitel 7 sammanfattar resultat och slutsatser. Hur parametrar ska väljas för de olika algoritmerna diskuteras utifrån slutsatserna. Ett förslag på hur sot-ning skulle kunna bli mer behovsstyrd presenteras. Slutligen presenteras förslag på fortsatt arbete.

Appendix Ainnehåller implementerad Matlab-kod för skattning av resonans-frekvenser och målföljningsalgoritmen.

Kapitel 2

Skattningsmetoder

I detta kapitel presenteras teorin bakom de metoder som har använts för att de-tektera resonansfrekvenser i signaler från de givare som monterats på överhettar-tuberna. Mer om hur givarna monterades går att läsa i avsnitt 4.1.

När en signal y samplas avgör samplingstiden T hur mycket av signalinne-hållet som bevaras i den samplade signalen. Höga frekvenser kan misstolkas som låga vilket kallas för aliaseffekten. Den kritiska frekvensen som samplingen måste överstiga bestäms av samplingsteoremet nedan, se [6].

1

T = 2fmax (2.1)

2.1

Fouriertransform och periodogram

Den diskreta fouriertransformen för en samplad signal y(t) är definierad i [6] som Y (ei2πf T_{) = T}

N −1 X

k=0

h(k)y(k)e−i2πf kT _(2.2) där N är antal data, f är frekvensen i Hz, h(k) är en fönsterfunktion, y(k) är signalen och T är samplingstiden. Exempel på olika fönsterfunktioner ges nedan. Det finns flera olika effektiva beräkningsalgoritmer för diskret fouriertransform. En sådan algoritm kallas för FFT1_.

Effektspektraltätheten, eller frekvensinnehållet, för y kan skattas med hjälp av ett periodogram som definieras nedan.

ˆ ΩN(f) = 1 N T  Y (ei2πf T)   2 (2.3) Periodogrammet detekterar dominerande frekvenser som toppar hos ˆΩN. För

att få ut en eller flera frekvenser måste periodogrammet avläsas ”manuellt”.

1

Fast Fourier Transform

2.1.1

Welchs metod

Genom att dela upp ett tidsintervall av en signal y i flera delintervall och medel-värdesbilda periodogrammen för alla delintervallen går det att minska påverkan av brus i skattningen. Nackdelen blir dock en försämrad upplösning då det används ett färre antal sampel per periodogram, se [6]. Intervallen kan vara överlappade. Att medelvärdesbilda flera periodogram kallas för Welchs metod. Om intervallen inte överlappar fås ett specialfall som också kallas Bartletts metod.

2.2

Fönstring

Att välja ett delintervall av en mätsignal y(t) och vikta värdet vid varje tidpunkt med en funktion h(t) kallas för att fönstra signalen. Storleken (antal sampel) och utseendet på fönstret (vikten hos varje sampel) påverkar upplösningen när signa-len fouriertransformeras. Fouriertransformen beskriver endast frekvensinnehållet för det valda tidsintervallet av signalen som analyseras. Information om när nå-got hänt i tidsintervallet går inte att utläsa. Genom att välja olika tidsintervall och beräkna fouriertransformen för vart och ett av dem kan förändringar i fre-kvensinnehållet jämföras. För att veta exakt när något inträffar i en signal skulle därför ett kort tidsintervall för varje fönster vara önskvärt. Innehåller fönstret för få sampel uppstår istället problem med frekvensupplösning som gör det svårare att särskilja närliggande frekvenstoppar. Således måste fönsterstorleken väljas så att en tillfredsställande avvägning mellan tidsupplösning och frekvensupplösning uppnås.

Det enklaste fönstret är ett rektangulärt fönster. Med det viktas alla mätningar i det valda tidsintervallet lika mycket och ger en upplösning i frekvensdomänen på

1

N T Hz. Upplösningen beskriver hur tätt frekvenskomponenter kan vara placerade

och samtidigt vara urskiljbara från varandra. Ett större fönster, det vill säga ett större val av N , gör att tätare frekvenskomponenter kan urskiljas. Fönstring av en signal ger också sidolober i frekvensplanet. Frekvensinnehåll i en signal som sprids ut i frekvensplanet när signalen fönstras kallas för läckage.

Ett annat vanligt fönster är Hammingfönstret. Hammingfönstret ger en dubbelt så bred huvudlob vilket ger halva upplösningen jämfört med ett rektangelfönster med samma fönsterlängd, se Figur 2.1, men ger också lägre sidolober. Hamming-fönstret kan enligt [6] beskrivas i tidsplanet av ekvationen

h(t) = 0.54 − 0.46 cos  _2πt N − 1  (2.4)

2.3

Koherens

Koherens är en funktion som studerar förhållandet mellan två signaler y1 och y2. Koherensfunktionen, γ2

y₁y₂(f), beskrivs nedan och definieras i [7] som

γy2₁y₂(f) =

|Gy₁y₂(f)|2

Gy1y1(f)Gy2y2(f)

, 0 ≤ γ2

2.4 Heuristisk metod för gradientanalys av frekvensinnehåll 13 10 20 30 40 50 60 70 80 90 100 110 120 10−2 100 102 Frekvens Läckage för olika fönster

Rektangelfönster Hammingfönster

Figur 2.1.Rektangelfönster och hammingfönster i frekvensplanet.

där Gxy(f) är korsspektraltäthetsfunktionen och Gxx(f) är

autospektraltäthets-funktionen och definieras som Gy₁y₂ = ∞ X k=−∞ Ry₁y₂(m)e −i2πf m Gyiyi = ∞ X k=−∞ Ryiyi(m)e −i2πf m_{, i = 1, 2 (2.6)}

Rxy och Rxx definieras som korskorrelationsfunktionen respektive

autokorrela-tionsfunktionen enligt (2.7). Ry₁y₂(m) = Ey1,n+my ∗ 2,n  Ryiyi(m) = Eyi,n+my ∗ i,n  (2.7) Koherensfunktionen har för varje frekvens ett värde i intervallet [0, 1] som säger hur väl frekvensinnehållet för y1och y2överensstämmer. Värden som ligger mellan noll och ett kan betyda att relationen inte är linjär eller att det finns brus i signalerna.

2.4

Heuristisk metod för gradientanalys av

fre-kvensinnehåll

Den algoritm som tidigare har använts av Qring för att detektera frekvensskift i vibrationsmätningarna från överhettartuberna beskrivs i [4]. I princip handlar metoden om att identifiera lokala maxima med hjälp av gradientanalys av FFT för olika tidsintervall hos signalen och sedan beräkna den procentuella frekvens-förändringen hos dessa maxima mellan tidsintervallen. Genom att följa en reso-nansfrekvens över tiden går den procentuella förändringen att omvandlas till en procentuell massförändring hos överhettartuben.

Svårigheten med denna metod är att minska påverkan från lokala maxima i frekvensplanet som inte är resonansfrekvenser. Antag att en maxima som följs inte är en resonansfrekvens. Då kommer förändringen för den frekvenstoppen att ge ett felaktig resultat för den faktiska massförändringen. I en brusig signal finns

många fler maxima än resonansfrekvenser. Därför behövs en statistisk metod för gradientanalys där inverkan av enskilda felaktigt tolkade lokala maxima minskas på det beräknade resultatet.

2.4.1

Manuell analysmetod

Qrings förslag, se [4], är att plocka ut ett stort antal lokala maxima och använda deras positioner i frekvensplanet för olika tidsintervall för att beräkna frekvens-skiftet. Principen handlar om att lokala maxima som beror på brus uppkommer slumpvis och är jämnt fördelade i frekvensplanet. Sammantaget antas alla frekvens-förändringarna som beräknas för dessa maxima då ta ut varandra. Den slutliga frekvensförändringen kommer endast bero på de maxima som tillhör resonansfre-kvenserna.

Att metodiskt identifiera lokala maxima och beräkna förändringar mellan olika tidsintervall av signalen har gjorts för hand. Ska en sådan metod kunna användas i drift behöver detta ske automatiskt. Metoden kan beskrivas i följande steg.

1. Beräkna frekvensinnehållet med FFT, periodogram eller någon annan lik-nande metod, för ett begränsat tidsintervall av signalen.

2. Bestäm frekvensintervall och normera frekvenserna i det intervallet med ro-ten av det kvadratiska medelvärdet för motsvarande intervall av tidssignalen. Plocka sedan ut de lokala maxima som överstiger ett tröskelvärde som väljs ovanför brusnivån.

3. Identifiera maxima i ett tidsintervall med motsvarande maxima från ett an-nat tidsintervall.

4. Beräkna den procentuella förflyttningen hos dessa maxima och beräkna me-delvärdet för alla förflyttningar.

2.4.2

Automatiserad statistisk metod

Ett alternativt förslag till en automatiserad metod för att beräkna den procentu-ella förändringen antar istället att majoriteten av de lokala maxima över ett valt tröskelvärde är resonansfrekvenser. Därmed är den vanligaste procentuella föränd-ringen mellan alla kombinationer av dessa maxima också den korrekta procentuella förändringen.

Antag att frekvensinnehållet har beräknats med hjälp av FFT för olika delin-tervall av en mätsignal y. För att hantera att amplituden hos signalen varierar normeras frekvensinnehållet inom ett valt frekvensintervall [fmin, fmax] utifrån

den högsta toppen. Frekvenser för de toppar inom frekvensintervallet som över-stiger ett förbestämt tröskelvärde sparas i z(t) och z(t + 1) där varje t avser ett bestämt delintervall av mätsignalen. Därefter skattas parametern µ i

z(t) = µz(t − 1) (2.8)

2.5 MUSIC-algoritmen 15

Antag också att det finns åtminstone M mätningar i z(t − 1) respektive z(t) och att det behövs N mätningar från varje t för att kunna skatta µ. Algoritmen plockar för varje iteration M mätningar från z(t − 1) respektive z(t), sorterar dessa i storleksordning och anpassar sedan dessa till (2.8) med minstakvadratan-passning. Alla skattningar av µ sorteras in i intervall med längd η där det vanligast förekommande intervallet blir den slutliga skattningen. Genomförs iterationen till-räckligt många gånger antas de flesta skattningarna av µ ligga nära den korrekta procentuella förändringen med en minsta möjliga precision som bestäms av η.

Exempel 2.1: Exempel med automatiserad statistisk metod

Två uppsättningar med frekvenser beskrivs i tabell 2.1. Gemensamma frekvenser är markerade med fet stil och är 10 % lägre i uppsättning två. För att beräkna parametern µ plockas två toppar slumpvis från respektive uppsättning och sorteras i storleksordning. Därefter skattas µ enligt (2.8). Skattade µ sorteras i intervall om 0,005 %. Försöket görs om 20000 gånger och fördelningen av skattade µ visas i Figur 2.2. µ skattas oftast som 0,9. Detta motsvarar en förändring på -10 % vilket stämmer.

Tabell 2.1.Exempel med automatiserad statistisk metod. Markerade frekvenser är 10 % lägre i uppsättning 2. Uppsättning Frekvenser [Hz] 1 6 6,5 7 7,5 8 9 2 5,4 5,85 6,3 6,75 7,2 5 0.6 0.7 0.8 0.9 1 1.1 0 500 1000

Exempel med automatiserad statistisk metod. Histogram över skattade µ

Figur 2.2. Exempel med automatiserad statistisk metod. Histogram över skattningar av µ. Den högsta stapeln är µ = 0, 90 vilket är korrekt.

2.5

MUSIC-algoritmen

En sätt att hitta dominerande frekvenser i en signal är att använda sig av egenvär-dena till kovariansmatrisen för den uppmätta signalen. En metod som gör detta är MUSIC2_{-algoritmen. I [10] antas att mätdata y(t) kan beskrivas av modellen}

¯

y(t) = A¯x(t) + ¯e(t) (2.9)

2

där

¯

y(t) = y(t) y(t − 1) . . . y(t − m + 1)T ¯

x(t) = x(t) x1(t) . . . xn(t) T

¯e(t) = e(t) e(t − 1) . . . e(t − m + 1)T (2.10) I (2.9) är x(t) signalen, e(t) är normalfördelat vitt brus och A är en Vandermode-matris som har egenskaperna

A = a(ω1) . . . a(ωn)

T

a(ω) = 1 e−iω

. . . e−i(m−1)ωT

rank(A) = n om m ≥ n och ωk6= ωn då k 6= n (2.11)

Efter att signalmodellen är bestämd kan kovariansmatrisen definieras för den mät-ta signalen, ¯y(t), som

R = E[¯y(t)¯y∗ (t)] = AP A∗ + σ2I (2.12) P =    α21 0 . .. 0 α2n   

och om m > n, kommer rang(AP A∗

) = n ge att AP A∗

har n strikt positiva egenvärden och resterande egenvärden är lika med noll. Egenvärdena för AP A∗ sorteras i ickestigande ordning och kan skrivas som

λk > σ2, k = 1, . . . , n

λk = σ2, k = n + 1, . . . , m (2.13)

Låt S = [s1, . . . , sn] vara en mängd ortonormerade egenvektorer som hör ihop

med de n första egenvärdena och G = [g1, . . . , gm−n] är den mängd ortonormerade

egenvektorer som motsvarar de resterande egenvärdena. Ortonormerade vektorer är både ortogonala och har längden ett. I algoritmen bestämmer parametern n hur många frekvenskomponenter som skattas.

Det går att visa att G tillhör nollrummet för A∗ ty RG = G    λn+1 0 . .. 0 λm   = AP A ∗ G + σ2G (2.14) där AP A∗

G = 0. AP har full rang vilket ger att A∗

G = 0 (2.15)

Kolumnerna i G tillhör då nollrummet till A∗. Det går också på motsvarande sätt att visa att

RS = S    λ1 0 . .. 0 λn   = AP A ∗ S + σ2S (2.16)

2.6 ESPRIT-algoritmen 17

vilket kan skrivas om som

S = AC (2.17) där C = P A∗ S       λ1 0 . .. 0 λn    −σ 2    −1 (2.18) som också visar att S spänner upp ett underrum till A. Vidare visar S∗

G = 0 tillsammans med (2.15) att S∗

och A∗

har samma nollrum och därmed så spänner S och A upp samma rum då de också har samma rang.

Med sambanden ovan fås till slut ett resultat som visar att lösningarna {ωk}nk=1 till ekvationen

a∗

(ω)GG∗

a(ω) = 0 (2.19)

också blir de skattade frekvenserna. GG∗

är en ortogonal projektion på underrum-met som spänns upp av G vilket också är nollrumunderrum-met till A∗

. Lösningen kan då sägas innehålla de vektorer som är ortogonala mot bruset, det vill säga signalkom-ponenterna.

Algoritmen kan därefter beskrivas i två steg: 1. Beräkna kovariansmatrisen ˆ R = 1 N N X t=m ˜ y(t)˜y∗ (t) (2.20)

och låt ˆS och ˆG motsvara S och G men skapas från egenvektorerna till ˆR. 2. Skatta vinkelfrekvenserna som de n högsta topparna till

1

a∗(ω) ˆ_{G ˆG}∗_a(ω), ω ∈ [−π, π] (2.21)

2.5.1

Parameterval för MUSIC

MUSIC använder N sampel för att skatta n/2 frekvenskomponenter. I [10] förkla-ras att den statistiska noggrannheten hos estimeringarna i MUSIC-algoritmen ökar då parametern m väljs större, samtidigt ökar också komplexiteten i algoritmen som blir mer beräkningskrävande. Parametern m måste väljas större än parametern n. Om m = n + 1 fås ett specialfall av MUSIC som kallas för Pisarenkos algoritm, se [10].

2.6

ESPRIT-algoritmen

ESPRIT3_{är en algoritm som utnyttjar rotationsmatriser för att skatta frekvenser} i en signal. Algoritmen är hämtad ur [10] och signalmodellen (2.9) som används är

3

densamma som för MUSIC-algoritmen. Låt A1= [Im−1, 0]A och A2= [0, Im−1]A

där Im−1är enhetsmatrisen med dimension (m − 1) × (m − 1) och matriserna A1 och A2 har dimension(m − 1) × n.

Sambandet mellan matriserna, A2= A1D, bestäms av D som är en rotations-matris. D =    e−iω_{1 0} . .. 0 e−iωn    (2.22)

S1and S2definieras som

S1=Im−1 0 S

S2=0 Im−1 S (2.23)

Sambandet mellan S och A ges av S = AC där C är en ickesingulär matris som ges av (2.18). Detta kan användas för att visa att

S2= S1C−1DC = S1Φ där Φ = C−1DC (2.24) Det går att visa att Φ kan bestämmas som

Φ = (S∗ 1S1)

−1 S∗

1S2 (2.25)

ESPRIT-algoritmen beräknar, enligt [10], de dominerande frekvenserna {ωk}nk=1 i signalen som de negativa argumenten till egenvärdena för den skattade matrisen ˆ Φ. ˆ Φ = ( ˆS∗ 1Sˆ1) −1ˆ S∗ 1Sˆ2 (2.26)

där ˆS bestäms på motsvarande sätt som för MUSIC-algoritmen från kovarians-matrisen definierad i ekvation (2.20). Enligt [10] motsvarar noggrannheten hos ESPRIT-algoritmen den för MUSIC men ESPRIT är bättre på att separera rötter som motsvarar signaler och rötter som motsvarar brus.

2.6.1

Parameterval för ESPRIT

Parametervalet för ESPRIT-algoritmen är densamma som för MUSIC där N sam-pel används för att skatta n/2 frekvenskomponenter och parametern m påverkar den statistiska noggrannheten i skattningarna. Valet av parametrarna n och m bestäms på motsvarande sätt som för MUSIC.

2.7

AR-modellering

En annan metod för att analysera signalinnehåll enligt [8] är att försöka anpassa en så kallad konfektionsmodell till en signal y(t) = f (u(t), e(t)) och kan skrivas som

y(t) = 1

q−n+ a1_q−(n−1)+ . . . + a

n

2.7 AR-modellering 19

där q är förskjutningsoperatorn, det vill säga q−1

y(t) = y(t − T ) där T är samp-lingsintervallet. Då kan (2.27) skrivas som

any(t) + an−1y(t − T ) + . . . + y(t − nT ) = e(t) (2.28) Denna modell ger upphov till, som man ser i (2.27), en överföringsfunktion med ett antal poler. Är polerna komplexa polpar kan dessa användas för att skatta dominerande frekenskomponenter i signalinnehållet, se [6]. Komplexa polpar kan ses som naturliga frekvenser.

Genom att betrakta AR-modellen på tillståndsform, i detta fall på en styrbar kanonisk form, se [9], kan modellen för systemet skrivas som

x(t + T ) = Ax(t) + Be(t)

y(t) = Cx(t) (2.29)

Egenvärdena till matrisen A i (2.29) beskriver också polerna till överföringsfunk-tionen. De dominerande frekvenskomponenterna i skattningen kan då i [6] plockas ut som vinkeln θ för polerna om deras värden beskrivs som reiθ_.

För att skatta parametrarna aianvänds i [8] en minstakvadratanpassning.

Me-toden hittar parametrar som minimerar en kostnadsfunktion som i detta fall är definierad som kvadraten av prediktionsfelet ǫ. Om parametervektorn definieras som θ = [a1, a2, . . . , an]T kan det skattade värdet på y(t) skrivas som en funktion

av θ, ˆy(t|θ). Prediktionsfelet för skattningen kan beräknas som

ǫ(t, θ) = y(t) − ˆy(t|θ) (2.30)

Detta går sedan att använda över ett antal sampel t = 1, . . . , N för att skapa ett mått för hur bra skattningen är.

VN(θ) = 1 N N X t=1 ǫ2(t, θ) (2.31)

Parametern θ kan nu väljas till det värde som minimerar kostnadsfunktionen VN(θ).

2.7.1

Parameterval för AR

Parametervalen för AR är färre än för MUSIC och ESPRIT. AR-modellering an-vänder N sampel för att skatta upp till n/2 frekvenskomponenter beroende på antalet komplexa polpar i (2.29). Valet av parametern N bestämmer hur många sampel som modellen ska anpassas efter. Med ett tillräckligt stort val av N kan påverkan av störningar minimeras om signalen antas stationär under tidsinterval-let.

Kapitel 3

Utvärdering av

skattningsmetoder på

simulerade data

I detta kapitel analyseras de olika skattningsmetoderna på simulerade data. Ana-lyserna ger förståelse för hur skattningarna ser ut och vilka egenskaper de har. I avsnitt 3.1 testas AR, MUSIC och ESPRIT på olika genererade signaler. I av-snitt 3.2 testas den automatiserade statistiska metoden. Resultatet från de olika metoderna sammanfattas slutligen i avsnitt 3.3.

3.1

Skattning av simulerade data med AR,

MU-SIC och ESPRIT

För att kunna jämföra och utvärdera skattningarna gjorda med AR, MUSIC och ESPRIT har tester genomförts på signaler med känt innehåll. Målet har varit att få en uppfattning om hur bra metoderna kan skatta och följa, till exempel, olika sinuskomponenter i en brusig signal. Testerna gjordes på olika genererade signaler med samplingsfrekvens 51,3 Hz vilket ger en samplingstid på 0,0195 s. 51,3 Hz motsvarar den samplingsfrekvens som blir om mätdata som inhämtats från Ryaverket nedsamplas till var 40:e sampel.

Testerna är uppdelade i olika genererade signaler som testar olika egenskaper hos skattningsmetoderna och hur olika parameterval påverkar resultatet. De flesta av de genererade signalerna kan beskrivas av signalmodellen

y(t) = a1(t) sin (2πf1(t)t + φ1) + a2(t) sin (2πf2(t)t + φ2) +

+ . . . + ak(t) sin (2πfk(t)t + φk) + e(t) , e(t) ∼ N(0, σe2) (3.1)

där k är antalet sinuskomponenter, ai(t) är amplituden, fi(t) är frekvensen, φi är

fasförskjutningen för respektive komponent och e(t) är vitt brus med variansen σe2.

I avsnitt 3.1.4 analyseras en signal som innehåller en dämpad sinuskomponent. Signalen är genererad genom att vitt brus med variansen ett filtrerats med hjälp av överföringsfunktionen nedan.

Y (s) = s + 10

s2_{+ 15s + 25(2π)}2E(s) (3.2) Samplen som används vid varje skattning med AR, MUSIC och ESPRIT ligger i icke-överlappande intervall. I huvudsak används 10000 sampel per skattning i testerna vilket motsvarar ett tidsfönster på 195 s. Detta betyder att en ny skattning kommer att göras en gång var 195 sekund.

Brusnivån i de genererade signalerna har en varians σ2

e = 50 då (3.1) används

som signalmodell. För att få en uppfattning om de genererade signalerna finns i varje avsnitt nedan en figur som visar en sekund av tidssignalen och ett periodo-gram av de första 10000 samplen.

Avsikten med de olika testerna sammanfattas här:

1. Skatta sinuskomponenter med varierande frekvens - Hur klarar AR, MUSIC och ESPRIT att skatta sinuskomponenter med varierande frekvens? Hur påverkar olika parameterval skattningen?

2. Skatta sinuskomponenter med olika amplitud - Hur skattar AR, MU-SIC och ESPRIT sinuskomponenter med olika amplitud? Hur påverkar pa-rametern n i algoritmerna resultatet hos skattningarna?

3. Skatta sinuskomponenter med varierande amplitud och frekvens -Hur ser skattningarna ut om sinuskomponenterna varierar både i amplitud och frekvens?

4. Skatta en dämpad sinuskomponent - Testerna ovan innehåller ideala sinuskomponenter som skattas. Hur påverkas skattningen om sinuskompo-nenten är dämpad?

3.1.1

Skattning av sinuskomponenter med varierande

fre-kvens

Metoderna som presenterats i Kapitel 2 ska användas för att skatta frekvenskom-ponenter som förflyttar sig i frekvensplanet. Den genererade signalen innehåller två sinuskomponenter och kan beskrivas av (3.1) där parametrarna är definierade i ta-bell 3.1. Amplituden och fasförskjutningen är konstant och frekvensen förändras över tiden.

Figur 3.1 visar ett tidsintervall och periodogram av den genererade signalen. Signalens utseende domineras av bruset.

Parameterinställningarna för AR, MUSIC och ESPRIT är valda så att två frekvenskomponenter skattas. För MUSIC och ESPRIT är parametern m vald till 50.

Figur 3.2 visar resultatet från skattningarna. MUSIC och ESPRIT följer de korrekta frekvenserna men skattningarna hoppar mellan frekvenskomponenterna.

3.1 Skattning av simulerade data med AR, MUSIC och ESPRIT 23

Tabell 3.1.Data för de sinuskomponenter i en genererad signal i avsnitt 3.1.1. Figur 3.1 visar ett tidsintervall av signalen. Sinuskomponenterna har konstant amplitud och fas men varierande frekvens.

Parameter Sinus 1 Sinus 2

ai(t) 1 1 fi(t) 10-0,0001t 12+0,00005t φi 0 0,4 0 1 2 3 4 5 6 7 8 9 −20 −10 0 10 20 Tid [s] Amplitud

Genererad signal med två frekvensskiftande sinuskomponenter i vitt brus ~N(0,50)

0 5 10 15 20 25 10−2 100 102 Frekvens [Hz] Effekt/Frekvens

Genererad signal med två sinuskomponenter med varierande frekvens i vitt brus ~N(0,50)

Figur 3.1.Tidsintervall och periodogram av genererad signal med två sinuskomponenter med varierande frekvens i vitt brus ∼ N(0, 50). Signalen beskrivs av (3.1) då k = 2 och val av parametrar enligt tabell 3.1.

AR-modelleringen hittar frekvenskomponenterna men kan inte separera dem. Där-för är en av skattningarna ofta fel. Skattningarna gjorda med AR varierar mer i frekvens med samma antal sampel per skattning. Vid test med en brusfri signal skattar AR, MUSIC, och ESPRIT frekvenskomponenterna korrekt.

Ökas antalet sampel per skattning till 20000 skattar MUSIC och ESPRIT de korrekta frekvenerna bättre, se Figur 3.3. Således påverkar brus i signalen skattningarna mindre om fler sampel används per skattning. AR kan fortfarande inte separera båda sinuskomponenterna.

Parametern m för ESPRIT och MUSIC väljs genom att prestanda vägs mot beräkningstid. I Figurerna 3.4 och 3.5 visas hur valet av parametern m påverkar resultatet hos skattningen. Om parametern m väljs liten varierar skattningen mer än om m väljs större. Detta gäller både för MUSIC och ESPRIT. Således ger ett stort m en skattning som bättre följer den korrekta frekvensen med mindre variationer.

2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med MUSIC(4,20)−algoritmen.

Korrekt frekvens Skattning 1 Skattning 2 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med ESPRIT(4,20)−algoritmen.

Korrekt frekvens Skattning 1 Skattning 2 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med AR(4)−modellering.

Korrekt frekvens Skattning 1 Skattning 2

Figur 3.2.Skattningar av två sinuskomponenter med varierande frekvens i vitt brus ∼ N (0, 50). MUSIC och ESPRIT följer de korrekta frekvenserna bättre än AR. Skatt-ningarna med de olika algoritmerna hoppar mellan de korrekta frekvenserna.

3.1.2

Skattning av sinuskomponenter med olika amplitud

För att testa hur parametern n påverkar resultatet om fler eller färre frekvenskom-ponenter skattas än vad som finns i signalen jämförs algoritmerna på en genererad signal som består av tre sinuskomponenter med olika amplitud. Amplitud, fre-kvens och fas för de tre sinuskomponeterna finns i tabell 3.2 och Figur 3.6 visar ett tidsintervall och periodogram av signalen.

Tabell 3.2.Data för de sinuskomponenter i en genererad signal i avsnitt 3.1.2. Figur 3.6 visar ett tidsintervall av signalen. Sinuskomponenterna har olika amplitud, frekvens och fas för att se hur skattningarna med AR, MUSIC och ESPRIT påverkas.

Parameter Sinus 1 Sinus 2 Sinus 3

ai(t) 1 1,2 1,4

fi(t) 6 8 10

φi 0 0,1 0,2

Olika värden på parametern n testas i olika skattningar. Resultatet jämförs då två, tre respektive fyra frekvenskomponenter skattas, det vill säga parametern n = 4, 6, 8. Figur 3.7 visar skattningar gjorda med AR. Algoritmen klarar inte

3.1 Skattning av simulerade data med AR, MUSIC och ESPRIT 25 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med MUSIC(4,20)−algoritmen.

Korrekt frekvens Skattning 1 Skattning 2 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med ESPRIT(4,20)−algoritmen.

Korrekt frekvens Skattning 1 Skattning 2 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med AR(4)−modellering.

Korrekt frekvens Skattning 1 Skattning 2

Figur 3.3. Skattningar av två sinus med varierande frekvens i vitt brus ∼ N(0, 50). Skattningarna använder 20000 sampel per skattning vilket är dubbelt så många i jämfö-relse med resultatet i Figur 3.2. Skattningarna blir bättre med fler sampel per skattning. AR skattar de korrekta frekvenserna sämre i jämförelse med MUSIC och ESPRIT med samma antal sampel.

att separarera de tre sinuskomponenterna och när parametern n ökas så sprids skattningarna ut i frekvensplanet.

Figur 3.8 och Figur 3.9 visar skattningar gjorda med MUSIC och ESPRIT. När n = 4 sorteras skattningarna efter storleksordning hos sinuskomponenternas amplituder. Den första skattningen följer sinuskomponenten vid 10 Hz och den andra skattningen följer 8 Hz. Ökas n så att tre frekvenskomponenter skattas, vilket är korrekt antal sinuskomponenter, hoppar skattningarna 2 och 3 mellan sinuskomponenterna 6 Hz och 8 Hz. Om fler frekvenskomponenter skattas än vad som finns i signalen skattas en frekvenskomponent som beror på bruset. Således bör n väljas så att antalet skattningar överensstämmer med det verkliga antalet frekvenskomponenter.

3.1.3

Skattning av sinuskomponenter med varierande

amp-litud och frekvens

För att undersöka hur AR, MUSIC och ESPRIT hanterar sinuskomponenter med varierande amplitud och frekvens jämförs algoritmerna på en genererad signal som innehåller tre sinuskomponenter. Amplituden för varje sinuskomponent ändras

2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med MUSIC(4,5)−algoritmen.

Korrekt frekvens Skattning 1 Skattning 2 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med MUSIC(4,50)−algoritmen.

Korrekt frekvens Skattning 1 Skattning 2

Figur 3.4.Jämförelse av skattningar med MUSIC med två val av parametern m. I den övre bilden är m = 5 och i den undre m = 50. Ett val av en större parameter m ger en skattning med mindre variationer.

2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med ESPRIT(4,5)−algoritmen.

Korrekt frekvens Skattning 1 Skattning 2 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med ESPRIT(4,50)−algoritmen.

Korrekt frekvens Skattning 1 Skattning 2

Figur 3.5.Jämförelse av skattningar med ESPRIT med två val av parametern m. I den övre bilden är m = 5 och i den undre m = 50. Ett val av en större parameter m ger en skattning med mindre variationer.

med jämna tidsintervall.

Den genererade signalen visas i Figur 3.10 där frekvens och fas är definierade i tabell 3.3 och ampliderna i tabell 3.4 enligt (3.1).

Figur 3.11 visar skattningar gjorda med MUSIC, ESPRIT och AR. Algoritmer-na skattar i huvudsak de frekvenser där sinuskomponenterAlgoritmer-na har störst amplitud. Skulle två sinuskomponenter ligga nära varandra i frekvens tolkas de som en sinus. AR har, som tidigare tester visat, svårt att särskilja två närliggande frekvenskom-ponenter.

3.1 Skattning av simulerade data med AR, MUSIC och ESPRIT 27 0 1 2 3 4 5 6 7 8 9 −20 −10 0 10 20 Tid [s] Effekt/Frekvens

Genererad signal med tre sinuskomponenter med olika amplituder i vitt brus ~N(0,50)

0 5 10 15 20 25

100

Frekvens [Hz]

Effekt/Frekvens

Periodogram av genererad signal med tre sinuskomponenter med olika amplituder i vitt brus ~N(0,50)

Figur 3.6.Tidsintervall och periodogram av genererad signal med tre sinuskomponenter som olika amplitud i vitt brus ∼ N(0, 50). Signalen beskrivs av (3.1) då k = 3 och val av parametrar enligt tabell 3.2.

2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med AR(4)−modellering.

Korrekt frekvens Skattning 1 Skattning 2 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med AR(6)−modellering.

Korrekt frekvens Skattning 1 Skattning 2 Skattning 3 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med AR(8)−modellering.

Korrekt frekvens Skattning 1 Skattning 2 Skattning 3 Skattning 4

Figur 3.7.Skattningar med AR(n) mellan olika val av n. AR klarar inte av att skatta sinuskomponenterna vid 6, 8 och 10 Hz.

3.1.4

Skattning av en dämpad sinus

De tidigare simuleringarna innehöll rena sinuskomponenter för att undersöka egen-skaperna hos MUSIC, ESPRIT och AR. För att testa algoritmerna på en känd

2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med MUSIC(4,50)−algoritmen.

Korrekt frekvens Skattning 1 Skattning 2 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med MUSIC(6,50)−algoritmen.

Korrekt frekvens Skattning 1 Skattning 2 Skattning 3 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med MUSIC(8,50)−algoritmen.

Korrekt frekvens Skattning 1 Skattning 2 Skattning 3 skattning 4

Figur 3.8.Jämförelse av skattningar med MUSIC(n,50) mellan olika val av n. Vid för lågt val av n skattar algoritmen de n/2 starkaste frekvenserna. Väljs parametern n så att precis så många frekvenskomponenter som finns i signalen så skattas alla frekvens-komponenter korrekt. Väljs parametern n så att fyra frekvensfrekvens-komponenter skattas är de tre första skattningarna korrekta men den fjärde skattningen varierar.

Tabell 3.3.Frekvens och fas för de sinuskomponenter i en genererad signal med varieran-de amplitud och frekvens. Figur 3.10 visar ett tidsintervall av signalen. Signalen används för att visa hur skattningarna med AR, MUSIC och ESPRIT påverkas när amplituderna hos sinuskomponenterna varierar över tiden. Amplituden för sinuskomponenterna ändras med jämna mellanrum och dessa beskrivs i tabell 3.4

Parameter Sinus 1 Sinus 2 Sinus 3 fi(t) 10-0,0001t 12-0,00005t 14-0,0003t

φi 0 0,4 0,4

men mer svårestimerad signal betraktas en signal som innehåller en dämpad fre-kvenstopp vid 5 Hz, se Figur 3.12. Signalen har genererats enligt (3.2). I de verkliga signalerna är signalinnehållet betydligt mer komplicerat. Den genererade signalen ska visa hur resultatet från algoritmerna påverkas av att frekvenstopparna är mer dämpade. För att skatta frekvenstoppen i Figur 3.12 krävs fler sampel per skatt-ning än i de tidigare fallen för att minska variationerna i skattskatt-ningen.

Figur 3.13 visar resultatet från skattningarna gjorda med MUSIC(2,100) och ESPRIT(2,100) där 1000, 10000 och 40000 sampel använts vid varje skattning.

3.1 Skattning av simulerade data med AR, MUSIC och ESPRIT 29 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med ESPRIT(4,50)−algoritmen.

Korrekt frekvens Skattning 1 Skattning 2 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med ESPRIT(6,50)−algoritmen.

Korrekt frekvens Skattning 1 Skattning 2 Skattning 3 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 5 10 15 20 25 T [s] Hz

Estimering med ESPRIT(8,50)−algoritmen.

Korrekt frekvens Skattning 1 Skattning 2 Skattning 3 skattning 4

Figur 3.9.Jämförelse av skattningar med ESPRIT(n,50) mellan olika val av n. Vid för lågt val av n skattar algoritmen de n/2 starkaste frekvenserna. Väljs parametern n så att precis så många frekvenskomponenter som finns i signalen så skattas alla frekvens-komponenter korrekt. Väljs parametern n så att fyra frekvensfrekvens-komponenter skattas är de tre sista skattningarna korrekta men den första skattningen varierar i frekvensplanet och indikerar felaktigen att det finns en till frekvenskomponent i signalen.

Tabell 3.4.Amplituderna för sinuskomponenterna som skattas i Figur 3.11 med MU-SIC(4,50), ESPRIT(4,50) och AR(4). I intervall 1 är sinus 2 och 3 störst, i intervall 2 är 1 och 2 störst och så vidare.

ai(t) i intervall i Sinus 1 Sinus 2 Sinus 3

1 0,3 0,8 1,0

2 1,0 1,0 0,6

3 0,6 0,8 1,0

4 1,0 0,8 1,0

5 1,0 0,6 0,8

Således om frekvenstoppen är mer dämpad behövs det fler sampel för att få en mindre varierande skattning. Motsvarande tester med AR visar en avvikelse i storleksordningen 1 Hz och skattar istället frekvenstoppen vid 6 Hz.