Der Rechenfehler der BGR:
Nachrechnen mit Originaldaten

Seit März 2020 bin ich mit der BGR in einer Diskussion, weil die Schalldruckangaben der BGR weit oberhalb der Werte anderer Institutionen liegen. Auf einen Twitter-Tweet mit großer Resonanz wurde ich auf die Möglichkeit einer Frag-den-Staat-Anfrage hingewiesen. Die Anfage vom 21.12.2020 zur Offenlegung der Rohdaten und Pegelberechnungen wurde von der BGR am Freitag 22.01.2021 nachmittags kurz vor Ablauf der Auskunftfrist beantwortet. Die Herausgabe der Berechnung wurde abgelehnt. Ich bekam jedoch den Hinweis, dass die Rohdaten frei verfügbar wären und auch eine Anleitung für den Zugriff.

Damit gibt es jetzt für jedermann die Möglichkeit, die Angaben der BGR selbst nachzurechnen. Ärgerlich ist, dass die BGR auf diese Möglichkeit erst nach meiner offiziellen Anfrage hinweist. Anscheinend war das Interesse der BGR bezüglich einer Überprüfung ihrer Berechnungen nicht all zu groß.

Die Ausführungen dieser Seite beziehen sich auf die deutschsprachigen Publikation der BGR "Der unhörbare Lärm von Windkraftanlagen – Infraschallmessungen an einem Windrad nördlich von Hannover". Insbesondere werde ich die Abbildungen 2-4 einer Überprüfung unterziehen. Die exakten Berechnungen finden Sie unten als R-Skript. Im R-Skript ist auch eine Methodenverifikation mit einem künstlich erzeugten Signal (80dB re 20µPa) in einem Rauschhintergrund. Die Verifikation zeigt, dass die Amplitude des Signals (bis auf numerische Ungenauigkeiten) von allen eingesetzten Methoden korrekt detektiert wird.

Zugriff auf Rohdaten der BGR

Für die Auswertung wurden die Daten vom 08.07.2004 00:00 bis 11.07.2004 00:00 heruntergeladen. Der Download kann über die folgende URL erfolgen:

http://eida.bgr.de/fdsnws/dataselect/1/query?station=HUF01&channel=HDF&starttime=2004-07-08T00:00:00&endtime=2004-07-11T00:00:00

Die Daten sind im mseed-Format. Zur Verarbeitung benötigt man ein Paket, das in der Lage ist, dieses Dateiformat einzulesen. Alternativ kann man die Daten mit dem mseed2ascii-Tool in eine Textdatei umwandeln. Eine Online-Variante bis 100MB finden Sie unter dieser URL: https://bayconf.bayceer.uni-bayreuth.de/mseed2ascii.php

Nach dem Einlesen müssen die Daten skaliert werden. Der Skalierungsfaktor ist in den Metadaten angegeben.

Überprüfung des Skalierungsfaktors

In Abbildung 2 der BGR ist ein Ausschnitt aus dem Hochpass gefilterten Drucksignal geplottet. Der Abschnitt ist zeitlich nicht genau datiert. Er ist jedoch aus dem Bereich mit 26rpm.

BGR Abbildung 2: Drucksignal

Die folgende Abbildung zeigt einen 15-Sekundenausschnitt aus dem Drucksignal ab 2004-07-09 22:28:40. Zu dieser Zeit rotierte das Windrad mit 26rpm.

Drucksignale mit unterschiedlicher Filterung berechnet mit den unten angehängten R-Skript und den Rohdaten der BGR

Auch wenn die BGR von einem Hochpass gefilterten Signal spricht, kommt das Signal mit Bandpassfilter der BGR-Abbildung am nächsten. Grundsätzlich ist festzustellen, dass die Amplituden vergleichbar sind. Damit sollte das Einlesen der Rohdaten und die Umwandlung in ein Drucksignal korrekt sein.

Nachtrag 25.02.2022:
Inzwischen hat die BGR eingeräumt, dass das in Abbildung 2 gezeigte Drucksignal fälschlicher Weise mit "Standort 3 in etwa 200m Entfernung zum Windrad" bezeichnet war. Tatsächlich ist das Drucksignal vom Standort 2 in ca. 100m Abstand zum Windrad.

Überprüfung der spektralen Daten

In Abbildung 4 zeigt die BGR ein Spektrogramm über drei Tage. Dies sind exakt die drei Tage, für die Rohdaten heruntergeladen wurden.

BGR Abbildung 4: Spektrogramm

Führt man auf Basis der im ersten Teil des R-Skriptes beschriebenen Methode eine Spektralanalyse durch und normiert das Ergebnis wie für Schalldruckpegel (SPL) üblich auf 20µPa, erhält man das folgende Spektrogramm. Das Spektrogramm wurde mit dem specgram Befehl aus dem oce-Package berechnet. Zur Anwendung kam ein Hanningfenster mit Länge 8192. Die Bin-Breite beträgt 0,012Hz.

Spektrogramm berechnet mit den unten angehängten R-Skript und den Rohdaten der BGR

Grundsätzlich werden die wesentlichen Muster des BGR-Sprektrogramms reproduziert. Auffällig sind jedoch die deutlich niedrigeren SPL-Werte. Während in der BGR-Grafik die Werte bis über 90dB reichen, endet die Skala hier bei 55dB. Auch fehlen im BGR-Spektrogramm manche verrauschte Bereiche (z.B. ab 09.07. 00:00) vollständig. Dies ist nicht nachvollziehbar. Interessant ist, dass es neben den zwei Windradharmonischen (20 und 26rpm) noch eine weitere harmonische Frequenz gibt, die nicht diskutiert wurde (vgl. Schmalbandige, frequenzkonstante Signale).

Abbildung 3 der BGR zeigt ein Frequenzspektrum. Die genauen Zeitslots für die Spektren wurden nicht angegeben.

BGR Abbildung 3: Frequenzspektrum

In der folgenden Abbildung wurde versucht die BGR Abbildung 3 zu reproduzieren. Die Zeitfenster wurden anhand des obigen Spektrogramms ausgewählt. Die 20 rpm +-Reihe ist ein Zeitpunkt, an dem neben der Windradharmonischen ein weiteres harmonische Signal detektiert wurde.

Frequenzspektren berechnet mit den unten angehängten R-Skript und den Rohdaten der BGR

Wenig überraschend sind auch in diesem Diagramm die Schalldruckpegel deutlich niedriger. Die 2.BPH bei 26 rpm wird im Skript mit unterschiedlichen Berechnungsmethoden konsistent mit 59dB re 20µPa (inkl. Hintergrund) quantifiziert. Die BGR gibt für das gleiche Signal basierend auf den gleichen Daten 90dB re 20µPa an. Damit ergibt sich ein Rechenfehler >30dB.

Fazit

Die Auswertung der Rohdaten bestätigt das Bild des Rechenfehlers. Die Indizien sind in meinen Augen erdrückend:

1. Alle anderen Institutionen messen deutlich niedrigere Schalldruckpegel an WEA. Das ist physikalisch nicht erklärbar. Zu diesem Schluss kommt sogar Herr Hübner von Windwahn.com (vgl. Diskussionsseite BGR-Studie) und im Grunde auch die BGR (vgl. Antwort der BGR auf Fragenkatalog über BMWi-Bürgerdialog)
2. Das von der BGR publizierte Drucksignal passt nicht zu den publizierten Schalldruckpegeln (vgl. Warum die Schalldrücke der BGR falsch sind)
3. Rechnet man mit den Rohdaten der BGR kann man das Drucksignal jedoch nicht die Schalldruckpegel reproduzieren. Die Diskrepanz zwischen den publizierten Werten der BGR und den hier errechneten Werten beträgt über 30dB

Neben dem Rechenfehler geht in die Hochrechnung der BGR ein Normierungsfehler mit ein, da das vermessene Windrad (Vestas V47) tatsächlich 660kW Leistung hatte und nicht wie von der BGR angegeben nur 200kW. Dieser Normierungsfehler sorgt für weitere 5dB. Der Modellfehler beträgt somit über 35dB. Das ist mehr als knapp daneben. Als wissenschaftlich völlig inakzeptabel erachte ich, dass die BGR nie versucht hat, ihr Modell mit inzwischen vorliegenden Messdaten größerer Windräder zu validieren.

Ich würde nicht soviel Zeit und Mühen in dieses Thema investieren, wenn ich nicht erlebt hätte, wie mit dem Infraschall-Argument ein Windenergieprojekt in der Nähe meines Wohnorts verhindert wurde. In der Argumentation der Windkraftgegner spielte der ZDF-Film und die hohen Schalldruckpegel der BGR eine Schlüsselrolle.

Ich hoffe, dass die BGR mit Veröffentlichung dieser Seite endlich zur Erkenntnis gelangt, dass ihre Auswertungen nicht den anerkannten Methoden der Schalldruckberechnung entsprechen. Bei Eingeständnis des Fehlers wäre es in meinen Augen angebracht, sich öffentlich bei der LUBW für die irreführenden bis falschen Aussagen im ZDF-Film zu entschuldigen. Selbstverständlich bin auch ich bereit, mich bei der BGR zu entschuldigen, sollten die Ausführungen auf dieser Seite nachweislich falsch sein.

Es ist kein Problem, einen Fehler zu machen. Ein Problem wird es jedoch, wenn man nicht bereit ist, diesen trotz erdrückender Beweislast zuzugeben. In den letzten Jahren konnten wir in den USA sehr eindrücklich erleben, was es bedeutet, wenn staatliche Institutionen sich weigern, wissenschaftliche Fakten anzuerkennen. Gegen derartige Zustände sollten wir im Interesse unserer Demokratie kämpfen.

Anhang: R-Skript zur Berechnung

## Libraries (R-Packages)
library(oce)
library(signal)

###########################################################
## 1. General test if SPL-Calculations give correct results
###########################################################
## time vector (1000s)
tt=1:1e6/100
## signal at 5.021Hz with 0.2Pa == 80dB re 20µPa
s=sqrt(2)*0.2*sin(5.021*2*pi*tt)+0.1*rnorm(length(tt))
## 
nfft=8192
window=hanning(nfft)
spec = specgram(x = s, n = nfft, Fs = 100,
                window = window, overlap = nfft/2 )

## https://www.sjsu.edu/people/burford.furman/docs/me120/FFT_tutorial_NI.pdf
# convert to PA (rms), normalize to 20µPa and convert to dB
# 1.5 is the noise power bandwidth of hanning window
SPL = 20*log10(sqrt(2)*abs(spec$S)/sum(window)/sqrt(1.5)/20e-6)
plot(spec$f,SPL[,1],type='l',ylim=c(35,81),xlim=c(4.5,5.5),xlab='[Hz]',ylab='SPL [dB] re 20µPa')

## Sum up some adjacent bins +/- 0.2Hz (single FFT)
sel=spec$f>4.821 & spec$f<5.221
dbSum<-function(...){
  10*log10(sum(sapply(c(...),function(db) 10^(db/10))))
}
dbSum(SPL[sel,1]) ## Should give approximately 80dB

## Sum up some adjacent bins +/- 0.2Hz (rms averaged FFT)
SPL = 20*log10(sqrt(2)*sqrt(rowMeans(abs(spec$S)^2))/sum(window)/sqrt(1.5)/20e-6)
plot(spec$f,SPL,type='l',ylim=c(35,81),xlim=c(4.5,5.5),xlab='[Hz]',ylab='SPL [dB] re 20µPa')
dbSum(SPL[sel]) ## Should give approximately 80dB

psd_rect=pwelch(ts(s/20e-6,frequency = 100),window = rep(1,nfft),plot=F,detrend=F)
plot(psd_rect$freq,10*log10(psd_rect$spec),type='l',ylim=c(52,100),xlim=c(4.5,5.5),xlab='[Hz]',ylab='PSD [dB] re (20µPa)²/Hz')
psd_hann=pwelch(ts(s/20e-6,frequency = 100),window = hanning(nfft),plot=F,detrend=F)
plot(psd_hann$freq,10*log10(psd_hann$spec),type='l',ylim=c(52,100),xlim=c(4.5,5.5),xlab='[Hz]',ylab='PSD [dB] re (20µPa)²/Hz')
psd_hamm=pwelch(ts(s/20e-6,frequency = 100),window = hamming(nfft),plot=F,detrend=F)
plot(psd_hamm$freq,10*log10(psd_hamm$spec),type='l',ylim=c(52,100),xlim=c(4.5,5.5),xlab='[Hz]',ylab='PSD [dB] re (20µPa)²/Hz')
## Compared to SPL curves PSD is shifted to higher dB values 
## The difference is 10*log10(2*sample_frequency/fft_length) = -16.1
10*log10(2*100/nfft)

## Integrating the pwelch results using the bin width (100/nfft).
## Apparently we need a factor two. pwelch probably only gives the
## one sided information
10*log10(2*sum(psd_rect$spec[sel])*100/nfft)
10*log10(2*sum(psd_hann$spec[sel])*100/nfft)
10*log10(2*sum(psd_hamm$spec[sel])*100/nfft)
## All pwelch calculations yield approximately 80dB
#######################################################
## RESULT: Calculation is correct.
## both methods (specgram, pwelch) are able to quantify 
## the 80dB signal correctly
## 

#######################################################
### 2. Analysis of BGR-Data  
##  and compare to BGR results in
## https://www.bgr.bund.de/DE/Themen/Erdbeben-Gefaehrdungsanalysen/Seismologie/Downloads/infraschall_WKA.pdf;jsessionid=342295BDBE13463536E976ED361E12CF.2_cid284?__blob=publicationFile&v=2
#######################################################
### Get the BGR raw-data (available to the public) with URL:
## http://eida.bgr.de/fdsnws/dataselect/1/query?station=HUF01&channel=HDF&starttime=2004-07-08T00:00:00&endtime=2004-07-11T00:00:00
## 
## For metadata use
## http://eida.bgr.de/fdsnws/station/1/query?station=HUF01&channel=HDF&starttime=2004-07-08T00:00:00&endtime=2004-07-11T00:00:00&level=response
##  
##
## Transform into ASCII-data
## https://github.com/iris-edu/mseed2ascii/blob/master/doc/mseed2ascii.md
## or online (max. 100MB) 
## https://bayconf.bayceer.uni-bayreuth.de/mseed2ascii.php
## 
setwd("~/Dokumente/Projekte/Infraschall/BGR/Rohdaten/")
d=read.csv("GR.HUF01..HDF.D.2004-07-07T235958.432000.txt",skip = 1, header = FALSE)
##
## Transform INT32-Data to real pressure values (Pa) (-> metadata)
pa=d$V1/10487.68

######################################################################
## Check if amplitude of Signal is comparable to BGR (Abbildung 2)
######################################################################
bp=butter(4,c(0.5,10)/50,'pass') ## Band-Passfilter 0.5-10Hz
hp=butter(4,0.5/50,'high') ## Highpass 0.5Hz
pa_bp <- filter(bp, pa) # apply filter
pa_hp <- filter(hp, pa) # apply filter
## Select time section with 26rpm (see below)
## The exact time stamp is not given in the BGR paper
s=46.7*3600*100+1:1500-80000
as.POSIXct('2004-07-08 00:00')+s[1]/100
png("raw_signal.png",width = 7, height = 10, res = 300, units = 'in')
par(mfrow=c(3,1))
plot(1:1500/100,pa[s],type='l',xlab='Zeit [s]',ylab='[Pa]',main='Raw signal')
plot(1:1500/100,pa_hp[s],type='l',xlab='Zeit [s]',ylab='[Pa]',main='Signal high pass filtered 0.5Hz')
plot(1:1500/100,pa_bp[s],type='l',xlab='Zeit [s]',ylab='[Pa]',main='Signal band pass filtered 0.5-10Hz')
dev.off()
######################################################################
## RESULT: Amplitude of filtered signal is ins the same range +/- 0.1Pa as BGR
## Scaling of raw data seems ok but
## - BGR plot does not seem to include high frequency data
##   looks like a band pass filtered signal
## - I did not find a sample with dominant blade passing signal


######################################################################
## Check if spectrogram (Abbildung 4) can be reproduced
######################################################################
# number of points to use for the fft
nfft=8192
# window 
window=hanning(nfft)
# overlap (in points)
overlap=nfft/2
# frequency
fs=100
# create spectrogram
spec = specgram(x = pa, n = nfft, Fs = fs,
                window = window, overlap = overlap )

# scale to Pa (rms), normalize to 20µPa and convert to dB (page 5+6)
# 1.5 is the noise power bandwidth of hanning window
SPL = 20*log10(sqrt(2)*abs(spec$S)/sum(hanning(nfft)/sqrt(1.5))/20e-6)

# plot spectrogram
png("spectrogram.png",width = 10, height = 7, res = 300, units = 'in')
imagep(x = as.POSIXct('2004-07-08 00:00')+spec$t,
       y = spec$f,
       z = t(SPL),
       col = oce.colorsJet,
       ylab = 'Frequency [Hz]',
       xlab = 'Time [s]',
       zlab = paste('SPL [dB] re 20µPa, bin =',round(100/nfft,3),'Hz'),
       drawPalette = T,
       decimate = F,
       zlim = c(30,55),
       ylim = c(0,10)
  )
dev.off()
######################################################################
## RESULT: main patterns can be reproduced but
## - SPLs of BPHs are below 60dB (BGR 90dB)
## - Compared to BGR plot, some regions in the BGR-Plot
##   seem amplified or deadened. Reason is unclear


######################################################################
## Check if spectra (Abbildung 3) can be reproduced
######################################################################

labels=c('26 rpm','20 rpm', '0 rpm','20 rpm +')
times=list(c(60,63),c(12,15),c(22,25),c(52,55))

SPLs=list()
psd_rect=list()
psd_hann=list()
nfft=8192
for(i in 1:length(labels)){
  pa_s=pa[(times[[i]][1]*3600*100):(times[[i]][2]*3600*100)]
  spec = specgram(x = pa_s, n = nfft, Fs = 100,
                  window = hanning(nfft), overlap = nfft/2 )
  P=sqrt(rowMeans(abs(spec$S)^2))
  SPL = 20*log10(sqrt(2)*P/sum(hanning(nfft))/sqrt(1.5)/20e-6)
  SPLs[[i]]=SPL
  
## In additions to direct FFT-Calculations we use
## pwelch with hanning and rectangular window
## Note: pwelch gives the PSD, detrending (Default) is disabled
  psd_rect[[i]]=pwelch(ts(pa_s/20e-6,frequency = 100),window = rep(1,nfft),plot=F,detrend=F)
  psd_hann[[i]]=pwelch(ts(pa_s/20e-6,frequency = 100),window = hanning(nfft),plot=F,detrend=F)
  
}

## Test for overlap of different methods 
plot(spec$f,SPLs[[1]],type='l',log='x',xlab='[Hz]',ylab=paste('SPL [dB] re 20µPa, bin =',round(100/nfft,3),'Hz'),ylim=c(20,80),xlim=c(0.05,50),col='blue')
lines(psd_hann[[1]]$freq,10*log10(2*psd_hann[[1]]$spec*100/nfft),col='red')
lines(psd_rect[[1]]$freq,10*log10(2*psd_rect[[1]]$spec*100/nfft),col='green')
legend('topright',c('specgram','pwelch hanning','pwelch rectangular'),lty=1,col=c('blue','red','green'))


png("frequenzspektrum.png",width = 8, height = 5, res = 300, units = 'in')
plot(spec$f,SPLs[[1]],type='l',log='x',xlab='[Hz]',ylab=paste('SPL [dB] re 20µPa, bin =',round(100/nfft,3),'Hz'),ylim=c(20,80),xlim=c(0.05,50),col='blue')
lines(spec$f,SPLs[[2]],col='red')
lines(spec$f,SPLs[[3]],col='green')
lines(spec$f,SPLs[[4]],col='orange')
legend('topright',c('0 rpm','20 rpm','26 rpm','20 rpm +'),lty=1,col=c('green','red','blue','orange'))
dev.off()

## Calculate dB value for 2.BPH at 26rpm
s=spec$f>2.4 & spec$f<2.8
dbSum(SPLs[[1]][s])
10*log10(2*sum(psd_hann[[1]]$spec[s])*100/nfft)
10*log10(2*sum(psd_rect[[1]]$spec[s])*100/nfft)

## RESULT: main patterns can be reproduced but
## - SPLs of BPHs are around 59dB (BGR ~ 90dB)
## - The BGR curves are in total >30dB over our calculations