Sơ đồ chi tiết phân tích điều hòa thủy triều

Đến nay ở nước ta có nhiều khả năng thu thập những chuỗi quan trắc mực nước biển có độ dài nhiều chục năm tại những trạm khí tượng hải văn ven bờ và hải đảo. Từ những chuỗi mực nước đó có thể phân tích để thu được bộ các hằng số điều hòa thủy triều với nhiều phân triều hơn, đồng thời tăng độ chính xác của các hằng số điều hòa để phục vụ dự tính thủy triều tốt hơn và nhiều bài toán nghiên cứu và tính toán ứng dụng khác về mực nước và dòng chảy trong biển. Ngược lại, thực tế điều tra khảo sát tìm kiếm tại các điểm ngoài khơi và ven bờ thường cho những chuỗi quan trắc mực nước và dòng chảy ngắn một số ngày do điều kiện quan trắc khó khăn và tốn kém hoặc những gián đoạn bất thường trong công việc khảo sát trên biển. Được biết, những phương pháp phân tích truyền thống và phổ biến hiện nay thường kèm theo những quy định khắt khe về độ dài chuỗi, tính liên tục của chuỗi và độ phân giải về thời gian của quan trắc, đôi khi làm cho số liệu quan trắc trở thành vô dụng. Chúng tôi thử nghiệm xây dựng một chương trình máy tính phân tích điều hòa bằng phương pháp bình phương nhỏ nhất với một sơ đồ phân tích chi tiết, mềm dẻo, phân tích được những chuỗi quan trắc liên tục dài cỡ nhiều chục năm với mục đích thu được nhiều phân triều chính xác, tin cậy, đồng thời phân tích được những chuỗi ngắn, không liên tục về thời gian quan trắc, độ phân giải (bước gián đoạn thời gian) khác nhau. Trong mục 1 sẽ tóm tắt về bản chất lý thuyết của phương pháp, phân tích các chi tiết chứng tỏ những ưu việt của sơ đồ phân tích của chương trình. Mục 2 giới thiệu về chương trình máy tính xây dựng trên sơ đồ này và kết quả thử nghiệm phân tích để chứng minh tính hiệu quả của chương trình thông qua so sánh kết quả phân tích các chuỗi dòng chảy độ dài khác nhau, thông báo về bộ hằng số điều hòa thủy triều phân tích được cho các cảng chính của Việt Nam với chuỗi số liệu đầy đủ nhất.

Phân tích điều hòa thủy triều dựa trên những phương pháp truyền thống do các nhà hải dương học kinh điển thế giới đề xuất có tính đến đặc điểm về chu kỳ của các dao động thủy triều và tập quán quan trắc mực nước liên tục từng giờ một trong ngày tại các cảng biển (xem [1]). Các phần mềm phân tích thủy triều hiện đại trên thế giới, kể cả những phầm mềm chính thức dùng tại các trung tâm mực nước đại dương quốc tế (xem tổng quan trong [2]) hiện nay đều dựa trên phương pháp bình phương nhỏ nhất, cho phép phân tích ra bộ hằng số điều hòa đến nhiều chục phân triều tùy thuộc vào độ dài chuỗi mực nước quan trắc liên tục từng giờ trong thời kỳ một hoặc hai năm. Bài báo [3] có thể xem là một trong những thông báo sớm nhất về áp dụng phương pháp bình phương nhỏ nhất vào phân tích điều hòa thủy triều ở Việt Nam. Tuy nhiên, sơ đồ tính toán nằm trong cơ sở của tất cả các chương trình phân tích nói trên không có gì đổi mới về nguyên tắc so với các phương pháp phân tích truyền thống.

Độ cao mực nước thủy triều z size 12{z} {} tại thời điểm bất kỳ t size 12{t} {} là tổng của các dao động triều thành phần (gọi là các phân triều hay các sóng triều):

zt=A0+∑i=1rfiHicos[qit+(V0+u)i−gi] size 12{z rSub { size 8{t} } =A rSub { size 8{0} } + Sum cSub { size 8{i=1} } cSup { size 8{r} } {f rSub { size 8{i} } H rSub { size 8{i} } "cos"` [ q rSub { size 8{i} } t+ ( V rSub { size 8{0} } +u ) rSub { size 8{i} } - g rSub { size 8{i} } ] } } {}, (1)

trong đó: A0− size 12{A rSub { size 8{0} } - {}} {} độ cao mực nước trung bình, fi− size 12{f rSub { size 8{i} } - {}} {} hệ số suy biến biên độ của phân triều i size 12{i} {}, Hi− size 12{H rSub { size 8{i} } - {}} {} hằng số điều hòa biên độ của phân triều i size 12{i} {}, qi− size 12{q rSub { size 8{i} } - {}} {} tốc độ góc không đổi của phân triều i size 12{i} {}, (V0+u)i− size 12{ ( V rSub { size 8{0} } +u ) rSub { size 8{i} } - {}} {} những phần pha thiên văn của phân triều i size 12{i} {} biểu diễn các góc giờ của những tinh tú giả định tại thời điểm t size 12{t} {}, gi− size 12{g rSub { size 8{i} } - {}} {} hằng số điều hòa về pha của phân triều i size 12{i} {}, r− size 12{r - {}} {} số lượng các phân triều. fi size 12{f rSub { size 8{i} } } {} và (V0+u)i size 12{ ( V rSub { size 8{0} } +u ) rSub { size 8{i} } } {} phụ thuộc thời gian t size 12{t} {}. Khi có n size 12{n} {} độ cao mực nước quan trắc zt size 12{z rSub { size 8{t} } } {}, nhiệm vụ của phân tích thủy triều là xác định bộ gồm r size 12{r} {} cặp hằng số điều hòa không đổi H size 12{H} {} và g size 12{g} {} cho từng phân triều của trạm nghiên cứu.

Để thuận tiện áp dụng phương pháp bình phương nhỏ nhất, người ta thường biến đổi phương trình (1) thành

zt=A0+∑i=1r(Aicosqit+Bisinqit) size 12{z rSub { size 8{t} } =A rSub { size 8{0} } + Sum cSub { size 8{i=1} } cSup { size 8{r} } { ( A rSub { size 8{i} } "cos"q rSub { size 8{i} } t+B rSub { size 8{i} } "sin"q rSub { size 8{i} } t ) } } {}, (2)

trong đó

Ai=fiHicosgi−(V0+u)i size 12{A rSub { size 8{i} } =f rSub { size 8{i} } H rSub { size 8{i} } "cos" left [g rSub { size 8{i} } - ( V rSub { size 8{0} } +u ) rSub { size 8{i} } right ]} {}, Bi=fiHisingi−(V0+u)i size 12{B rSub { size 8{i} } =f rSub { size 8{i} } H rSub { size 8{i} } "sin" left [g rSub { size 8{i} } - ( V rSub { size 8{0} } +u ) rSub { size 8{i} } right ]} {}. (3)

Biết mực nước tại n size 12{n} {} giờ, người ta có n size 12{n} {} phương trình đại số dạng (2) đối với các ẩn số Ai size 12{A rSub { size 8{i} } } {} và Bi size 12{B rSub { size 8{i} } } {} để giải bằng phương pháp bình phương nhỏ nhất. Từ mỗi cặp ẩn Ai size 12{A rSub { size 8{i} } } {} và Bi size 12{B rSub { size 8{i} } } {} tìm được sẽ tính ra

Hi=Ai2+Bi2fi size 12{H rSub { size 8{i} } = { { sqrt {A rSub { size 8{i} } rSup { size 8{2} } +B rSub { size 8{i} } rSup { size 8{2} } } } over {f rSub { size 8{i} } } } } {}, gi=arctgBiAi+(V0+u)i size 12{g rSub { size 8{i} } ="arctg" { {B rSub { size 8{i} } } over {A rSub { size 8{i} } } } + ( V rSub { size 8{0} } +u ) rSub { size 8{i} } } {}. (4)

Chuỗi quan trắc càng dài, số phương trình dạng (2) càng nhiều, thì A0 size 12{A rSub { size 8{0} } } {} và số cặp hằng số điều hòa H size 12{H} {} và g size 12{g} {} nhận được càng nhiều, càng chính xác. Với một năm quan trắc có thể xác định được khoảng 60-68 cặp hằng số điều hòa H size 12{H} {} và g size 12{g} {} của điểm quan trắc.

Nhược điểm cơ bản của các phương trình dạng (2) là những đại lượng thiên văn biến thiên với thời gian f size 12{f} {} và (V0+u) size 12{ ( V rSub { size 8{0} } +u ) } {} của mỗi dao động thành phần i size 12{i} {} đã bị xem là không đổi trong thời gian quan trắc và bị đưa vào trong các ẩn số của các phương trình, do đó từng phương trình ở dạng (2) trở thành không chính xác, bởi vì trong thực tế mỗi dao động phân triều ở công thức (1) là một dao động điều biến biên độ, f size 12{f} {} biến đổi với thời gian và phần phụ pha (V0+u) size 12{ ( V rSub { size 8{0} } +u ) } {}cũng biến đổi với thời gian một cách đáng kể. Khi tính Hi size 12{H rSub { size 8{i} } } {} và gi size 12{g rSub { size 8{i} } } {} theo các công thức (4) người ta phải dùng giá trị trung bình của fi size 12{f rSub { size 8{i} } } {} tại thời điểm giữa thời kỳ quan trắc và giá trị của (V0+u)i size 12{ ( V rSub { size 8{0} } +u ) rSub { size 8{i} } } {} tại thời điểm đầu thời kỳ quan trắc. Điều này lại gây nên những mâu thuẫn kỹ thuật như: chuỗi quan trắc càng dài thì sai số càng tăng, chuỗi không liên tục (ví dụ 2 năm quan trắc không kế tiếp, mà cách xa nhau) thì không thể có thời điểm giữa quan trắc...

Các chương trình phân tích điều hòa thủy triều bằng phương pháp bình phương nhỏ nhất hiện nay xuất phát từ công thức (2) và mang những nhược điểm cơ bản như vậy.

Trong sơ đồ phân tích của chương trình do chúng tôi xây dựng phương trình độ cao mực nước triều (1) đã được biến đổi theo một kiểu khác, do Peresipkin [4] đề xuất, cho phép tính tới sự biến đổi của các đại lượng thiên văn f size 12{f} {} và (V0+u) size 12{ ( V rSub { size 8{0} } +u ) } {} với thời gian. Nếu nhóm riêng biệt các đại lượng biến thiên với thời gian và không biến thiên với thời gian bằng các ký hiệu:

ai=ficos[qit+(V0+u)i],bi=fisin[qit+(V0+u)i],Xi=Hicosgi,Yi=Hisingi,alignc { stack { size 12{a rSub { size 8{i} } =f rSub { size 8{i} } "cos" [ q rSub { size 8{i} } t+ ( V rSub { size 8{0} } +u ) rSub { size 8{i} } ] ,``````b rSub { size 8{i} } =f rSub { size 8{i} } "sin" [ q rSub { size 8{i} } t+ ( V rSub { size 8{0} } +u ) rSub { size 8{i} } ] ,} {} # X rSub { size 8{i} } =H rSub { size 8{i} } "cos"g rSub { size 8{i} } ,`````Y rSub { size 8{i} } =H rSub { size 8{i} } "sin"g rSub { size 8{i} } , {} } } {} (5)

phương trình độ cao mực nước (1) trở thành:

zt=A0+∑i=1r[(ai)tXi+(bi)tYi] size 12{z rSub { size 8{t} } =A rSub { size 8{0} } + Sum cSub { size 8{i=1} } cSup { size 8{r} } {` [ ( a rSub { size 8{i} } ) rSub { size 8{t} } X rSub { size 8{i} } + ( b rSub { size 8{i} } ) rSub { size 8{t} } Y rSub { size 8{i} } ] } } {}. (6)

Thấy rằng những đại lượng không phụ thuộc thời gian bây giờ nằm trong các ẩn số X size 12{X} {} và Y size 12{Y} {}. Còn những đại lượng phụ thuộc thời gian nằm trong các hệ số ai size 12{a rSub { size 8{i} } } {} và bi size 12{b rSub { size 8{i} } } {} của mỗi phương trình, do đó chúng được tính đến đầy đủ khi lập ra hệ n size 12{n} {} phương trình ứng với n size 12{n} {} độ cao mực nước quan trắc tại những thời điểm khác nhau. Vì vậy gọi là sơ đồ chi tiết. Giải những phương trình này bằng phương pháp bình phương nhỏ nhất, tìm được các ẩn số A0,Xi,Yi size 12{A rSub { size 8{0} } ,``X rSub { size 8{i} } ,``Y rSub { size 8{i} } } {}, từ đó tính các cặp hằng số điều hòa:

Hi=Xi2+Yi2,gi=arctgYiXi size 12{H rSub { size 8{i} } = sqrt {X rSub { size 8{i} } rSup { size 8{2} } +Y rSub { size 8{i} } rSup { size 8{2} } } ,`````g rSub { size 8{i} } ="arctg" { {Y rSub { size 8{i} } } over {X rSub { size 8{i} } } } } {}. (7)

Phương pháp bình phương nhỏ nhất cho phép xác định các ẩn số của những phương trình (6) sao cho

∑t=t1tnzt−A0−∑i=1r[(ai)tXi+(bi)tYi]2→min size 12{ Sum cSub { size 8{t=t rSub { size 6{1} } } } cSup {t rSub { size 6{n} } } { left lbrace z rSub { size 8{t} } - A rSub { size 8{0} } - Sum cSub { size 8{i=1} } cSup { size 8{r} } { [ ( a rSub { size 8{i} } ) rSub { size 8{t} } X rSub { size 8{i} } + ( b rSub { size 8{i} } ) rSub { size 8{t} } Y rSub { size 8{i} } ] } right rbrace rSup { size 8{2} } rightarrow "min"} } {}.

Điều kiện cực tiểu này sẽ cho một hệ gồm 2r+1 size 12{2r+1} {} phương trình đại số tuyến tính (hệ phương trình chính tắc), trong đó r− size 12{r - {}} {} số các phân triều được phân tích (từ phân triều M2 size 12{M rSub { size 8{2} } } {} đến phân triều cuối cùng được quy ước ký hiệu là W size 12{W} {}):

ở đây ký hiệu [ ] dùng để chỉ phép lấy tổng theo thời gian từ t1 đến tn.

Rõ ràng sơ đồ chi tiết khắc phục được những nhược điểm của các phương pháp phân tích truyền thống. Thực tế các đại lượng f size 12{f} {} và (V0+u) size 12{ ( V rSub { size 8{0} } +u ) } {}, do đó các hệ số a size 12{a} {} và b size 12{b} {}, trong sơ đồ này có thể tính chi tiết, tỉ mỉ ứng với từng thời điểm quan trắc độ cao mực nước z size 12{z} {}. Độ cao mực nước zt size 12{z rSub { size 8{t} } } {} có thể lấy tại thời điểm bất kỳ. Ta có thể ghép các độ cao mực nước quan trắc lẻ tẻ ở các năm tháng khác nhau thành một chuỗi để phân tích, do đó làm tăng số phương trình dạng (6), tăng độ chính xác phân tích. Ưu điểm này đặc biệt quan trọng đối với quan trắc dòng chảy; dòng chảy thường khó quan trắc dài ngày, nhưng được ghi với bước thời gian khác nhau, thường bé hơn một giờ, một điểm trên biển có thể có vài lần quan trắc dòng chảy vào các năm khác nhau, nếu ghép tất cả các số đo lại với nhau ta được nhiều phương trình dạng (6), tức tận dụng được thông tin.

Khi tính các hằng số điều hòa đối với những chuỗi quá ngắn, không đủ để tách những phân triều chính, thì một số phân triều có thể được xác định gần đúng dựa trên cơ sở các mối quan hệ lý thuyết giữa các phân triều có tần số (hay chu kỳ) gần bằng nhau.

Trong mỗi cặp các phân triều với tần số dao động gần nhau ( K2−S2 size 12{K rSub { size 8{2} } - S rSub { size 8{2} } } {}, P1−K1 size 12{P rSub { size 8{1} } - K rSub { size 8{1} } } {}, Q1−O1 size 12{Q rSub { size 8{1} } - O rSub { size 8{1} } } {}, N2−M2 size 12{N rSub { size 8{2} } - M rSub { size 8{2} } } {}) mà để tách được chúng đáng lẽ cần phải có chuỗi quan trắc đủ dài, ta có thể biểu diễn một phân triều (ít quan trọng hơn) theo các yếu tố của phân triều kia xuất phát từ những mối quan hệ lý thuyết giữa chúng. Như vậy, tùy thuộc vào độ dài quan trắc có thể biểu diễn được từ một đến bốn phân triều và kết quả là số ẩn trong hệ các phương trình (6) sẽ giảm đi 2, 4, 6 hoặc 8 ẩn. Khi thay thế tất cả bốn phân triều (gọi là “phương án 1”) độ dài chuỗi quan trắc theo điều kiện tách phân triều n≥360°qi−qj size 12{n >= { {"360" rSup { size 8{ circ } } } over {q rSub { size 8{i} } - q rSub { size 8{j} } } } } {} phải không ít hơn 15 ngày, khi thay thế các phân triều trong hai cặp K2−S2 size 12{K rSub { size 8{2} } - S rSub { size 8{2} } } {} và P1−K1 size 12{P rSub { size 8{1} } - K rSub { size 8{1} } } {} (“phương án 2”) - độ dài chuỗi không ít hơn 30. Trong trường hợp đầu có thể phân tích các hằng số điều hòa của 10 phân triều cơ bản ( M2 size 12{M rSub { size 8{2} } } {}, S2 size 12{S rSub { size 8{2} } } {}, K2 size 12{K rSub { size 8{2} } } {}, N2 size 12{N rSub { size 8{2} } } {}, K1 size 12{K rSub { size 8{1} } } {}, O1 size 12{O rSub { size 8{1} } } {}, P1 size 12{P rSub { size 8{1} } } {}, Q1 size 12{Q rSub { size 8{1} } } {}, M4 size 12{M rSub { size 8{4} } } {}, M6 size 12{M rSub { size 8{6} } } {}), trong trường hợp thứ hai － 11 phân triều (tính thêm được phân triều MS4 size 12{ ital "MS" rSub { size 8{4} } } {}). Trên thực tế với những chuỗi quan trắc ngắn hơn nữa vẫn nhận được những kết quả đủ thoả mãn [4].

Những quan hệ lý thuyết giữa các hằng số điều hòa của các phân triều với tần số gần nhau dựa trên những lập luận sau [4]: Tỉ số của các biên độ trung bình của các phân triều được chấp nhận bằng tỉ số của các hệ số trung bình của các phân triều đó trong khai triển chuỗi hàm thế vị lực tạo triều. Các hằng số điều hòa về pha của các phân triều tần số gần nhau chấp nhận là bằng nhau:

HK2=13,67HS2,gK2=gS2 size 12{H rSub { size 8{K rSub { size 6{2} } } } = { {1} over {3,"67"} } H rSub {S rSub { size 6{2} } } size 12{,`````g rSub {K rSub { size 6{2} } } } size 12{ {}=g rSub {S rSub { size 6{2} } } }} {}, HP1=13HK1,gP1=gK1 size 12{H rSub { size 8{P rSub { size 6{1} } } } = { {1} over {3} } H rSub {K rSub { size 6{1} } } size 12{,`````g rSub {P rSub { size 6{1} } } } size 12{ {}=g rSub {K rSub { size 6{1} } } }} {},

HQ1=15HO1,gQ1=gO1 size 12{H rSub { size 8{Q rSub { size 6{1} } } } = { {1} over {5} } H rSub {O rSub { size 6{1} } } size 12{,`````g rSub {Q rSub { size 6{1} } } } size 12{ {}=g rSub {O rSub { size 6{1} } } }} {}, HN2=15HM2,gN2=gM2 size 12{H rSub { size 8{N rSub { size 6{2} } } } = { {1} over {5} } H rSub {M rSub { size 6{2} } } size 12{,`````g rSub {N rSub { size 6{2} } } } size 12{ {}=g rSub {M rSub { size 6{2} } } }} {}. (8)

Với những quan hệ này, phương trình độ cao thủy triều dạng (6) gồm 11 phân triều có thể viết lại cụ thể như sau:

+(aMS4)tXMS4+(bMS4)tYMS4 size 12{+ ( a rSub { size 8{ ital "MS" rSub { size 6{4} } } } ) rSub {t} size 12{X rSub { ital "MS" rSub { size 6{4} } } } size 12{+ ( b rSub { ital "MS" rSub { size 6{4} } } } size 12{ ) rSub {t} } size 12{Y rSub { ital "MS" rSub { size 6{4} } } }} {} (9)

với các ký hiệu

a M 2 N 2 = f M 2 cos [ q M 2 t + ( V 0 + u ) M 2 ] + + 1 5 f N 2 cos [ q N 2 t + ( V 0 + u ) N 2 ] ; alignl { stack { size 12{a rSub { size 8{M rSub { size 6{2} } N rSub { size 6{2} } } } =f rSub {M rSub { size 6{2} } } size 12{"cos"` [ q rSub {M rSub { size 6{2} } } } size 12{t+ ( V rSub {0} } size 12{+u ) rSub {M rSub { size 6{2} } } } size 12{ ] +{}}} {} # size 12{+ { {1} over {5} } f rSub { size 8{N rSub { size 6{2} } } } "cos"` [ q rSub {N rSub { size 6{2} } } size 12{t+ ( V rSub {0} } size 12{+u ) rSub {N rSub { size 6{2} } } } size 12{ ] ;}} {} } } {} b M 2 N 2 = f M 2 sin [ q M 2 t + ( V 0 + u ) M 2 ] + + 1 5 f N 2 sin [ q N 2 t + ( V 0 + u ) N 2 ] ; alignl { stack { size 12{b rSub { size 8{M rSub { size 6{2} } N rSub { size 6{2} } } } =f rSub {M rSub { size 6{2} } } size 12{"sin"` [ q rSub {M rSub { size 6{2} } } } size 12{t+ ( V rSub {0} } size 12{+u ) rSub {M rSub { size 6{2} } } } size 12{ ] +{}}} {} # size 12{+ { {1} over {5} } f rSub { size 8{N rSub { size 6{2} } } } "sin"` [ q rSub {N rSub { size 6{2} } } size 12{t+ ( V rSub {0} } size 12{+u ) rSub {N rSub { size 6{2} } } } size 12{ ] ;}} {} } } {}

a S 2 K 2 = f S 2 cos [ q S 2 t + ( V 0 + u ) S 2 ] + + 1 3, 67 f K 2 cos [ q K 2 t + ( V 0 + u ) K 2 ] ; alignl { stack { size 12{a rSub { size 8{S rSub { size 6{2} } K rSub { size 6{2} } } } =f rSub {S rSub { size 6{2} } } size 12{"cos"` [ q rSub {S rSub { size 6{2} } } } size 12{t+ ( V rSub {0} } size 12{+u ) rSub {S rSub { size 6{2} } } } size 12{ ] +{}}} {} # size 12{+ { {1} over {3,"67"} } f rSub { size 8{K rSub { size 6{2} } } } "cos"` [ q rSub {K rSub { size 6{2} } } size 12{t+ ( V rSub {0} } size 12{+u ) rSub {K rSub { size 6{2} } } } size 12{ ] ;}} {} } } {} b S 2 K 2 = f S 2